import ROOT, glob, math, sys, os from ROOT import std # {{{ helper functions # {{{ draw prelim def drawPrelim() : prelim=ROOT.TLatex() prelim.SetTextColor(ROOT.kGray) prelim.SetTextSize(0.07) prelim.SetTextAngle(20) prelim.SetNDC() prelim.DrawLatex(0.38,0.45,"preliminary") return # }}} # {{{ draw text def drawtext(data): text=ROOT.TLatex() text.SetTextColor(ROOT.kBlack) text.SetTextSize(0.05) text.SetNDC() y=0.8 ystep=0.04 for tex in data: print tex text.DrawLatex(0.15,y,tex) y-=ystep return # }}} # {{{ draw residual data def drawresdata(data): text=ROOT.TLatex() text.SetTextColor(ROOT.kBlack) text.SetTextSize(0.05) text.SetNDC() y=0.8 ystep=0.05 text.DrawLatex(0.15,y,"#sigma_{narrow}= "+str(round(data['nsig'],1))+" #mum") y-=ystep text.DrawLatex(0.15,y,"#sigma_{broad}="+str(round(data['bsig'],1))+" #mum") y-=ystep text.DrawLatex(0.15,y,"#sigma_{mean}= "+str(round(data['csig'],1))+" #mum") y-=ystep text.DrawLatex(0.15,y,"#chi^{2}/ndf="+str(round(data['chindf'],1))) y-=ystep if (data['bsig']*data['bconst'])>0: text.DrawLatex(0.15,y,"S/B="+str(round((data['nsig']*data['nconst'])/(data['bsig']*data['bconst']),2))) return # }}} # {{{ calc Mean er 2 gauss version def calcMeanErr2(s1,ds1,h1,dh1,s2,ds2,h2,dh2) : divider=(h1*s1+h2*s2) divider *= divider ddh1 = math.fabs(h2*s1*s1*s2-h2*s2*s2*s1) ddh1 /= divider ddh1 *= dh1 # ddh1 *= ddh1 dds1 = math.fabs(h1*h1*s1*s1+2*h1*h2*s1*s2-h2*h1*s2*s2) dds1 /=divider dds1 *= ds1 # dds1 *= dds1 ddh2 = math.fabs(h1*s1*s2*s2-h1*s1*s1*s2*s2) ddh2 /= divider ddh2 *= dh2 # ddh2 *= ddh2 dds2 = math.fabs(h2*h2*s2*s2+2*h2*h1*s1*s2 - h1*h2*s1*s1) dds2 /= divider dds2 *= ds2 # dds2 *= dds2 meanerr = math.sqrt(ddh1+ddh2+dds1+dds2) # meanerr = (ddh1+ddh2+dds1+dds2) return meanerr # }}} # {{{ calc mean err 3 gauus version def ddh1s1(s1,s2,s3,h1,h2,h3,dh1,ds1): divider = h1*s1+h2*s2+h3*s3 divider *= divider err1 = (h1*s1+h2*s2+h3*s3)*s1*s1 - (h1*s1*s1+h2*s2*s2+h3*s3*s3)*s1 err1 /= divider err1 *= dh1 err1 *= err1 err2 = (h1-s1+h2-s2+h3-s3)*2*h1*s1-(h1*s1*s1+h2*s2*s2+h3*s3*s3)*h1 err2 /= divider err2 *= ds1 err2 *= err2 print "error from h:",err1," error from s",err2," tot:",err1+err2 return err1+err2 def calcMeanErr3(s1,ds1,h1,dh1,s2,ds2,h2,dh2,s3,ds3,h3,dh3): meanerr = 0 print "smallest" meanerr += ddh1s1(s1,s2,s3,h1,h2,h3,dh1,ds1) print "broadest" meanerr += ddh1s1(s1,s2,s3,h1,h2,h3,dh2,ds2) print "middle" meanerr += ddh1s1(s1,s2,s3,h1,h2,h3,dh3,ds3) print "meanerr:", meanerr meanerr = math.sqrt(meanerr) return meanerr # }}} # {{{ Fitf Residual function fitres(histo)->sig,meansig def fitres(histo,ngaus,df,factors=[1.,1.,2.,12.]) : # histo.SetStats(ROOT.kFALSE) data={} print '***************i=',i,colors[i],'*********************' #histo.SetLineColor(colors[i]) histo.SetFillColor(ROOT.kAzure-8) testfit = ROOT.TF1("testfitY"+str(1),"gaus",-10000.,10000.) testfit.SetLineColor(2) histo.Fit(testfit, "QN+", "", -10000,10000) if ngaus == 3 : rfit = ROOT.TF1("fitfuncY"+str(1),"gaus + gaus(3) + gaus(6)",-10000.,10000.) elif ngaus == 2: rfit = ROOT.TF1("fitfuncY"+str(1),"gaus + gaus(3)",-10000.,10000.) rfit1 = ROOT.TF1("gaus1","gaus",-10000.,10000.) rfit2 = ROOT.TF1("gaus2","gaus",-10000.,10000.) rfit.SetNpx(1000) rfit.SetParName(0,"const1") rfit.SetParName(1,"mean1") rfit.SetParName(2,"sigma1") rfit.SetParName(3,"const2") rfit.SetParName(4,"mean2") rfit.SetParName(5,"sigma2") if ngaus == 3: rfit.SetParName(6,"const3") rfit.SetParName(7,"mean3") rfit.SetParName(8,"sigma3") tconst = [None]*3 tmean = [None]*3 tsigma = [None]*3 tconst[0] = tconst[1] = tconst[2] = testfit.GetParameter(0) tmean[0] = tmean[1] = tmean[2] = 0 tsigma[0] = tsigma[1] = tsigma[2] = testfit.GetParameter(2) #medium gaus rfit.SetParameter(0,tconst[0]/1.5) rfit.SetParameter(1,0.) # rfit.SetParLimits(1,0.,0.) # rfit.FixParameter(1,0.) rfit.SetParameter(2,tsigma[0]*2.5) if ngaus == 3 : rfit.SetParLimits(2,tsigma[0]*0.5, tsigma[0]*2) if ngaus == 2 : rfit.SetParLimits(2,tsigma[0]*factors[2], tsigma[0]*factors[3]) #broadest gaus # if ng=2 small gaus rfit.SetParameter(3, tconst[1]) rfit.SetParameter(4, 0.) # rfit.FixParameter(4,0.) rfit.SetParameter(5, tsigma[1]*0.5) if ngaus == 3 : rfit.SetParLimits(5, tsigma[1]*0.8, tsigma[1]*2) if ngaus == 2 : rfit.SetParLimits(5, 0.0*factors[0], tsigma[1]*factors[1]) if ngaus == 3: #smallest gaus rfit.SetParameter(6, tconst[2]*4) rfit.SetParameter(7,0.) rfit.FixParameter(7,0.) rfit.SetParameter(8, tsigma[2]*0.2) rfit.SetParLimits(8, 50,tsigma[2]*0.6) rfit.SetLineWidth(1) print "now fitting in ranges:",-3*histo.GetRMS(),3*histo.GetRMS() histo.Fit(rfit, "MIQ+","", -3*histo.GetRMS(),3*histo.GetRMS()) # histo.Fit(rfit, "+","", -10000,10000) # histo.Fit(rfit, "+","", -100,100) # {{{ second fit if rfit.GetNDF()>0: scaleu=rfit.GetChisquare()/rfit.GetNDF() scaled=1/scaleu else: scaleu=2. scaled=0.5 if df==1: for j in range(ngaus) : tconst[j]=rfit.GetParameter(j*ngaus) tmean[j]=rfit.GetParameter(j*ngaus+1) tsigma[j]=rfit.GetParameter(j*ngaus+2) #medium gaus rfit.SetParameter(0,tconst[0]) # rfit.SetParLimits(0,tconst[0]*0.5, tconst[0]*2) rfit.SetParameter(2,tsigma[0]) rfit.SetParLimits(2,tsigma[0]*scaled, tsigma[0]*scaleu) #broadest gaus rfit.SetParameter(3, tconst[1]) # rfit.SetParLimits(3, 0, tconst[1]/1.5) rfit.SetParameter(5, tsigma[1]) rfit.SetParLimits(5, tsigma[1]*scaled, tsigma[1]*scaleu) #smallest gaus if ngaus==3: rfit.SetParameter(6, tconst[2]) # rfit.SetParLimits(6, tconst[2]*2, tconst[2]*5) rfit.SetParameter(8, tsigma[2]) rfit.SetParLimits(8, tsigma[2]*scaled, tsigma[2]*scaleu) rfit.SetLineWidth(2) histo.Fit(rfit, "+", "", -10000,10000) # }}} # {{{ debug fit plotting if 0: # canv=ROOT.TCanvas("bla","bla") # leg=ROOT.TLegend(0.7,.7,0.9,.9,"","NDC") # canv.cd() # histo.Draw() rfit1.SetParameter(0,rfit.GetParameter(0)) rfit1.SetParameter(1,rfit.GetParameter(1)) rfit1.SetParameter(2,rfit.GetParameter(2)) rfit1.SetLineColor(ROOT.kRed) # leg.AddEntry(rfit1,"Broad") rfit1.DrawCopy("same") rfit2.SetParameter(0,rfit.GetParameter(3)) rfit2.SetParameter(1,rfit.GetParameter(4)) rfit2.SetParameter(2,rfit.GetParameter(5)) rfit2.SetLineColor(ROOT.kBlue) # leg.AddEntry(rfit2,"Narrow") rfit2.DrawCopy("same") # testfit.SetLineColor(ROOT.kGreen) # testfit.Draw("same") # leg.AddEntry(testfit,"testfit") # leg.Draw() # print "*******************************************" # print tsigma[0]*10.2, tsigma[1]*0.05, tsigma[1]*1.5 # print "*******************************************" # drawtext(["broad"+str(rfit.GetParameter(2)),"small"+str(rfit.GetParameter(5))]) # nix=raw_input("blahorst") # }}} if ngaus == 3: int1 = rfit.GetParameter(6)*rfit.GetParameter(8) # smallest gaus elif ngaus ==2: int1 = rfit.GetParameter(0)*rfit.GetParameter(2) #medium gaus int2 = rfit.GetParameter(3)*rfit.GetParameter(5) #broadest gaus int3 = rfit.GetParameter(0)*rfit.GetParameter(2) #medium gaus meanerr = 0 if int1 > 0 : ratio = int2/int1 #mean resolution: if ngaus == 3: data['csig'] = (int1*rfit.GetParameter(8)+int2*rfit.GetParameter(5)+int3*rfit.GetParameter(2))/(int1+int2+int3) data['csigerr'] = calcMeanErr3(rfit.GetParameter(8), rfit.GetParError(8), rfit.GetParameter(6), rfit.GetParError(6), rfit.GetParameter(5), rfit.GetParError(5), rfit.GetParameter(3), rfit.GetParError(3), rfit.GetParameter(2), rfit.GetParError(2), rfit.GetParameter(0), rfit.GetParError(0)) # print "***********************************" # print "small=",rfit.GetParameter(8)," medium=",rfit.GetParameter(2), "broad",rfit.GetParameter(5)," mean=",meanRes # print "tsigma",tsigma," tmean",tmean," tconst",tconst # print "***********************************" elif ngaus ==2: data['csig'] = (int1*rfit.GetParameter(2)+int2*rfit.GetParameter(5))/(int1+int2) # 2 gaus version data['csigerr'] = calcMeanErr2(rfit.GetParameter(2), rfit.GetParError(2), rfit.GetParameter(0), rfit.GetParError(0), rfit.GetParameter(5), rfit.GetParError(5), rfit.GetParameter(3), rfit.GetParError(3)) # print "***********************************" # print "small=",rfit.GetParameter(5)," medium=",rfit.GetParameter(2)," mean=",meanRes # print "tsigma",tsigma," tmean",tmean," tconst",tconst # print "***********************************" else: data['csig'] = -9999 data['csigerr']=0. if ngaus == 3: data['nsig'] = rfit.GetParameter(8) data['nsigerr'] = rfit.GetParError(8) elif ngaus == 2: data['nsig'] = rfit.GetParameter(5) data['nsigerr'] = rfit.GetParError(5) data['bconst'] = rfit.GetParameter(0) data['nconst'] = rfit.GetParameter(3) data['bmean'] = rfit.GetParameter(1) data['nmean'] = rfit.GetParameter(4) if rfit.GetNDF()>0: data['chindf']=rfit.GetChisquare()/rfit.GetNDF() else: data['chindf']=0 data['bsig']=rfit.GetParameter(2) data['bsigerr']=rfit.GetParError(2) return data # }}} # {{{ fit sample ads distribution def fitadc(histo): expfit=ROOT.TF1("expfit","expo",9,36) expfit.SetLineColor(2) histo.Fit(expfit,"N+","",9,36) const=expfit.GetParameter(0) slope=expfit.GetParameter(1) slopeerr=expfit.GetParError(1) return slope, slopeerr # }}} # {{{ histo options def histopt(hist,hname): if hopt['scale'].get(hname,0)==1: if hist.GetEntries()>0: hist.Scale(1/hist.GetEntries()) if hopt['scale'].get(hname,0)==2: if hist.Integral()>0: hist.Scale(1/hist.Integral()) hist.GetXaxis().SetTitle(hopt['xtit'].get(hname,"")) hist.GetYaxis().SetTitle(hopt['ytit'].get(hname,"")) hist.GetYaxis().SetTitleOffset(1.3) if hopt['ymax'].get(hname,-9999)!=-9999: if hopt['ymin'].get(hname,-9999)!=-9999: hist.GetYaxis().SetRangeUser(hopt['ymin'][hname],hopt['ymax'][hname]*1.2) else: hist.GetYaxis().SetRangeUser(0.,hopt['ymax'][hname]*1.2) if hopt['xmax'].get(hname,-9999)!=-9999: if hopt['xmin'].get(hname,-9999)!=-9999: hist.GetXaxis().SetRangeUser(hopt['xmin'][hname],hopt['xmax'][hname]) else: hist.GetXaxis().SetRangeUser(0.,hopt['xmax'][hname]) return # }}} # {{{ set y max def setymax(hlist): print "set ymax" for hd in hists: # print hlist for hname in hlist: # print hname if type(hd[hname])==list: histos=hd[hname] else: histos=[hd[hname]] for hist in histos: factor=1 if hopt['scale'].get(hname,0)==1: if hist.GetEntries()>0: factor=hist.GetEntries() if hopt['scale'].get(hname,0)==2: if hist.Integral()>0: factor=hist.Integral() value=hist.GetBinContent(hist.GetMaximumBin()) hopt['ymax'][hname]=max([hopt['ymax'].get(hname,0),value/factor]) return # }}} # }}} # {{{ default values infilen="sampleAnaOut3282.root" fitn = "gaus" fitnn = 0 files = [] prelim = 0 ng = 2 quietfit = 1 derr=1 noarg=0 noB=0 ar=0 ne=0 double=0 pdir="" colors = [ROOT.kBlack,ROOT.kRed+1,ROOT.kGreen+3,ROOT.kBlue+1,ROOT.kMagenta+1,ROOT.kOrange+1,ROOT.kViolet-7] # }}} # {{{ argument parsing #argument parsing: print len(sys.argv)-1 for iarg in range(len(sys.argv)): arg = sys.argv[iarg] if arg == "-pdir" and noarg==0: pdir = sys.argv[iarg+1] if arg == "-h" and noarg==0: print "-pdir\tdirectory for saving the pictures" print "-files\tlist of files used for plotting. no further args allowed after this" print "-3g\tdo the fit with 3 gaussians instead of 2" print "-noB\tDiffusion curv for no B-field" print "-ar\tDiffusion curve for argon" print "-ne\tDiffusion curve for neon" exit(0) if arg =="-3g": ng=3 if arg=="-noB": noB=1 if arg=="-ar": ar=1 ne=0 if arg=="-ne": ar=0 ne=1 if arg=="-prel": prelim=1 if noarg == 1: files.append(arg) if arg == "-files": noarg = 1 print "picture directory:\t\t",pdir print "Number of gaussians for fit:\t",ng print "No magnetic field:\t\t",noB print "Gas, Ar/Ne:\t\t\t",ar,"/",ne print "Prelim:\t\t\t\t",prelim print files nix=raw_input("Happy with this settings?") # print sys.argv[iarg], iarg # }}} # {{{ get infile and root options ROOT.gROOT.ProcessLine('gStyle->SetPalette(1)') ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') #ROOT.gROOT.ProcessLine('gStyle->SetOptFit(1111)') infile = [None]*len(files) for ifile in range(len(files)): print ifile infile[ifile] = ROOT.TFile(files[ifile]) # }}} # {{{ build legend, set gains nlegend = [None]*len(files) gains = [None]*len(files) simar=[None]*len(files) for j in range(len(files)): nlegend[j]="something something" gains[j] = 9999 if files[j].find("Out2453")>=0: nlegend[j]="gain 1562" gains[j]=1560 if files[j].find("Out555")>=0: nlegend[j]="gain 1562 (sim)" gains[j]=1560 if files[j].find("Out2400")>=0: nlegend[j]="gain 2110" gains[j]=2000 if files[j].find("Out2464")>=0: nlegend[j]="gain 1155" gains[j]=1150 if files[j].find("Out333")>=0: nlegend[j]="gain 2110 (sim)" gains[j]=2110 if files[j].find("Out2451")>=0: nlegend[j]="gain 1562" gains[j]=1560 if files[j].find("Out666")>=0: nlegend[j]="gain 1155 (sim)" gains[j]=1150 if files[j].find("Out2621")>=0: nlegend[j]="gain 632" gains[j]=632 if files[j].find("Out2284")>=0: nlegend[j]="gain 3864" gains[j]=3864 if files[j].find("Out2280")>=0: nlegend[j]="gain 2857" gains[j]=2857 if files[j].find("cosmics_70") >=0: nlegend[j]="gain 632" gains[j]=632 if files[j].find("cosmics_72") >=0: nlegend[j]="gain 1150" gains[j]=1150 if files[j].find("cosmics_73") >=0: nlegend[j]="gain 1560" gains[j]=1560 if files[j].find("cosmics_74") >=0: nlegend[j]="gain 2113" gains[j]=2113 if files[j].find("cosmics_75") >=0: nlegend[j]="gain 2857" gains[j]=2857 if files[j].find("cosmics_76") >=0: nlegend[j]="gain 3864" gains[j]=3864 if files[j].find("cosmics_82") >=0: nlegend[j]="gain 1512" gains[j]=1512 if files[j].find("sim") >=0: nlegend[j]+=" sim" if files[j].find("newreco")>=0: nlegend[j]+=" new clustering" gains[j]+=.3 if files[j].find("smoothed")>=0: nlegend[j]+=" unbiased" gains[j]+=.2 if files[j].find("biased")>=0: nlegend[j]+=" biased" gains[j]+=.1 if files[j].find("oldres")>=0: nlegend[j]+=" old residual" gains[j]+=.4 if files[j].find("newcluster")>=0: nlegend[j]+=" new clustering" gains[j]+=.3 if files[j].find("3775")>=0: nlegend[j]="gain 1512" gains[j]=1512 if files[j].find("newreco")>=0: nlegend[j]+=" new clustering" gains[j]+=.3 if files[j].find("smoothed")>=0: nlegend[j]+=" unbiased" gains[j]+=.2 if files[j].find("biased")>=0: nlegend[j]+=" biased" gains[j]+=.1 if files[j].find("cosmics_84") >=0: nlegend[j]="gain 2781" gains[j]=2781 if files[j].find("cosmics_85") >=0: nlegend[j]="gain 3771" gains[j]=3771 if files[j].find("cosmics_86") >=0: nlegend[j]="gain 5113" gains[j]=5113 if files[j].find("65ar81") >0: nlegend[j]="gain 1100 (81), field 234 (65)" gains[j]=1100 if files[j].find("65ar85") >0: nlegend[j]="gain 3700 (85), field 234 (65)" gains[j]=3700 if files[j].find("65ar83") >0: nlegend[j]="gain 2000 (83), field 234 (65)" gains[j]=2000 if files[j].find("85ar81") >0: nlegend[j]="gain 1100 (81), field 306 (85)" gains[j]=1100 if files[j].find("86ar81") >0: nlegend[j]="gain 1100 (81), field 309 (86)" gains[j]=1100 if files[j].find("86ar82")>0: nlegend[j]="gain 1500 (82), field 309 (86)" gains[j]=1500 if files[j].find("86ar83")>0: nlegend[j]="gain 2000 (83), field 309 (86)" gains[j]=2000 if files[j].find("86ar84")>0: nlegend[j]="gain 2700 (84), field 309 (86)" gains[j]=2700 if files[j].find("86ar85")>0: nlegend[j]="gain 3700 (85), field 309 (86)" gains[j]=3700 if files[j].find("86ar86")>0: nlegend[j]="gain 5100 (86), field 309 (86)" gains[j]=5100 if files[j].find("100ar85") >0: nlegend[j]="gain 3700 (85), field 360 (100)" gains[j]=3700 if files[j].find("100ar83") >0: nlegend[j]="gain 2000 (83), field 360 (100)" gains[j]=2000 if files[j].find("100ne70")>0: nlegend[j]="gain 600 (70), field 360 (100)" gains[j]=600 if files[j].find("100ne71")>0: nlegend[j]="gain 800 (71), field 360 (100)" gains[j]=800 if files[j].find("100ne72")>0: nlegend[j]="gain 1100 (72), field 360 (100)" gains[j]=1100 if files[j].find("100ne73")>0: nlegend[j]="gain 1500 (73), field 360 (100)" gains[j]=1500 if files[j].find("100ne74")>0: nlegend[j]="gain 2100 (74), field 360 (100)" gains[j]=2100 if files[j].find("100ne75")>0: nlegend[j]="gain 2800 (75), field 360 (100)" gains[j]=2800 if files[j].find("100ne76")>0: nlegend[j]="gain 3800 (76), field 360 (100)" gains[j]=3800 if files[j].find("no_thetaCut")>0: nlegend[j]+=" no th Cut" if files[j].find("noRadCut")>0: nlegend[j]+=" no rad Cut" if files[j].find("nov")>0: nlegend[j]+=" Nov" if files[j].find("apr")>0: nlegend[j]+=(" Apr") if files[j].find("old")>0: nlegend[j]+=(" old") if files[j].find("cp20")>0: nlegend[j]+=(" cp20") if files[j].find("MC")>0: nlegend[j]+=(" MC") simar[j]=1 else: simar[j]=0 # }}} # {{{ define and get hists, define tgraphs if noB==1: if ar==1: DIFFX = 0.0227*10000 #(mu m/ mu s) without B-Field else: if ar==1: DIFFX = 0.020534*10000 #(mu m/ mu s)Argon with B-Field if ne==1: DIFFX = 0.019833*10000#(mu m/ mu s)Neon with B-Field diffT = ROOT.TF1("diffT","[0]*TMath::Sqrt(x)",0,75) diffT.SetNpx(1000) diffT.SetLineWidth(2) diffT.SetLineStyle(4) diffT.SetParameter(0,DIFFX) diffT.SetTitle("") hists=[] zSlices=[] dummyh=ROOT.TH1D("d","d",1,0.,1.) graphs=[] hopt={'scale':{},'same':{},'legend':{},'xtit':{},'ytit':{},'ymax':{},'ymin':{},'xmin':{},'xmax':{}} topt={'marker':{},'xtit':{},'ytit':{},'ymax':{},'ymin':{},'xmax':{},'xmin':{}} rSlices = (6.6,8.2,9.8,11.4,13,15) xySlices = (-10,-5,0,5,10,20) cpSlices = (7,14,21,28,35,42) thSlices = (90,95,100,105,110,120) print "files:",files for i in range(len(files)): hdict={} tdict={} hdict['clpertrack']=infile[i].Get("hclusterpertrack") hopt['scale']['clpertrack']=2 hopt['xtit']['clpertrack']="Cluster per Track" hopt['ytit']['clpertrack']="Normalized Counts (Integral)" hdict['clusteramp']=infile[i].Get("hClusterAmp5") hopt['scale']['clusteramp']=2 hopt['xtit']['clusteramp']="Cluster Amplitude" hopt['ytit']['clusteramp']="Normalized Counts (Integral)" hdict['someinfo']=infile[i].Get("someinfo") hdict['clustermax']=infile[i].Get("hClusterMax3") hopt['scale']['clustermax']=2 hopt['xtit']['clustermax']="Max Cluster Amplitude" hopt['ytit']['clustermax']="Normalized Counts (Integral)" hdict['csnr']=infile[i].Get("csnr") hopt['scale']['csnr']=2 hopt['xtit']['csnr']="SNR (cluster)" hopt['ytit']['csnr']="Normalized Counts (Integral)" hdict['snr']=infile[i].Get("hsnr") hopt['scale']['snr']=2 hopt['xtit']['snr']="SNR" hopt['ytit']['snr']="Normalized Counts (Integral)" hdict['sigma']=infile[i].Get("pedsigma") hdict['mean']=infile[i].Get("pedmean") hdict['clustersize all']=infile[i].Get("hpadspercluster2D") hdict['clustersize all 3D']=infile[i].Get("hpadspercluster") hopt['scale']['clustersize all']=2 hopt['scale']['clustersize all 3D']=2 hdict['hresX first']=infile[i].Get("hresXfirst") hdict['hresX last']=infile[i].Get("hresXlast") # {{{ get zSlices hdict['slices']=infile[i].Get("hslices") if type(hdict['slices'])!=ROOT.TH1D: slicef= ROOT.TFile("anafiles/dummyslice.root") hdict['slices']=slicef.Get("hslices") print hdict['slices'].GetEntries() temp=[] if hdict['slices'].GetEntries()>6: hdict['slices'].Scale(1./float(hdict['slices'].GetEntries()/6)) for j in range(1,7): temp.append(hdict['slices'].GetBinContent(j)) zSlices.append(temp) print "file",i,"slices:",zSlices[i] # zSlices[i].pop(5) # }}} # {{{ loop over zslices hdict['clustersize']=[] hopt['scale']['clustersize']=1 hopt['same']['clustersize']=1 hopt['legend']['clustersize']=1 hdict['clustersize2d']=[] hopt['scale']['clustersize2d']=1 hopt['same']['clustersize2d']=1 hopt['legend']['clustersize2d']=1 hdict['residualX'+str(gains[i])+"_"+str(i)]=[] hopt['same']['residualX'+str(gains[i])+"_"+str(i)]=0 hopt['xmin']['residualX'+str(gains[i])+"_"+str(i)]=-5000. hopt['xmax']['residualX'+str(gains[i])+"_"+str(i)]=5000. hdict['residualZ'+str(gains[i])+"_"+str(i)]=[] hopt['same']['residualZ'+str(gains[i])+"_"+str(i)]=0 hopt['xmin']['residualZ'+str(gains[i])+"_"+str(i)]=-5000. hopt['xmax']['residualZ'+str(gains[i])+"_"+str(i)]=5000. hdict['snr slice']=[] hopt['scale']['snr slice']=2 hopt['xtit']['snr slice']="SNR" hopt['ytit']['snr slice']="Normalized Counts (Integral)" hdict['samperpad']=[] hopt['scale']['samperpad']=1 hdict['digiperpad']=[] hopt['scale']['digiperpad']=1 hdict['csnr slice']=[] hopt['scale']['csnr slice']=2 hopt['xtit']['csnr slice']="SNR (cluster)" hopt['ytit']['csnr slice']="Normalized Counts (Integral)" hdict['clustermax slice']=[] hopt['scale']['clustermax slice']=2 hopt['xtit']['clustermax slice']="Max Amplitude in Cluster" hopt['ytit']['clustermax slice']="Normalized Counts (Integral)" hdict['clusteramp slice']=[] hopt['scale']['clusteramp slice']=2 hopt['xtit']['clusteramp slice']="Cluster Amplitude" hopt['ytit']['clusteramp slice']="Normalized Counts (Integral)" hdict['RresidualsX'+str(gains[i])+"_"+str(i)]=[] hopt['same']['RresidualsX'+str(gains[i])+"_"+str(i)]=0 hopt['xmin']['RresidualsX'+str(gains[i])+"_"+str(i)]=-5000. hopt['xmax']['RresidualsX'+str(gains[i])+"_"+str(i)]=5000. for j in range(len(zSlices[i])): hdict['samperpad'].append(infile[i].Get("nsampperpad"+str(zSlices[i][j]))) hdict['digiperpad'].append(infile[i].Get("digiperpad"+str(zSlices[i][j]))) hdict['clustersize'].append(infile[i].Get("npadspercluster3D"+str(zSlices[i][j]))) hdict['clustersize2d'].append(infile[i].Get("npadspercluster2D"+str(zSlices[i][j]))) hdict['residualX'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResX Z:"+str(zSlices[i][j]))) hdict['residualZ'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResZ Z:"+str(zSlices[i][j]))) hdict['snr slice'].append(infile[i].Get("nsnr"+str(zSlices[i][j]))) hdict['csnr slice'].append(infile[i].Get("csnr"+str(zSlices[i][j]))) hdict['clustermax slice'].append(infile[i].Get("clustermax"+str(zSlices[i][j]))) hdict['clusteramp slice'].append(infile[i].Get("clusteramp"+str(zSlices[i][j]))) # }}} # {{{ loop over xy and cp slices hdict['XResidualsX'+str(gains[i])+"_"+str(i)]=[] hopt['same']['XResidualsX'+str(gains[i])+"_"+str(i)]=0 hopt['xmin']['XResidualsX'+str(gains[i])+"_"+str(i)]=-5000. hopt['xmax']['XResidualsX'+str(gains[i])+"_"+str(i)]=5000. hdict['YResidualsX'+str(gains[i])+"_"+str(i)]=[] hopt['same']['YResidualsX'+str(gains[i])+"_"+str(i)]=0 hopt['xmin']['YResidualsX'+str(gains[i])+"_"+str(i)]=-5000. hopt['xmax']['YResidualsX'+str(gains[i])+"_"+str(i)]=5000. hdict['cpResidualsX'+str(gains[i])+"_"+str(i)]=[] hopt['same']['cpResidualsX'+str(gains[i])+"_"+str(i)]=0 hopt['xmin']['cpResidualsX'+str(gains[i])+"_"+str(i)]=-5000. hopt['xmax']['cpResidualsX'+str(gains[i])+"_"+str(i)]=5000. # hdict['thResidualsX'+str(gains[i])+"_"+str(i)]=[] # hopt['same']['thResidualsX'+str(gains[i])+"_"+str(i)]=0 # hopt['xmin']['thResidualsX'+str(gains[i])+"_"+str(i)]=-2000. # hopt['xmax']['thResidualsX'+str(gains[i])+"_"+str(i)]=2000. for j in range(len(xySlices)): hdict['XResidualsX'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResX X:"+str(xySlices[j]))) hdict['YResidualsX'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResX Y:"+str(xySlices[j]))) for j in range(len(cpSlices)): hdict['cpResidualsX'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResX cp:"+str(cpSlices[j]))) # }}} for j in range(len(rSlices)): hdict['RresidualsX'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResX R:"+str(rSlices[j]))) # for j in range(len(thSlices)): # hdict['thResidualsX'+str(gains[i])+"_"+str(i)].append(infile[i].Get("StatsResX th:"+str(thSlices[j]))) # for his in hdict: # if hdict[his]!=ROOT.TH1: # hdict[his]=dummyh # print hdict.keys() hists.append(hdict) tdict['residunarrow']=ROOT.TGraphErrors(len(zSlices[i])) tdict['residunarrow'].SetTitle("Narrow Gaussian") topt['marker']['residunarrow']=21 topt['xtit']['residunarrow']="Drift Length Z (cm)" topt['ytit']['residunarrow']="Resolution (#mum)" tdict['residucomb'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['residucomb'].SetTitle("Mean") topt['marker']['residucomb']=20 topt['xtit']['residucomb']="Drift Length Z (cm)" topt['ytit']['residucomb']="Resolution (#mum)" tdict['residucomb Z'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['residucomb Z'].SetTitle("Mean") topt['marker']['residucomb Z']=20 topt['xtit']['residucomb']="Drift Length Z (cm)" topt['ytit']['residucomb']="Resolution (#mum)" tdict['snrmpv'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['snrmpv'].SetTitle("Signal to noise mpv") topt['xtit']['snrmpv']="Drift Length Z (cm)" topt['ytit']['snrmpv']="Most probable SNR" tdict['csnrmpv'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['csnrmpv'].SetTitle("Signal to noise mpv (cluster)") topt['xtit']['csnrmpv']="Drift Length Z (cm)" topt['ytit']['csnrmpv']="Most probable SNR" tdict['clustersize'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['clustersize'].SetTitle("Clustersize") topt['xtit']['clustersize']="Drift Length Z (cm)" topt['ytit']['clustersize']="Mean Clustersize" tdict['clusteramp'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['clusteramp'].SetTitle("Clusteramplitude (mpv)") topt['xtit']['clusteramp']="Drift Length Z (cm)" topt['ytit']['clusteramp']="Cluster Amplitude (mpv)" tdict['clustermax'] = ROOT.TGraphErrors(len(zSlices[i])) tdict['clustermax'].SetTitle("Max Amp in Cluster") topt['xtit']['clustermax']="Drift Length Z (cm)" topt['ytit']['clustermax']="Max Amp in Cluster (mpv)" tdict['resolpergain'] = ROOT.TGraphErrors(1) tdict['resolpergain'].SetTitle("Resolution per Gain") topt['xtit']['resolpergain']="Gain" topt['ytit']['resolpergain']="Resolution" topt['xmax']['resolpergain']=8000 tdict['XResidualsX'] = ROOT.TGraphErrors(len(xySlices)) tdict['XResidualsX'].SetTitle("Resolution in X Slices") topt['xtit']['XResidualsX']="X-Slice (cm)" topt['ytit']['XResidualsX']="Resolution (#mum)" tdict['YResidualsX'] = ROOT.TGraphErrors(len(xySlices)) tdict['YResidualsX'].SetTitle("Resolution in Y Slices") topt['xtit']['YResidualsX']="Y-Slice (cm)" topt['ytit']['YResidualsX']="Resolution (#mum)" tdict['cpResidualsX'] = ROOT.TGraphErrors(len(cpSlices)) tdict['cpResidualsX'].SetTitle("Resolution in Tracklength Slices") topt['xtit']['cpResidualsX']="Tracklength (Cluster)" topt['ytit']['cpResidualsX']="Resolution (#mum)" tdict['RresidualsX'] = ROOT.TGraphErrors(len(cpSlices)) tdict['RresidualsX'].SetTitle("Resolution Radial Slices") topt['xtit']['RresidualsX']="Radius (cm)" topt['ytit']['RresidualsX']="Resolution (#mum)" # tdict['thResidualsX'] = ROOT.TGraphErrors(len(thSlices)) # tdict['thResidualsX'].SetTitle("Resolution Theta Slices") # topt['xtit']['thResidualsX']="Theta (#circ)" # topt['ytit']['thResidualsX']="Resolution (#mum)" tdict['snrmpvpergain']=ROOT.TGraphErrors(1) tdict['snrmpvpergain'].SetTitle("SNR mpv") topt['xtit']['snrmpvpergain']="Gain" topt['ytit']['snrmpvpergain']="SNR mpv" topt['xmax']['snrmpvpergain']=8000 graphs.append(tdict) # }}} # {{{ fill tgraphs landaufit = ROOT.TF1("landaufit","landau") resolutions={} for i in range(len(files)): resolutions['residualX'+str(gains[i])+"_"+str(i)]=[] resolutions['residualZ'+str(gains[i])+"_"+str(i)]=[] onegainresol=0 totentries=0 resolutions['hresX first']=fitres(hists[i]['hresX first'],ng,double,[1.,.8,1.5,20.]) resolutions['hresX last']=fitres(hists[i]['hresX last'],ng,double,[1.,.8,1.5,20.]) # {{{ z slices for j in range(len(zSlices[i])): if j==0: offset = abs(-62.5-zSlices[i][j])/2#-zSlices[i][0] elif j+1==len(zSlices[i]) and zSlices[i][j]>10.78: # offset = offset+(2*zSlices[i][j]-72.8-zSlices[i][0]-zSlices[i][j-1])/2 # offset=abs(72.8-zSlices[i][j])/2 offset=72.8-abs(10.78-zSlices[i][j-1])/2 else: offset = offset+abs(zSlices[i][j-1]-zSlices[i][j]) print "**************************************************************************" print "offset=",offset,j,zSlices[i][j],zSlices[i][j-1] print "**************************************************************************" # {{{ residuals # {{{ X print 'fit:',hists[i]['residualX'+str(gains[i])+"_"+str(i)][j],"gain:",gains[i],j resdata=fitres(hists[i]['residualX'+str(gains[i])+"_"+str(i)][j],ng,double) resolutions['residualX'+str(gains[i])+"_"+str(i)].append(resdata) if resdata['csig']>0: onegainresol+=resdata['csig']*(hists[i]['residualX'+str(gains[i])+"_"+str(i)][j].GetEntries()) totentries+=(hists[i]['residualX'+str(gains[i])+"_"+str(i)][j].GetEntries()) # topt['ymax']['residunarrow']=max([topt['ymax'].get('residunarrow',0),resdata['nsig']]) topt['ymin']['residunarrow']=0. topt['xmax']['residunarrow']=80. topt['xmin']['residunarrow']=0. topt['ymax']['residucomb']=max([topt['ymax'].get('residucomb',0),resdata['csig']]) topt['ymax']['residunarrow']=topt['ymax']['residucomb'] topt['ymin']['residucomb']=0. topt['xmax']['residucomb']=80. topt['xmin']['residucomb']=0. graphs[i]['residunarrow'].SetPoint(j,offset,resdata['nsig']) graphs[i]['residucomb'].SetPoint(j,offset,resdata['csig']) if derr == 1: graphs[i]['residunarrow'].SetPointError(j,0.,resdata['nsigerr']) graphs[i]['residucomb'].SetPointError(j,0.,resdata['csigerr']) # }}} # {{{ Z resdata=fitres(hists[i]['residualZ'+str(gains[i])+"_"+str(i)][j],ng,double,[1.,0.8,1.5,20.]) resolutions['residualZ'+str(gains[i])+"_"+str(i)].append(resdata) topt['ymax']['residucomb Z']=max([topt['ymax'].get('residucomb',0),resdata['csig']]) topt['ymin']['residucomb Z']=0. topt['xmax']['residucomb Z']=80. topt['xmin']['residucomb Z']=0. graphs[i]['residucomb Z'].SetPoint(j,offset,resdata['csig']) if derr == 1: graphs[i]['residucomb Z'].SetPointError(j,0.,resdata['csigerr']) # }}} # }}} # {{{ snr hists[i]['snr slice'][j].Fit(landaufit,"NQM") print '****************************************************' print i,j,landaufit.GetParameter(1) print '****************************************************' graphs[i]['snrmpv'].SetPoint(j,offset,landaufit.GetParameter(1)) topt['ymax']['snrmpv']=max([topt['ymax'].get('snrmpv',0), landaufit.GetParameter(1)]) hists[i]['csnr slice'][j].Fit(landaufit,"NQM") graphs[i]['csnrmpv'].SetPoint(j,offset,landaufit.GetParameter(1)) topt['ymax']['csnrmpv']=max([topt['ymax'].get('csnrmpv',0), landaufit.GetParameter(1)]) # }}} # {{{ clustersize hist=hists[i]['clustersize'][j] if hist.GetMean()>0: graphs[i]['clustersize'].SetPoint(j,offset,hist.GetMean()) topt['ymax']['clustersize']=max([hist.GetMean(), topt['ymax'].get('clustersize',0)]) topt['ymin']['clustersize']=0. topt['xmax']['clustersize']=80. topt['xmin']['clustersize']=0. hist=hists[i]['clusteramp slice'][j] hist.Fit(landaufit,"NQM") graphs[i]['clusteramp'].SetPoint(j,offset,landaufit.GetParameter(1)) topt['ymax']['clusteramp']=max([topt['ymax'].get('clusteramp',0.),landaufit.GetParameter(1)]) topt['ymin']['clusteramp']=min([topt['ymin'].get('clusteramp',0.),landaufit.GetParameter(1)]) # graphs[i]['clusteramp'].SetPointError(j,0,landaufit.GetParError(1)) hist=hists[i]['clustermax slice'][j] hist.Fit(landaufit,"NQM") graphs[i]['clustermax'].SetPoint(j,offset,landaufit.GetParameter(1)) topt['ymax']['clustermax']=max([topt['ymax'].get('clustermax',0.),landaufit.GetParameter(1)]) topt['ymin']['clustermax']=min([topt['ymin'].get('clustermax',0.),landaufit.GetParameter(1)]) # graphs[i]['clustermax'].SetPointError(j,0,landaufit.GetParError(1)) # }}} hists[i]['snr'].Fit(landaufit,"NQM") graphs[i]['snrmpvpergain'].SetPoint(0,gains[i],landaufit.GetParameter(1)) topt['ymax']['snrmpvpergain']=max([topt['ymax'].get('snrmpvpergain',0),landaufit.GetParameter(1)]) graphs[i]['resolpergain'].SetPoint(0,gains[i],onegainresol/totentries) topt['ymax']['resolpergain']=max([topt['ymax'].get('resolpergain',0),onegainresol/totentries]) # }}} # {{{ X Slices resolutions['XResidualsX'+str(gains[i])+"_"+str(i)]=[] resolutions['YResidualsX'+str(gains[i])+"_"+str(i)]=[] resolutions['RresidualsX'+str(gains[i])+"_"+str(i)]=[] for j in range(len(xySlices)): offset=xySlices[j]-2.5 if j==len(xySlices): offset=xySlices[j]-7.5 resdata=fitres(hists[i]['XResidualsX'+str(gains[i])+"_"+str(i)][j],ng,double,[1.,0.8,1.5,20.]) resolutions['XResidualsX'+str(gains[i])+"_"+str(i)].append(resdata) graphs[i]['XResidualsX'].SetPoint(j,offset,resdata['csig']) graphs[i]['XResidualsX'].SetPointError(j,2.5,resdata['csigerr']) resdata=fitres(hists[i]['YResidualsX'+str(gains[i])+"_"+str(i)][j],ng,double) resolutions['YResidualsX'+str(gains[i])+"_"+str(i)].append(resdata) graphs[i]['YResidualsX'].SetPoint(j,offset,resdata['csig']) graphs[i]['YResidualsX'].SetPointError(j,2.5,resdata['csigerr']) # }}} # {{{ cp Slices resolutions['cpResidualsX'+str(gains[i])+"_"+str(i)]=[] for j in range(len(cpSlices)): offset=cpSlices[j]-3.5 resdata=fitres(hists[i]['cpResidualsX'+str(gains[i])+"_"+str(i)][j],ng,double) resolutions['cpResidualsX'+str(gains[i])+"_"+str(i)].append(resdata) graphs[i]['cpResidualsX'].SetPoint(j,offset,resdata['csig']) graphs[i]['cpResidualsX'].SetPointError(j,2.5,resdata['csigerr']) topt['ymax']['cpResidualsX']=max([topt['ymax'].get('resolpergain',0),onegainresol/totentries]) # }}} # {{{ th Slices # resolutions['thResidualsX'+str(gains[i])+"_"+str(i)]=[] # for j in range(len(thSlices)): # offset=thSlices[j] # resdata=fitres(hists[i]['thResidualsX'+str(gains[i])+"_"+str(i)][j],ng,double) # resolutions['thResidualsX'+str(gains[i])+"_"+str(i)].append(resdata) # graphs[i]['thResidualsX'].SetPoint(j,offset,resdata['csig']) # graphs[i]['thResidualsX'].SetPointError(j,2.5,resdata['csigerr']) # }}} for j in range (len(rSlices)): resdata=fitres(hists[i]['RresidualsX'+str(gains[i])+"_"+str(i)][j],ng,double,[1.,0.8,1.5,20.]) resolutions['RresidualsX'+str(gains[i])+"_"+str(i)].append(resdata) graphs[i]['RresidualsX'].SetPoint(j,rSlices[j],resdata['csig']) # graphs[i]['RresidualsX'].SetPointError(j,2.5,resdata['csigerr']) # }}} topt['ymax']['clustersize']*=1.2 topt['ymax']['residucomb']*=1.2 topt['ymax']['residunarrow']*=1.2 topt['xmax']['snrmpv']=80 for i in range(len(files)): if hists[i]['someinfo'].GetEntries()>0: hists[i]['sigma'].Scale(1/hists[i]['someinfo'].GetEntries()) hists[i]['mean'].Scale(1/hists[i]['someinfo'].GetEntries()) if len(files)==2: nx=hists[0]['sigma'].GetNbinsX() hists[0]['sigmadiff']=ROOT.TH1D("sigmadiff","difference of noise sigmas",nx,0.,nx) hists[1]['sigmadiff']=ROOT.TH1D("sigmadiff2","difference of noise sigmas fpn",nx,0.,nx) hopt['ymax']['sigmadiff']=300 hopt['ymin']['sigmadiff']=-300 hopt['xtit']['sigmadiff']="Arbitrary pad number" hopt['ytit']['sigmadiff']="Relative difference in %" hopt['legend']['sigmadiff']=0 hists[0]['allsigma']=ROOT.TH1D("all sigma 0","all sigma",500,0.,20) hists[1]['allsigma']=ROOT.TH1D("all sigma 1","all sigma",500,0.,20) hopt['xtit']['allsigma']="Sigma of Noise" hists[0]['meandiff']=ROOT.TH1D("meandiff","difference if noise means",nx,0.,nx) hists[1]['meandiff']=ROOT.TH1D("meandiff2","difference if noise means (fpn)",nx,0.,nx) hopt['ymax']['meandiff']=300 hopt['ymin']['meandiff']=-300 hopt['xtit']['meandiff']="Arbitrary pad number" hopt['ytit']['meandiff']="Relative difference in %" hopt['legend']['meandiff']=0 fpn=13 fpncounter=0 for i in range(hists[0]['sigma'].GetNbinsX()): fpncounter+=1 if hists[0]['sigma'].GetBinContent(i+1)>0: sdiff=(hists[0]['sigma'].GetBinContent(i+1)-hists[1]['sigma'].GetBinContent(i+1))/hists[0]['sigma'].GetBinContent(i+1) mdiff=(hists[0]['mean'].GetBinContent(i+1)-hists[1]['mean'].GetBinContent(i+1))/hists[0]['mean'].GetBinContent(i+1) else: sdiff=0 mdiff=0 if fpncounter != fpn: hists[0]['sigmadiff'].Fill(i,sdiff*100) # hists[1]['sigmadiff'].Fill(i,100) hists[0]['meandiff'].Fill(i,mdiff*100) # hists[1]['meandiff'].Fill(i,200) hists[0]['allsigma'].Fill(round(hists[0]['sigma'].GetBinContent(i+1),2)) hists[1]['allsigma'].Fill(round(hists[1]['sigma'].GetBinContent(i+1),2)) else: # hists[1]['sigmadiff'].Fill(i,sdiff*100+100) # hists[1]['meandiff'].Fill(i,mdiff*100+200) fpncounter=0 if fpn==13 and i > 20: fpn=25 else: fpn=13 # {{{ plot hists setymax(['clustermax','snr','csnr','snr slice','csnr slice','clustermax slice','clusteramp slice','clusteramp', 'clustersize','clustersize2d','clpertrack']) canvases={} legends={} canvascounter=0 fcounter=0 hcounter=0 rfit1 = ROOT.TF1("gaus1","gaus",-10000.,10000.) rfit2 = ROOT.TF1("gaus2","gaus",-10000.,10000.) for hd in hists: hcounter=0 for h in hd: if fcounter==0 or hopt['same'].get(h,1)==0: canvases[h]=ROOT.TCanvas("c"+str(canvascounter),h) legends[h]=ROOT.TLegend(0.45,.9-len(nlegend)*0.04,0.9,.9,"","NDC") legends[h].SetFillColor(0) canvascounter+=1 # {{{ if type(hd[h])==list : if type(hd[h])==list : if fcounter==0 or hopt['same'].get(h,1)==0: canvases[h].Divide(3,2) for i in range(len(hd[h])): canvases[h].cd(i+1) hd[h][i].SetLineColor(colors[fcounter]) histopt(hd[h][i],h) if i==0 and hopt['legend'].get(h,0)==1: legends[h].AddEntry(hd[h][i],nlegend[fcounter]) if fcounter==0 or hopt['same'].get(h,1)==0: hd[h][i].Draw() if resolutions.get(h,None)!=None: drawresdata(resolutions[h][i]) rfit1.SetParameter(0,resolutions[h][i]['bconst']) rfit1.SetParameter(1,resolutions[h][i]['bmean']) rfit1.SetParameter(2,resolutions[h][i]['bsig']) rfit1.SetLineColor(ROOT.kRed) rfit1.DrawCopy("same") rfit2.SetParameter(0,resolutions[h][i]['nconst']) rfit2.SetParameter(1,resolutions[h][i]['nmean']) rfit2.SetParameter(2,resolutions[h][i]['nsig']) rfit2.SetLineColor(ROOT.kBlue) rfit2.DrawCopy("same") if hopt['legend'].get(h,0)==1: legends[h].Draw() else: hd[h][i].Draw("same") if prelim==1: drawPrelim() # }}} # {{{ else else: canvases[h].cd() legends[h].AddEntry(hd[h],nlegend[fcounter]) hd[h].SetLineColor(colors[fcounter]) histopt(hd[h],h) if fcounter==0: hd[h].Draw() if resolutions.get(h,None)!=None: drawresdata(resolutions[h]) rfit1.SetParameter(0,resolutions[h]['bconst']) rfit1.SetParameter(1,resolutions[h]['bmean']) rfit1.SetParameter(2,resolutions[h]['bsig']) rfit1.SetLineColor(ROOT.kRed) rfit1.DrawCopy("same") rfit2.SetParameter(0,resolutions[h]['nconst']) rfit2.SetParameter(1,resolutions[h]['nmean']) rfit2.SetParameter(2,resolutions[h]['nsig']) rfit2.SetLineColor(ROOT.kBlue) rfit2.DrawCopy("same") if hopt['legend'].get(h,1)==1: legends[h].Draw() else: hd[h].Draw("same") if prelim==1: drawPrelim() # }}} hcounter+=1 fcounter+=1 if len(files)==2: print '*****************sigmas****************************' print hists[0]['allsigma'].GetMean(),hists[1]['allsigma'].GetMean() print '*****************sigmas****************************' # }}} # {{{ plot graphs fcounter=0 #canvascounter=0 for gd in graphs: gcounter=hcounter for g in gd: if fcounter==0: canvases[g]=ROOT.TCanvas("c"+str(canvascounter),g) legends[g]=ROOT.TLegend(0.45,.9-len(nlegend)*0.04,0.9,.9,"","NDC") legends[g].SetFillColor(0) canvascounter+=1 canvases[g].cd() legends[g].AddEntry(gd[g],nlegend[fcounter],"P") gd[g].SetMarkerStyle(topt['marker'].get(g,20)) gd[g].GetXaxis().SetTitle(topt['xtit'].get(g,"")) gd[g].GetXaxis().SetTitleOffset(1.3) gd[g].GetYaxis().SetTitle(topt['ytit'].get(g,"")) gd[g].GetYaxis().SetTitleOffset(1.3) gd[g].SetMarkerColor(colors[fcounter]) if topt['ymax'].get(g,-99999)!=-99999: gd[g].GetYaxis().SetRangeUser(topt['ymin'].get(g,0.),topt['ymax'][g]*1.5) if topt['xmax'].get(g,-99999)!=-99999: gd[g].GetHistogram().SetBins(500,0.,topt['xmax'][g]) if fcounter==0: gd[g].Draw("AP") legends[g].Draw() if prelim==1: drawPrelim() else: gd[g].Draw("P") gcounter+=1 fcounter+=1 # }}} fcounter=0 for i in range(len(files)): canvases['resol']=ROOT.TCanvas("c"+str(canvascounter),"NarrowAndMean"+str(gains[i])) legends['resol'+str(i)]=ROOT.TLegend(0.45,.9-len(nlegend)*0.08,0.9,.9,"","NDC") legends['resol'+str(i)].SetFillColor(0) canvascounter+=1 canvases['resol'].cd() graphs[i]['residunarrow'].Draw("AP") graphs[i]['residucomb'].Draw("P") diffT.Draw("same") legends['resol'+str(i)].AddEntry(graphs[i]['residunarrow'],"Narrow","P") legends['resol'+str(i)].AddEntry(graphs[i]['residucomb'],"Combined","P") legends['resol'+str(i)].Draw() canvases['residucomb'].cd() diffT.Draw("same") canvases['residunarrow'].cd() diffT.Draw("same") canvases['residucomb Z'].cd() diffT.Draw("same") # {{{ print for c in canvases: # print pdir+"/"+c.GetTitle()+".pdf" if pdir!="": canvases[c].SaveAs(pdir+"/"+canvases[c].GetTitle()+".pdf") canvases[c].SaveAs(pdir+"/"+canvases[c].GetTitle()+".C") fitfile=open(pdir+"/"+"fitvalues.txt",'w') for hd in hists: for h in hd: if resolutions.get(h,None)!=None: fitfile.write(h+"\n") if type (hd[h])==list: for i in range(len(hd[h])): fitfile.write(str(round(resolutions[h][i]['bsig'],2))+"&") fitfile.write(str(round(resolutions[h][i]['nsig'],2))+"&") fitfile.write(str(round(resolutions[h][i]['csig'],2))+"&") fitfile.write(str(round(resolutions[h][i]['chindf'],2))+"&") fitfile.write(str(round((resolutions[h][i]['nsig']*resolutions[h][i]['nconst'])/(resolutions[h][i]['bsig']*resolutions[h][i]['bconst']),2))+"\\\\\\hline\n") else: fitfile.write(str(round(resolutions[h]['bsig'],2))+"&") fitfile.write(str(round(resolutions[h]['nsig'],2))+"&") fitfile.write(str(round(resolutions[h]['csig'],2))+"&") fitfile.write(str(round(resolutions[h]['chindf'],2))+"&") fitfile.write(str(round((resolutions[h]['nsig']*resolutions[h]['nconst'])/(resolutions[h]['bsig']*resolutions[h]['bconst']),2))+"\\\\\\hline\n") fitfile.write("\n") # }}} # {{{ now some extra stuff if len(files)==1: tcanv=ROOT.TCanvas("temp","temp") tcanv.cd() h='residualX'+str(gains[0])+"_"+str(0) i=0 hdict['residualX'+str(gains[0])+"_"+str(0)][0].Draw() drawresdata(resolutions['residualX'+str(gains[0])+"_"+str(0)][0]) rfit1.SetParameter(0,resolutions[h][i]['bconst']) rfit1.SetParameter(1,resolutions[h][i]['bmean']) rfit1.SetParameter(2,resolutions[h][i]['bsig']) rfit1.SetLineColor(ROOT.kRed) rfit1.DrawCopy("same") rfit2.SetParameter(0,resolutions[h][i]['nconst']) rfit2.SetParameter(1,resolutions[h][i]['nmean']) rfit2.SetParameter(2,resolutions[h][i]['nsig']) rfit2.SetLineColor(ROOT.kBlue) rfit2.DrawCopy("same") # }}} if prelim==1: drawPrelim() nix=raw_input('Press enter to finish')