# Create bookie plots (TpcQATask output) for files # or automatically for a whole run # Author: Felix Boehmer import argparse, ROOT, glob, string, random, subprocess, array #settings for dEdxVsMom projection PMIN = 0.300 PMAX = 0.350 #settings for VelVsMom projection PMINVEL = 0.500 PMAXVEL = 0.550 PMINVELBAR = 0.250 PMAXVELBAR = 0.300 PMINVELRPC = 0.400 PMAXVELRPC = 0.450 TEMP="/tmp/" #DEFINE SPECIAL DRAWING OPTIONS HERE drawOptions = {} padOptions = {} #MAKE SURE THE TREATMENT OF WHAT YOU WANT TO PUT IN PAD OPTIONS IS IMPLEMENTED FURTHER DOWN drawOptions["TpcdEdxVsMomHist"] = ["colz"] padOptions["TpcdEdxVsMomHist"] = ["logz"] drawOptions["CdcdEdxVsMomHist"] = ["colz"] padOptions["CdcdEdxVsMomHist"] = ["logz"] drawOptions["TpcXYSPHits"] = ["colz"] padOptions["TpcXYSPHits"] = ["logz"] drawOptions["TpcLastSPHits"] = ["colz"] padOptions["TpcLastSPHits"] = ["logz"] drawOptions["TpcCdcPhaseSpace"] = ["colz"] padOptions["TpcCdcPhaseSpace"] = ["logz"] drawOptions["RpcMomVsVelHist"] = ["colz"] padOptions["RpcMomVsVelHist"] = ["logz"] drawOptions["BarMomVsVelHist"] = ["colz"] padOptions["BarMomVsVelHist"] = ["logz"] drawOptions["TpcTotAmps2D"] = ["colz"] padOptions["TpcTotAmps2D"] = ["logz"] drawOptions["TpcPadSig2D"] = ["colz"] #padOptions["TpcPadSig2D"] = ["logz"] drawOptions["CdcHitMap"] = ["colz"] padOptions["CdcHitMap"] = ["logz"] drawOptions["CdcMomDiffVsTheta"] = ["colz"] padOptions["CdcMomDiffVsTheta"] = ["logz"] drawOptions["TpcCdcDeltaZVsZ"] = ["colz"] padOptions["TpcCdcDeltaZVsZ"] = ["logz"] drawOptions["TpcCdcDeltaThetaVsZ"] = ["colz"] padOptions["TpcCdcDeltaThetaVsZ"] = ["logz"] drawOptions["TpcCdcDeltaPhiVsZ"] = ["colz"] padOptions["TpcCdcDeltaPhiVsZ"] = ["logz"] padOptions["TpcdEdxSlice"] = ["logy"] padOptions["CdcdEdxSlice"] = ["logy"] padOptions["TpcSampleAmps"] = ["logy"] padOptions["TpcPval"] = ["logy"] padOptions["CdcPval"] = ["logy"] padOptions["TpcCdcPval"] = ["logy"] #argument parsing parser=argparse.ArgumentParser(description='') group = parser.add_mutually_exclusive_group() group.add_argument("-f","--files",type=str,help="Input ROOT files (output of TPC reconstruction)\nMutually exclusive with the -r option") group.add_argument("-r","--run",type=str,help="Run number to process. Will try to find all matching files from folder given by -F") parser.add_argument("-F","--folder",type=str,help="Folder to look in for reco files when option -r was given.") parser.add_argument("-o","--outfile",type=str,help="Output ROOT file for storing histograms") parser.add_argument("-n","--nplots",type=int,help="Plot n histograms per Canvas; Default: 6",default=6) parser.add_argument("-t","--tempdir",type=str,help="Temporary pdfs go here; Default: %s"%TEMP,default=TEMP) parser.add_argument("--pdf",help="Flag: Write pdf file (same name as ROOT output)",action="store_true") parser.add_argument("--pretty",help="FLAG: Pretty Bernhard-Style plots", action="store_true") parser.add_argument("--batch",help="FLAG: Batch mode; close after draw", action="store_true") args = parser.parse_args() if args.tempdir[-1] != '/' : TEMP = args.tempdir + '/' else : TEMP = args.tempdir ROOT.gROOT.ProcessLine(".x rootlogon.C") if args.pretty : ROOT.gROOT.ProcessLine(".x rootlogon_Bernhard_New.C") ROOT.gROOT.ProcessLine('gROOT->SetStyle("col")') ROOT.gROOT.ProcessLine('gROOT->ForceStyle()') #check if files or a run was given if args.files != None : files = glob.glob(args.files) elif args.run != None : if args.folder == None : print "ERROR: Run option given but no folder specified (use -F option) " exit() folder = args.folder if not folder.endswith("/") : folder += "/" files = glob.glob(folder+("*%s*.reco.root"%args.run)) rfiles = [] #ROOT TFile references histos = {} #output histos lhistos = [] mhistos = [] for f in sorted(files) : print f rfile = ROOT.TFile(f,"read") rfiles.append(rfile) dir = rfile.GetDirectory("TpcQA") if dir == None : print " ... something wrong with input file ... SKIPPING" continue if dir.GetName() != "TpcQA" : print " ... something wrong with input file ... SKIPPING" continue keys = dir.GetListOfKeys() if args.batch : ROOT.gROOT.SetBatch() for k in range(keys.GetSize()) : if keys.At(k).GetName() not in histos.keys() : khisto = dir.Get(keys.At(k).GetName()).Clone() classname = khisto.Class().GetName() if "TH1" in classname : khisto.SetFillColor(ROOT.kAzure-8) histos[keys.At(k).GetName()] = khisto if "TpcCdcMatching" in keys.At(k).GetName() : if "_matched" not in keys.At(k).GetName() : mhistos.append(khisto) else : lhistos.append(khisto) else : histos[keys.At(k).GetName()].Add(dir.Get(keys.At(k).GetName())) #produce sliced plot from dEdxVsMom plot: hcount = 0 for h in lhistos : #produce additional slice of dE/dx (CDC) if h.GetName() == "CdcdEdxVsMomHist" : binlow = h.GetXaxis().FindBin(PMIN) binhigh = h.GetXaxis().FindBin(PMAX) slicehist = h.ProjectionY("CdcdEdxSlice",binlow,binhigh,"e") slicehist.SetTitle("dEdx Slice %d < p < %d MeV/c (CDC)"%((int)(1.e3*PMIN),(int)(1.e3*PMAX))) slicehist.GetXaxis().SetTitle("dE/dx_{trunc} (a.u.)") slicehist.GetXaxis().SetRangeUser(0,14000) lhistos.insert(hcount+1,slicehist) #produce additional slice of dE/dx (TPC) if h.GetName() == "TpcdEdxVsMomHist" : binlow = h.GetXaxis().FindBin(PMIN) binhigh = h.GetXaxis().FindBin(PMAX) slicehist = h.ProjectionY("TpcdEdxSlice",binlow,binhigh,"e") slicehist.SetTitle("dEdx Slice %d < p < %d MeV/c (TPC)"%((int)(1.e3*PMIN),(int)(1.e3*PMAX))) slicehist.GetXaxis().SetTitle("dE/dx_{trunc} (a.u.)") slicehist.GetXaxis().SetRangeUser(0,6) lhistos.insert(hcount+1,slicehist) #produce additional slice of RPC velocity (TPC) if h.GetName() == "RpcMomVsVelHist" : binlow = h.GetXaxis().FindBin(PMINVELRPC) binhigh = h.GetXaxis().FindBin(PMAXVELRPC) slicehist = h.ProjectionY("RpcVelSlice",binlow,binhigh,"e") slicehist.SetTitle("RPC vel Slice %d < p < %d MeV/c (TPC)"%((int)(1.e3*PMINVELRPC), (int)(1.e3*PMAXVELRPC))) slicehist.GetXaxis().SetTitle("RPC velocity (cm/ns)") lhistos.insert(hcount+1,slicehist) #produce additional slice of BAR velocity (TPC) if h.GetName() == "BarMomVsVelHist" : binlow = h.GetXaxis().FindBin(PMINVELBAR) binhigh = h.GetXaxis().FindBin(PMAXVELBAR) slicehist = h.ProjectionY("BarVelSlice",binlow,binhigh,"e") slicehist.SetTitle("BAR vel Slice %d < p < %d MeV/c (TPC)"%((int)(1.e3*PMINVELBAR), (int)(1.e3*PMAXVELBAR))) slicehist.GetXaxis().SetTitle("Barrel velocity (cm/ns)") lhistos.insert(hcount+1,slicehist) #plot 2D version of pad amplitudes instead : if h.GetName() == "TpcTotAmpsPerPadHist" : print "_\_\ Building 2D version of total samples" padPlaneFile = "tpc/FOPI/par/padplane_FOPI.dat" shapeFile = "tpc/parfiles/PadShapePrototype.dat" pool = ROOT.TpcPadShapePool(shapeFile) plane = ROOT.TpcPadPlane(padPlaneFile,pool) polyHist = ROOT.TH2Poly("TpcTotAmps2D","Total Sample Amplitudes",-16,16,-16,16) polyHist.GetXaxis().SetTitle("x (cm)") polyHist.GetYaxis().SetTitle("y (cm)") for p in range(h.GetXaxis().GetNbins()) : val = h.GetBinContent(p+1) pad = plane.GetPad(p) centx = pad.x() centy = pad.y() x = [] y = [] for i in range(6) : boundVec = pad.GetBoundaryPoint(i) x.append(boundVec[0]) y.append(boundVec[1]) polyHist.AddBin(6,array.array('d',x),array.array('d',y)) polyHist.Fill(centx,centy,val) lhistos[hcount] = polyHist #plot 2D version of pad sigma instead : if h.GetName() == "TpcSampleSigPerPadHist" : print "_\_\ Building 2D version of PAD sigmas" padPlaneFile = "tpc/FOPI/par/padplane_FOPI.dat" shapeFile = "tpc/parfiles/PadShapePrototype.dat" if pool == None: pool = ROOT.TpcPadShapePool(shapeFile) if plane is None : plane = ROOT.TpcPadPlane(padPlaneFile,pool) polyHist = ROOT.TH2Poly("TpcPadSig2D","Tpc PAD Sigmas",-16,16,-16,16) polyHist.GetXaxis().SetTitle("x (cm)") polyHist.GetYaxis().SetTitle("y (cm)") for p in range(h.GetXaxis().GetNbins()) : val = h.GetBinContent(p+1) pad = plane.GetPad(p) centx = pad.x() centy = pad.y() x = [] y = [] for i in range(6) : boundVec = pad.GetBoundaryPoint(i) x.append(boundVec[0]) y.append(boundVec[1]) polyHist.AddBin(6,array.array('d',x),array.array('d',y)) polyHist.Fill(centx,centy,val) lhistos[hcount] = polyHist hcount += 1 offp = len(lhistos) % args.nplots if offp > 0 : offp = 1 npages = len(lhistos)/args.nplots + offp canvases = [] if args.outfile != None : outfile = ROOT.TFile(args.outfile,"recreate") #lhistos = histos.values() for ic in range(npages) : c = ROOT.TCanvas("c%d"%(ic),"c%d"%(ic),800,1000) canvases.append(c) offset = args.nplots % 2 if offset > 0 : offset = 1 c.Divide(2,args.nplots/2 + offset) for h in range(args.nplots) : pad = c.cd(h+1) if h+ic*args.nplots >= len(lhistos) : break ihisto = lhistos[h+ic*args.nplots] gTitle = ihisto.GetTitle() gName = ihisto.GetName() print gName #SPECIAL: resize lashit z-range if gName == "TpcLastSPHits" : print " ... resizing TpcLastSPHits ... " zhist = ihisto.ProjectionX("testz") zax = zhist.GetXaxis() for i in range(zax.GetNbins()) : if zhist.GetBinContent(i+1) > 0. : minz = zax.GetBinCenter(i+1) break for i in range(zax.GetNbins(),0,-1) : if zhist.GetBinContent(i) > 0. : maxz = zax.GetBinCenter(i) break ihisto.GetXaxis().SetRangeUser(minz*1.2,maxz*1.2) print " ... done. " #SPECIAL: resize z-range if gName == "TpcZSPHits" : print " ... resizing TpcZSPHits ... " zax = ihisto.GetXaxis() for i in range(zax.GetNbins()) : if ihisto.GetBinContent(i+1) > 0. : minz = zax.GetBinCenter(i+1) break for i in range(zax.GetNbins(),0,-1) : if ihisto.GetBinContent(i) > 0. : maxz = zax.GetBinCenter(i) break ihisto.GetXaxis().SetRangeUser(minz*1.2,maxz*1.2) print " ... done. " if args.run != None : gTitle += " [RUN %s]"%args.run ihisto.SetTitle(gTitle) ihisto.Draw() if gName in drawOptions.keys() : for do in drawOptions[gName] : ihisto.SetDrawOption(do) print " ... Applying draw option %s" % do if gName in padOptions.keys() : for po in padOptions[gName] : if po == "logz" : pad.SetLogz() if po == "logy" : pad.SetLogy() if args.outfile != None : outfile.cd() ihisto.Write() offc = npages npages = len(mhistos)/args.nplots + offp lines = [] for ic in range(npages) : c = ROOT.TCanvas("c%d"%(ic+offc),"c%d"%(ic+offc),800,1000) canvases.append(c) offset = args.nplots % 2 if offset > 0 : offset = 1 c.Divide(2,args.nplots/2 + offset) for h in range(args.nplots) : pad = c.cd(h+1) if h+ic*args.nplots >= len(mhistos) : break ihisto = mhistos[h+ic*args.nplots] gTitle = ihisto.GetTitle() gName = ihisto.GetName() if "_matched" in gName : continue print gName if args.run != None : gTitle += " [RUN %s]"%args.run ihisto.SetTitle(gTitle) ihisto.Draw() imax = ihisto.GetMaximum() #draw lines if cuts were applied if "cutl" in gName : il = gName.find("cutl") parts = gName[il:].split("_") low = float(parts[1]) high = float(parts[3]) ll = ROOT.TLine(low,0,low,imax) lh = ROOT.TLine(high,0,high,imax) ll.SetLineColor(ROOT.kRed+1) lh.SetLineColor(ROOT.kRed+1) ll.Draw() lh.Draw() lines.append(ll) lines.append(lh) if gName in drawOptions.keys() : for do in drawOptions[gName] : ihisto.SetDrawOption(do) print " ... Applying draw option %s" % do if gName in padOptions.keys() : for po in padOptions[gName] : if po == "logz" : pad.SetLogz() if po == "logy" : pad.SetLogy() #get matched version of this plot mname = gName+"_matched" mh = None if mname in histos.keys() : mh = histos[mname] if mh == None : print "Could not find matched version of",gName else : mh.SetFillColor(ROOT.kBlack) mh.Draw("same") if args.outfile != None : outfile.cd() ihisto.Write() if mh != None : mh.Write() #write canvases to one .pdf if args.pdf : #generate random name for temp files pdfout = args.outfile.replace(".root",".pdf") rlst = [random.choice(string.ascii_letters + string.digits) for n in range(20)] rstr = "".join(rlst) template = TEMP+"%s_c%d.pdf" ccount = 0 fnames = [] for c in canvases : fname = template % (rstr,ccount) fnames.append(fname) c.SaveAs(fname) ccount+=1 #merge them fargs = "" for fn in fnames: fargs+=(fn+" ") #command = "pdftk %s cat output %s" % (fargs,pdfout) command = "gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=%s %s" %(pdfout,fargs) cseq = command.split() #print command #print cseq rcode = subprocess.call(cseq) if rcode == 0 : print "Written .pdf to",pdfout else : print "Could not merge .pdfs .... pdftk installed?" if not args.batch : raw_input()