Newer
Older
TB_Chris / TbUT / scripts / .svn / text-base / fastAnalysis.py.svn-base
  1. #!/usr/bin/env python
  2.  
  3. """
  4. fastAnalysis.py for use during UT testbeam data taking
  5.  
  6. To use:
  7. 1. Set up TbUT running environment (https://twiki.cern.ch/twiki/bin/view/LHCb/TbUT)
  8. 2. Edit this script to point to your current eos mount point and the data dir for current testbeam, or use the command line options
  9. 3. From Tb/TbUT directory, run scripts/fastAnalysis.py with the following required options:
  10. -b board : Letter+number combination, e.g. 'A8'
  11. -p pednum : run number of the pedestal run to use (UT run number)
  12. -r runnum : run number to analyze (UT run number) -- enter 0 for pedestal analysis only
  13. 4. These are optional arguments:
  14. -n nevmax : number of events to analyze
  15. -f force : Force overwrite of any output files existing
  16. -t senstype : sensor type ('PType' or 'NType')
  17. -e eosmount : eos mount point for data to use (uses as a simple path)
  18. -i indir : input path relative to eosmount point to find data
  19. -o outdir : Directory to put the monitoring in
  20. 5. Inside outdir it will make monitoring directories, and save the pedestal .dat, root files, and defined plots there
  21. 6. It will not re-run the analysis if the files exist already unless you use -f
  22. 7. After running it will pop open a window to view the new plots, and exit when that window is closed
  23. """
  24.  
  25. ############################################################
  26. ##Modify this to change defaults:
  27. nevdef=100000
  28. eosmount='/afs/cern.ch/work/m/mrudolph/private/nov_testbeam/eos'
  29. indir='lhcb/testbeam/ut/OfficialData/October2015'
  30. outdir='.'
  31. ############################################################
  32.  
  33.  
  34.  
  35.  
  36.  
  37.  
  38. ############################################################
  39. ##Define all the draw commands for simple monitoring plots from Cluster tuple here
  40. ##Noise mean and width are automatically created as the first 2 plots
  41. ##The last plot will be the ADC in the fiducial region (that gets defined based on previous draw command postprocessing)
  42. #Potential TODO: create a new module file with the plot configuration to make extension a bit easier
  43.  
  44. class draw_command():
  45. """draw command skeleton for a TTree
  46. Need at least command, cut, option, title for histogram to make
  47. Can also bind postprocessor functions to draw more stuff or
  48. to return to the list of fiducial cuts for the final plot
  49. """
  50. def __init__(self, command, cut, option, title, postprocessor=None):
  51. self.command = command
  52. self.cut = cut
  53. self.option = option
  54. self.title = title
  55. self.postprocessor = postprocessor
  56.  
  57.  
  58. ##Post processing functions for histograms defined in draw commands -- expect TPad, TTree, TH1 parameters
  59. ## Code below could be used for a couple things -- you have the pad so you can process the hist and draw more things,
  60. ## or you can return a cut string that will be used to define the fiducial region later on
  61.  
  62. def post_nostat( pad, tree, hist ):
  63. """Generic postprocessor to remove stat box"""
  64. hist.SetStats(False)
  65.  
  66. def post_clusterpos( pad, tree, hist ):
  67. """Process the cluster position histogram and find beam area"""
  68. post_nostat( pad, tree, hist)
  69.  
  70. maxbin = hist.GetMaximumBin()
  71. maxval = hist.GetBinContent(maxbin)
  72. thresh = 0.25
  73. #look to the left and then to the right
  74. currbin = maxbin
  75. while( hist.GetBinContent(currbin) > thresh*maxval ):
  76. currbin -= 1
  77. left = hist.GetXaxis().GetBinLowEdge(currbin)
  78. currbin = maxbin
  79. while( hist.GetBinContent(currbin) > thresh*maxval ):
  80. currbin += 1
  81. right = hist.GetXaxis().GetBinUpEdge(currbin)
  82. hist.GetYaxis().SetRangeUser(0,maxval*1.1)
  83. from ROOT import TLine
  84. l1 = TLine( left ,0, left,maxval*1.1)
  85. pad.posl1=l1
  86. l1.Draw()
  87. l2 = TLine( right,0,right,maxval*1.1)
  88. pad.posl2=l2
  89. l2.Draw()
  90. pad.Update()
  91. return "clustersPosition > {} && clustersPosition < {}".format(left,right)
  92.  
  93. def post_tdc_prof( pad, tree, hist ):
  94. """Process the charge v. TDC profile and find timing cut"""
  95. ymin = hist.GetMinimum()
  96. if ymin < 0:
  97. ymin *= 1.2
  98. else:
  99. ymin *= 0.8
  100. ymax = hist.GetMaximum()
  101. if ymax < 0:
  102. ymax *= 0.8
  103. else:
  104. ymax *= 1.2
  105. hist.GetYaxis().SetRangeUser(ymin,ymax)
  106. isneg = (ymax < 0) or ( ymin < 0 and abs(ymin) > abs(ymax) )
  107. if isneg:
  108. b = hist.GetMinimumBin()
  109. else:
  110. b = hist.GetMaximumBin()
  111.  
  112. from ROOT import TLine
  113. l1 = TLine( b - 1.5 ,ymin, b - 1.5,ymax)
  114. pad.tdcl1=l1
  115. l1.Draw()
  116. l2 = TLine( b+1.5,ymin,b+1.5,ymax)
  117. pad.tdcl2=l2
  118. l2.Draw()
  119. return "abs(clustersTDC - {}) <= 1.5".format(b)
  120.  
  121. ############################################################
  122. ##The list of draw commands to run
  123. ##Make sure for each command if you specify anything about the histogram to draw to, give it a unique name!
  124. ## Fiducial adc name is "fidadc", noise are called "meanNoise" and "widthNoise"
  125. drawcommands = [ draw_command("clusterNumberPerEvent","","","Number of clusters;N_{clusters};Events"),
  126. draw_command("clustersCharge>>charge(200,-1200,0)","clustersSize > 0","","Cluster charge;Tot. charge [ADC];Clusters"),
  127. draw_command("clustersPosition","clustersSize > 0","","Cluster position;Channel;Clusters",postprocessor=post_clusterpos),
  128. draw_command("clustersSize","clustersSize > 0", "","Cluster size;N_{strips};Clusters"),
  129. draw_command("clustersTDC>>tdc(10,0.5,10.5)","","","TDC;TDC;Events",postprocessor=post_nostat),
  130. draw_command("clustersCharge:clustersTDC>>tdcp(10,0.5,10.5)","clustersSize > 0","prof","Charge v TDC;TDC;<charge> [ADC]",postprocessor=post_tdc_prof),
  131. ]
  132.  
  133. ##Should be possible to hook into an external file with TTree::MakeProxy using these commands as well if needed, though it hasn't been tested.
  134. ##If you have other histograms to calculate from available inputs, there are two lists (hlist_pre and hlist_post) in the plotting loop they can be added to;
  135. ## see where the noise histograms are created. Pre goes before the ttree command and post afterwards
  136. ############################################################
  137.  
  138.  
  139. import argparse
  140. parser = argparse.ArgumentParser(
  141. description = "Fast analysis of UT testbeam data for quality monitoring",
  142. formatter_class=argparse.ArgumentDefaultsHelpFormatter,)
  143. parser.add_argument('-b','--board',type=str,required=True,
  144. help="Letter+number combination, e.g. 'A8'")
  145. parser.add_argument('-p','--pednum',type=int,required=True,
  146. help="Run number of the pedestal run to use (UT run number)")
  147. parser.add_argument('-r','--runnum',type=int,required=True,
  148. help="Run number to analyze (UT run number); use 0 to analyze pedestal only")
  149. parser.add_argument('-n','--nevmax',type=int,required=False,default=nevdef,
  150. help="Max number of events to analyze")
  151. parser.add_argument('-f','--force',required=False,action='store_true',
  152. help="Force overwrite of any output files existing")
  153. parser.add_argument('-t','--senstype',type=str,required=False,default='PType',
  154. help="Sensor type ('PType' or 'NType')")
  155. parser.add_argument('-e','--eosmount',type=str,required=False,default=eosmount,
  156. help="Eos mount point")
  157. parser.add_argument('-i','--indir',type=str,required=False,default=indir,
  158. help='Input data directory path relative to eosmount')
  159. parser.add_argument('-o','--outdir', type=str, required=False,default=outdir,
  160. help="Output directory for root files and plots")
  161. args = parser.parse_args()
  162.  
  163. import sys
  164. import os
  165. import subprocess
  166.  
  167. #make sure I have places for output
  168. pedestaldir = args.outdir + '/Pedestal-Board{}-{}'.format( args.board, args.pednum )
  169. monitordir = args.outdir + '/Monitoring-Board{}-{}'.format( args.board, args.runnum )
  170. def ensuredir(dir):
  171. try:
  172. os.makedirs(dir)
  173. except OSError:
  174. if not os.path.isdir(dir):
  175. raise
  176. ensuredir(pedestaldir)
  177. if(args.runnum!=0):
  178. ensuredir(monitordir)
  179.  
  180. #Run pedestal if it does not already exist
  181. didPedestal = False
  182. if(args.force or not os.path.isfile("{}/Fast-Pedestal-Board{}-{}.dat".format( pedestaldir, args.board, args.pednum ) ) ):
  183. print "================================================================"
  184. print "==ANALYZE PEDESTAL FOR BOARD {} RUN {}=========================".format( args.board, args.pednum )
  185. print "================================================================"
  186.  
  187. didPedestal = True
  188.  
  189. ##Set up the gaudi config for running the pedestal
  190. pedCode = ("import sys\n"
  191. "sys.path.append( 'options/python' )\n"
  192. "from TbUTPedestalRunner import TbUTPedestalRunner\n"
  193. "app=TbUTPedestalRunner()\n"
  194. "# set parameter\n"
  195. "app.inputData= '{}/{}/Board{}/RawData/Pedestal-B{}-{}-{}-.dat'\n"
  196. "app.isAType={}\n"
  197. "# have to be more than 4k (~10k)\n"
  198. "app.eventMax={}\n"
  199. "# keep the pedestals files in $KEPLERROOT/../TbUT/options/UT/ directory !!!!!\n"
  200. "app.pedestalOutputData ='{}/Fast-Pedestal-Board{}-{}.dat'\n"
  201. "app.runPedestals()\n").format(args.eosmount,args.indir, args.board, args.board[1:], args.board[:1],args.pednum, (args.board[:1]=='A'), args.nevmax, pedestaldir, args.board,args.pednum)
  202.  
  203.  
  204. with open('myTempPedRun.py','w') as ftarget:
  205. ftarget.write(pedCode)
  206.  
  207. ret = subprocess.call(['gaudirun.py','myTempPedRun.py'])
  208.  
  209. os.remove('myTempPedRun.py')
  210.  
  211. if(ret!=0):
  212. sys.exit("Bad return code from pedestal run")
  213.  
  214. #Determine if we have a data run to process
  215. if( args.runnum != 0):
  216.  
  217. #find the data file by run number
  218. dir_to_search = "{}/{}/Board{}/RawData/".format(args.eosmount,args.indir,args.board)
  219. paths = subprocess.check_output("find {} -iname '*-{}-*.dat'".format(dir_to_search,args.runnum), shell=True).splitlines()
  220.  
  221. if( len(paths)==0):
  222. sys.exit("ERROR: no data file found for run number {}".format(args.runnum))
  223. elif( len(paths) > 1):
  224. print "WARNING: more than one file matching run number, using",paths[0]
  225.  
  226. inpath = paths[0]
  227. #These are the output names used by TbUTClusterizator
  228. outnames = [ os.path.basename(inpath).replace('.dat','.root'), os.path.basename(inpath).replace('.dat','_Tuple.root') ]
  229.  
  230. #Skip running if files exist unless forced
  231. if( args.force or not (os.path.isfile(monitordir+'/'+outnames[0]) and os.path.isfile(monitordir+'/'+outnames[1])) ):
  232.  
  233. print "================================================================"
  234. print "==ANALYZE DATA FOR BOARD {} RUN {}=============================".format( args.board, args.runnum )
  235. print "================================================================"
  236.  
  237. runCode = ("import sys\n"
  238. "sys.path.append('options/python')\n"
  239. "from TbUTClusterizator import TbUTClusterizator\n"
  240. "app = TbUTClusterizator()\n"
  241. "# set parameters\n"
  242. "app.inputData = '{}'\n"
  243. "app.isAType = {}\n"
  244. "app.sensorType = '{}'\n"
  245. "app.eventMax = {}\n"
  246. "app.pedestalInputData = '{}/Fast-Pedestal-Board{}-{}.dat'\n"
  247. "app.eventNumberDisplay = 1000\n"
  248. "app.runClusterization()\n").format(inpath,(args.board[:1]=='A'),args.senstype,args.nevmax,pedestaldir,args.board,args.pednum)
  249.  
  250.  
  251.  
  252. with open('myTempRun.py','w') as ftarget:
  253. ftarget.write(runCode)
  254.  
  255. #run twice because of noise input issue...
  256. ret = subprocess.call(['gaudirun.py','myTempRun.py'])
  257. ret = subprocess.call(['gaudirun.py','myTempRun.py'])
  258.  
  259. #Move them since currently can't control output location from the run
  260. if( monitordir != '.' ):
  261. subprocess.call(['mv',outnames[0],monitordir+'/'])
  262. subprocess.call(['mv',outnames[1],monitordir+'/'])
  263.  
  264. os.remove('myTempRun.py')
  265.  
  266. if(ret!=0):
  267. sys.exit("Bad return code from analysis run")
  268.  
  269.  
  270.  
  271.  
  272. ########################################################
  273. ##Analysis part
  274. ########################################################
  275.  
  276.  
  277. print "================================================================"
  278. print "==PLOT DATA FOR BOARD {} RUN {}================================".format( args.board, args.runnum )
  279. print "================================================================"
  280.  
  281. #Define what extensions to save plots with
  282. def saveall( c, outname ):
  283. for ext in '.png', '.pdf', '.C':
  284. c.SaveAs(outname+ext)
  285.  
  286. #setup ROOT
  287. from ROOT import gROOT, TH1F, TCanvas, TFile, TTree, gStyle, TF1,gDirectory,gPad
  288. gROOT.SetBatch()
  289. gROOT.SetStyle("Plain")
  290.  
  291. #keep track of plots we made to open at the end
  292. plotlist = []
  293.  
  294. if( didPedestal ):
  295. c = TCanvas()
  296. #analyze the pedestal only
  297. #dat file is space separated list of pedestals by channel number
  298. peds = open("{}/Fast-Pedestal-Board{}-{}.dat".format(pedestaldir,args.board,args.pednum), "r")
  299. try:
  300. pedvals = [float(p) for p in peds.read().split()]
  301. except:
  302. sys.exit("Bad pedestal data format")
  303. nc = len(pedvals)
  304. hped = TH1F("hped","Pedestal v. channel;Channel;Pedestal [ADC]",nc,-0.5, nc-0.5)
  305. for i,v in enumerate(pedvals):
  306. hped.SetBinContent(i+1, v)
  307.  
  308. hped.Draw();
  309. saveall(c,"{}/pedestal".format(pedestaldir))
  310. plotlist += ["{}/pedestal.png".format(pedestaldir)]
  311.  
  312. #analyze the run
  313. if( args.runnum != 0 ):
  314. fhists = TFile("{}/{}".format(monitordir,outnames[0]))
  315. ftrees = TFile("{}/{}".format(monitordir,outnames[1]))
  316.  
  317. t = ftrees.Get("TbUT/Clusters")
  318. #make the rough noise plots -- use the range -120 to 120 ADC and ignore most of the signal contribution
  319. noise2d = fhists.Get("TbUT/CMSData_vs_channel")
  320. noise2d.GetYaxis().SetRangeUser(-120,120)
  321. nc = noise2d.GetNbinsX()
  322.  
  323. meanNoise = TH1F("meanNoise","Mean noise;Channel;<noise> [ADC]",nc,-0.5,nc-0.5)
  324. widthNoise = TH1F("widthNoise","Width noise;Channel;#sigma_{noise} [ADC]",nc,-0.5,nc-0.5)
  325. meanNoise.SetStats(False)
  326. widthNoise.SetStats(False)
  327.  
  328. for b in range(1, nc+1):
  329. noise2d.GetXaxis().SetRange(b,b)
  330. meanNoise.SetBinContent(b, noise2d.GetMean(2))
  331. meanNoise.SetBinError(b, noise2d.GetMeanError(2))
  332. widthNoise.SetBinContent(b, noise2d.GetStdDev(2))
  333. widthNoise.SetBinError(b, noise2d.GetStdDevError(2))
  334.  
  335. #hlist_pre is a list of histograms to draw
  336. hlist_pre = [meanNoise,widthNoise]
  337. npre = len(hlist_pre)
  338.  
  339. #hlist_post, however, is a list of functions that create histograms; define some functions for post-draw histograms first:
  340. def fiducialADC( t, fidcut ):
  341. print "Use fiducial cut:", fidcut
  342. t.Draw("clustersCharge>>fidcharge(200,-1200,0)", fidcut)
  343. h = t.GetHistogram()
  344. h.SetTitle("Fiducial cluster charge;Tot. charge [ADC];Clusters")
  345. return h
  346.  
  347. #the output ntuple has empty clusters saved in every event, so always cut those out
  348. fidcut = "clustersSize > 0"
  349. hlist_post = [fiducialADC]
  350. npost = len(hlist_post)
  351.  
  352. #create the canvases
  353. #keep track of how many plots left to do
  354. ntreedraw = len(drawcommands)
  355. nplotall = npre + ntreedraw + npost
  356.  
  357. #9 plots per canvas max
  358. ncanvas = nplotall/9
  359. #if not an exact divisor we need an extra one for the leftovers:
  360. if( nplotall % 9 != 0 ):
  361. ncanvas += 1
  362. #how many columns and rows we want for each number of plots per canvas
  363. colnums = [1,2,2,2,3,3,3,3,3]
  364. rownums = [1,1,2,2,2,2,3,3,3]
  365. clist = [TCanvas('c{}'.format(i)) for i in range(ncanvas)]
  366. plotidx = 0 #count the plots
  367. for c in clist:
  368. #how many plots on this canvas?
  369. #we'll use nplotall as a countdown
  370. if( nplotall >= 9):
  371. nplot=9
  372. nplotall -= 9
  373. else:
  374. nplot=nplotall
  375. #setup the canvas based on number of plots to draw
  376. ncols = colnums[nplot - 1]
  377. nrows = rownums[nplot - 1]
  378. c.SetCanvasSize( 600*ncols, 300*nrows )
  379. c.Divide(ncols,nrows)
  380.  
  381. for i in range(nplot):
  382. #plotidx is overall iterator, i is current canvas iterator
  383. if (plotidx < npre):
  384. c.cd(i+1)
  385. hlist_pre[plotidx].Draw()
  386. elif(plotidx < npre + ntreedraw):
  387. c.cd(i+1)
  388. t.Draw( drawcommands[plotidx-npre].command, drawcommands[plotidx-npre].cut, drawcommands[plotidx-npre].option )
  389. h = t.GetHistogram()
  390. gPad.h = h
  391. h.SetTitle(drawcommands[plotidx-npre].title)
  392. if(drawcommands[plotidx-npre].postprocessor):
  393. newcut = drawcommands[plotidx-npre].postprocessor( gPad, t, h )
  394. if type(newcut) == str:
  395. fidcut += " && " + newcut
  396. else:
  397. c.cd(i+1)
  398. hlist_post[plotidx - npre - ntreedraw](t,fidcut).Draw()
  399.  
  400. plotidx+=1
  401. plotpath = "{}/monitoring-{}".format( monitordir, c.GetName()[1:] )
  402. saveall(c, plotpath)
  403. plotlist += ["{}.png".format(plotpath)]
  404.  
  405. #open up an eye of gnome window to look at the plots; use the next and previous arrows to navigate between them
  406. if( plotlist ):
  407. subprocess.call(["eog",]+plotlist)
  408.