import ROOT, glob, math, sys, os from ROOT import std 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[0],1))+" #mum") y-=ystep text.DrawLatex(0.15,y,"#sigma_{broad}="+str(round(data[4],1))+" #mum") y-=ystep text.DrawLatex(0.15,y,"#sigma_{mean}= "+str(round(data[1],1))+" #mum") y-=ystep text.DrawLatex(0.15,y,"#chi^{2}/ndf="+str(round(data[6],1))) 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) : # histo.SetStats(ROOT.kFALSE) data={} histo.SetLineColor(colors[i]) histo.SetFillColor(ROOT.kAzure-8) testfit = ROOT.TF1("testfitY"+str(1),"gaus",-10000.,10000.) testfit.SetLineColor(2) histo.Fit(testfit, "N+", "", -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] = testfit.GetParameter(1) tmean[0] = tmean[1] = tmean[2] = 0 tsigma[0] = tsigma[1] = tsigma[2] = testfit.GetParameter(2) #medium gaus rfit.SetParameter(0,tconst[0]) rfit.SetParameter(1,0.) # rfit.SetParLimits(1,tmean[0]-0.03, tmean[0]+0.03) rfit.SetParLimits(1,0.,0.) rfit.FixParameter(1,0.) rfit.SetParameter(2,tsigma[0]) if ngaus == 3 : rfit.SetParLimits(2,tsigma[0]*0.5, tsigma[0]*2) if ngaus == 2 : # rfit.SetParLimits(2,0.01, tsigma[0]*1.2) rfit.SetParLimits(2,tsigma[0]*2., tsigma[0]*12.) #broadest gaus # if ng=2 small gaus rfit.SetParameter(3, tconst[1]/1.5) rfit.SetParameter(4, 0.) # rfit.SetParLimits(4, tmean[1]-0.03, tmean[1]+0.03) rfit.FixParameter(4,0.) # rfit.SetParLimits(4, 0., 0.) rfit.SetParameter(5, tsigma[1]) if ngaus == 3 : rfit.SetParLimits(5, tsigma[1]*0.8, tsigma[1]*2) if ngaus == 2 : rfit.SetParLimits(5, 0.0, tsigma[1]*0.8) if ngaus == 3: #smallest gaus rfit.SetParameter(6, tconst[2]*4) rfit.SetParameter(7,0.) rfit.FixParameter(7,0.) # rfit.SetParLimits(7, tmean[2]-0.03, tmean[2]+0.03) # rfit.SetParLimits(7, 0., 0.) rfit.SetParameter(8, tsigma[2]*0.2) rfit.SetParLimits(8, 50,tsigma[2]*0.6) rfit.SetLineWidth(1) histo.Fit(rfit, "+", "", -10000,10000) # {{{ 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) # }}} 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.Draw("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.Draw("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: # meanRes = (int1*rfit.GetParameter(8)+int2*rfit.GetParameter(5))/(int1+int2) # 3 gaus version meanRes = (int1*rfit.GetParameter(8)+int2*rfit.GetParameter(5)+int3*rfit.GetParameter(2))/(int1+int2+int3) meanerr = 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: meanRes = (int1*rfit.GetParameter(2)+int2*rfit.GetParameter(5))/(int1+int2) # 2 gaus version meanerr = 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: meanRes = -9999 if ngaus == 3: rsig = rfit.GetParameter(8) sigerr = rfit.GetParError(8) elif ngaus == 2: rsig = rfit.GetParameter(5) sigerr = rfit.GetParError(5) if rfit.GetNDF()>0: data['chindf']=rfit.GetChisquare()/rfit.GetNDF() else: data['chindf']=0 return rsig,meanRes,sigerr,meanerr, rfit.GetParameter(2), rfit.GetParError(2), 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 # }}} infilen="sampleAnaOut3282.root" fitn = "gaus" fitnn = 0 files = [] prelim = 1 ng = 2 quietfit = 1 derr=1 noarg=0 noB=0 ar=0 ne=0 double=0 pdir="pics/cosmic/comb" # {{{ argument parsing #argument parsng: 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 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 files nix=raw_input("Happy with this settings?") # print sys.argv[iarg], iarg # }}} 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]) #zSlices = (-50, -40, -30, -20, -10, 0, 20 ) #zSlices = ( -52, -41, -29, -16, -2, 20 ) # cosmic data #zSlices = ( 20, 30, 40, 49, 57, 70 ) #beam data zSlices = [] # {{{ define histo arrays and canvas hdigiperpad = [None]*len(files) #hnsamperpad = [None]*len(files) hsnr = [None]*len(files) #hsamperdigi = [None]*len(files) hClusterMax = [None]*len(files) hresidu = [None]*len(files) hslices = [None]*len(files) hsampletot = [None]*len(files) nsnr = [None]*len(files) nsamperdigi = [None]*len(files) ndigiperpad = [None]*len(files) nsamperpad = [None]*len(files) nsamperpadpsa = [None]*len(files) nresidu = [None]*len(files) nsamperclust = [None]*len(files) legend = [None]*16 canvas = [None]*16 rescanvas = [None]*len(files) # }}} # {{{ define tgraphs snrPerGain = ROOT.TGraphErrors(len(files)) snrPerGain.SetName("snrpergain") snrPerGain.SetTitle("") snrPerGainmpv = ROOT.TGraphErrors(len(files)) snrPerGainmpv.SetName("snrpergainmpv") snrPerGainmpv.SetTitle("") residupergain = ROOT.TGraphErrors(len(files)) residupergain.SetName("residupergain") residupergain.SetTitle("") residupergaincomb = ROOT.TGraphErrors(len(files)) residupergaincomb.SetName("residupergaincomb") residupergaincomb.SetTitle("") adcvsgain = ROOT.TGraphErrors(len(files)) adcvsgain.SetName("adcvsgain") adcvsgain.SetTitle("") snrPerGainSlice = [None]*len(files) snrPerGainSlicempv = [None]*len(files) residuSlice = [None]*len(files) residuSlicecomb = [None]*len(files) 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("") # }}} cnames = [] cnames.append("digis_per_pad") # 0 cnames.append("clustersize") # 1 cnames.append("snr") # 2 cnames.append("clustersize2d") # 3 cnames.append("max_sampl_amp") # 4 cnames.append("snr_sliced") # 5 cnames.append("samples_per_digi") # 6 cnames.append("digis_per_pad") # 7 cnames.append("samples_per_pad2") # 8 cnames.append("samples_per_pad_psa") # 9 cnames.append("empty") # 10 cnames.append("empty") # 11 cnames.append("residuals") # 12 cnames.append("residuals_narrow_vs_z") # 13 cnames.append("residuals_combined_vs_z") # 14 cnames.append("samples_per_cluster") # 15 hists=[] dummyh=ROOT.TH1D("d","d",1,0.,1.) for i in range(len(files)): hdict={} # hdigiperpad[i] = infile[i].Get("hdigiperpad") # hnsamperpad[i] = infile[i].Get("hnsamperpad") hsnr[i] = infile[i].Get("hsnr") # hsamperdigi[i] = infile[i].Get("hsamperdigi") hClusterMax[i] = infile[i].Get("hClusterMax") if type(hClusterMax[i])!=ROOT.TH1D: hClusterMax[i] = infile[i].Get("hClusterMax3") hresidu[i] = infile[i].Get("hresX") if type(hresidu[i])!=ROOT.TH1D: hresidu[i]=dummyh hslices[i] = infile[i].Get("hslices") hsampletot[i] = infile[i].Get("hsampletot") hslices[i].Draw() temp = [] for k in range(1,7): print hslices[i].GetBinContent(k) temp.append(hslices[i].GetBinContent(k)) print temp zSlices.append(temp) print "slices:",zSlices, zSlices[i] snrPerGainSlice[i] = ROOT.TGraphErrors(len(zSlices[i])) snrPerGainSlicempv[i]= ROOT.TGraphErrors(len(zSlices[i])) residuSlice[i] = ROOT.TGraphErrors(len(zSlices[i])) residuSlicecomb[i] = ROOT.TGraphErrors(len(zSlices[i])) snrPerGainSlice[i].SetTitle("") snrPerGainSlicempv[i].SetTitle("") residuSlice[i].SetTitle("Narrow Gaussian") residuSlicecomb[i].SetTitle("Mean") if len(files) == 1 : residuSlice[i].SetTitle("") residuSlicecomb[i].SetTitle("") nsnr[i] = [None]*len(zSlices[i]) nsamperpad[i] = [None]*len(zSlices[i]) ndigiperpad[i] = [None]*len(zSlices[i]) nsamperdigi[i] = [None]*len(zSlices[i]) nsamperpadpsa[i] = [None]*len(zSlices[i]) nresidu[i] = [None]*len(zSlices[i]) nsamperclust[i] = [None]*len(zSlices[i]) hdict['clustersize']=[] hdict['clustersize2d']=[] for j in range(len(zSlices[i])): print zSlices[i][j] nsnr[i][j] =(infile[i].Get("nsnr"+str(zSlices[i][j]))) nsamperpad[i][j] =(infile[i].Get("nsampperpad"+str(zSlices[i][j]))) ndigiperpad[i][j] =(infile[i].Get("digiperpad"+str(zSlices[i][j]))) nsamperdigi[i][j] =(infile[i].Get("samperdigi"+str(zSlices[i][j]))) nsamperpadpsa[i][j] =(infile[i].Get("nsampperpad psa"+str(zSlices[i][j]))) nresidu[i][j] =(infile[i].Get("StatsResX Z:"+str(zSlices[i][j]))) nsamperclust[i][j] =(infile[i].Get("samperclust"+str(zSlices[i][j]))) if type(nsamperclust[i][j])!=ROOT.TH1D: nsamperclust[i][j]=dummyh hdict['clustersize'].append(infile[i].Get("npadspercluster3D"+str(zSlices[i][j]))) hdict['clustersize2d'].append(infile[i].Get("npadspercluster2D"+str(zSlices[i][j]))) if type(nresidu[i][j]) != ROOT.TH1D: print "didnt find residual histo, try another name, file",i nresidu[i][j] =(infile[i].Get("StatsResX"+str(zSlices[i][j]))) nresidu[i][j].SetTitle("") nresidu[i][j].GetXaxis().SetTitle("Residual X (#mum)") nresidu[i][j].GetXaxis().SetTitleSize(0.05) hists.append(hdict) #zSlices = (-50, -40, -30, -20, -10, 0, 10 ) #zSlices = ( -52, -41, -29, -16, -2, 10 ) #cosmic data #zSlices = ( 10, 30, 40, 49, 57, 67 ) # beam data c99 = ROOT.TCanvas() c99.cd() nlegend = [None]*len(files) gains = [None]*len(files) # {{{ build legend, set gains for j in range(len(files)): nlegend[j]=" " 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("86ar81") >0: nlegend[j]="gain 1100 (81), field 309 (86)" gains[j]=1100 if files[j].find("65ar81") >0: nlegend[j]="gain 1100 (81), field 234 (65)" gains[j]=1100 if files[j].find("85ar81") >0: nlegend[j]="gain 1100 (81), field 306 (85)" gains[j]=1100 if files[j].find("100ar85") >0: nlegend[j]="gain 3770 (85), field 360 (100)" gains[j]=3770 if files[j].find("86ar85")>0: nlegend[j]="gain 3770 (85), field 309 (86)" gains[j]=3770 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") cn = "rc"+str(j) ct = "res"+str(gains[j])+"_"+str(j) rescanvas[j]=ROOT.TCanvas(cn,ct) rescanvas[j].Divide(3,2) for i in range(len(canvas)): canvas[i] = ROOT.TCanvas() canvas[i].SetTitle(cnames[i]) if i == 13: legend[i] = ROOT.TLegend(0.7,.85-len(nlegend)*0.05,0.9,.9,"","NDC") else: legend[i] = ROOT.TLegend(0.65,.9-len(nlegend)*0.05,0.9,.9,"","NDC") legend[i].SetFillColor(0) # }}} fit = ROOT.TF1("fitfunc","landau") csnrmpvslice = ROOT.TCanvas("csnrmpvslice","snr mpv sliced") csnrslice = ROOT.TCanvas("csnrslice","snr mean sliced") c99.cd() # {{{ work on hists colors = [ROOT.kBlack,ROOT.kRed+1,ROOT.kGreen+3,ROOT.kBlue+1,ROOT.kMagenta+1,ROOT.kOrange+1,ROOT.kViolet-7] # {{{ set up residual data array resdata = [] for k in range(len(files)): resdata.append([]) print resdata for l in range(len(zSlices[k])): print resdata resdata[k].append([]) for m in range(7): resdata[k][l].append(k+l+m) print resdata # }}} maxmpvsnr=0 maxmeansnr=0 for i in range(len(files)): # {{{ hdigiperpad # hdigiperpad[i].SetStats(ROOT.kFALSE) # if hdigiperpad[i].GetEntries() > 0 : # hdigiperpad[i].Scale(1/hdigiperpad[i].GetEntries()) # hdigiperpad[i].SetLineColor(colors[i]) # }}} # {{{ hnsamperpad # # hnsamperpad[i].SetStats(ROOT.kFALSE) # if hnsamperpad[i].GetEntries() > 0 : # hnsamperpad[i].Scale(1/hnsamperpad[i].GetEntries()) # hnsamperpad[i].SetLineColor(colors[i]) # }}} # {{{ hsnr fit and snrPerGain and snrPerGainmpv hsnr[i].Fit(fit,"NOQ") hsnr[i].SetStats(ROOT.kFALSE) if hsnr[i].Integral() > 0: hsnr[i].Scale(1/hsnr[i].Integral()) hsnr[i].SetLineColor(colors[i]) snrPerGain.SetPoint(i,gains[i],hsnr[i].GetMean()) print "**********************************************" print "SNR:",fit.GetParameter(1),"+-",fit.GetParError(1)," gain:", gains[i] print "**********************************************" # snrPerGain.SetPointError(i,0,hsnr[i].GetMeanError()) snrPerGain.SetPointError(i,0,0) # hsnr[i].Fit("landau") snrPerGainmpv.SetPoint(i,gains[i],fit.GetParameter(1)) snrPerGainmpv.SetPointError(i,0,fit.GetParError(1)) print "singeldings:",fit.GetParameter(1),fit.GetParError(1) # snrPerGainmpv.SetPointError(i,0,0) # }}} # {{{ hsamperdigi # hsamperdigi[i].SetStats(ROOT.kFALSE) # hsamperdigi[i].SetLineColor(colors[i]) # if hsamperdigi[i].GetEntries() > 0 : # hsamperdigi[i].Scale(1/hsamperdigi[i].GetEntries()) # }}} # {{{ hclustermax hClusterMax[i].SetStats(ROOT.kFALSE) if hClusterMax[i].Integral() > 0: hClusterMax[i].Scale(1/hClusterMax[i].Integral()) hClusterMax[i].SetLineColor(colors[i]) # }}} # {{{ residuals # sig, meansig = fitres(hresidu[i],ng,0) # residupergain.SetPoint(i,gains[i],sig) # residupergaincomb.SetPoint(i,gains[i],meansig) # }}} # {{{ sample distribution # }}} # {{{ zslices csig = 0 csige = 0 cmsig = 0 cmsige = 0 nslices = 0 totentries = 0 bsig=0 bsige=0 print zSlices[i] for j in range(len(zSlices[i])): totentries+=nresidu[i][j].GetEntries() for j in range(len(zSlices[i])): # if nsnr[i][j].__class__ == ROOT.TF1: nsnr[i][j].Fit(fit,"NOQ") # print i,j,fit.GetParameter(1), nsnr[i][j].GetMean() #fitting of residual slice histos banane sig, meansig, sige, msige ,bsig, bsige, data= fitres(nresidu[i][j],ng,double) resdata[i][j][0]=sig resdata[i][j][1]=meansig resdata[i][j][2]=sige resdata[i][j][3]=msige resdata[i][j][4]=bsig resdata[i][j][5]=bsige resdata[i][j][6]=data['chindf'] if j==0: offset = abs(-62.5-zSlices[i][j])/2#-zSlices[i][0] elif j+1==len(zSlices[i]): offset = offset+(2*zSlices[i][j]-72.8-zSlices[i][0]-zSlices[i][j-1])/2 else: offset = offset+abs(zSlices[i][j-1]-zSlices[i][j]) snrPerGainSlicempv[i].SetPoint(j,offset,fit.GetParameter(1)) snrPerGainSlice[i].SetPoint(j,offset,nsnr[i][j].GetMean()) maxmpvsnr=max([maxmpvsnr,fit.GetParameter(1)]) maxmeansnr=max([maxmeansnr,nsnr[i][j].GetMean()]) print"****************************************************************" print "offset=",offset,j,"pos=",zSlices[i][j]+offset,len(zSlices[i]) print"****************************************************************" residuSlice[i].SetPoint(j,offset,sig) residuSlicecomb[i].SetPoint(j,offset,meansig) if derr == 1: residuSlicecomb[i].SetPointError(j,0.,msige) residuSlice[i].SetPointError(j,0.,sige) if sig > 0 and meansig > 0: csig += sig*nresidu[i][j].GetEntries() csige += pow(nresidu[i][j].GetEntries()*sige/totentries,2) cmsig += meansig*nresidu[i][j].GetEntries() cmsige += pow(nresidu[i][j].GetEntries()*msige/totentries,2) nslices+=1 nsnr[i][j].SetStats(ROOT.kFALSE) if nsnr[i][j].Integral() > 0 : nsnr[i][j].Scale(1/nsnr[i][j].Integral()) nsamperpad[i][j].SetStats(ROOT.kFALSE) if nsamperpad[i][j].GetEntries() > 0 : nsamperpad[i][j].Scale(1/nsamperpad[i][j].GetEntries()) nsamperpad[i][j].SetLineColor(colors[i]) ndigiperpad[i][j].SetStats(ROOT.kFALSE) if ndigiperpad[i][j].GetEntries() > 0 : ndigiperpad[i][j].Scale(1/ndigiperpad[i][j].GetEntries()) ndigiperpad[i][j].SetLineColor(colors[i]) nsamperclust[i][j].SetStats(ROOT.kFALSE) if nsamperclust[i][j].GetEntries() > 0 : nsamperclust[i][j].Scale(1/nsamperclust[i][j].GetEntries()) nsamperclust[i][j].SetLineColor(colors[i]) # nsamperpadpsa[i][j].SetStats(ROOT.kFALSE) # if nsamperpadpsa[i][j].GetEntries() > 0 : # nsamperpadpsa[i][j].Scale(1/nsamperpadpsa[i][j].GetEntries()) # nsamperpadpsa[i][j].SetLineColor(colors[i]) csig /= totentries cmsig /= totentries print "*************************************************************************************************" print "combined:" print "sigma:",csig,"+-",csige,"msigma:",cmsig,"+-",cmsige print "*************************************************************************************************" residupergain.SetPoint(i,gains[i],csig) residupergaincomb.SetPoint(i,gains[i],cmsig) if derr==1: residupergain.SetPointError(i,0.,math.sqrt(csige)) residupergaincomb.SetPoint(i,0.,math.sqrt(cmsige)) # }}} # hsampletot[i].SetLineColor(colors[i]) # slope,slopeerr=fitadc(hsampletot[i]) # adcvsgain.SetPoint(i,gains[i],slope) # adcvsgain.SetPointError(i,0,slopeerr) # }}} # snrmg.Add(snrPerGainSlice[i]) # snrmgmpv.Add(snrPerGainSlicempv[i]) # {{{ set tgraph titles snrPerGainSlice[0].GetXaxis().SetTitle("Drift Length Z (cm)") snrPerGainSlicempv[0].GetXaxis().SetTitle("Drift Length Z (cm)") snrPerGainSlice[0].GetYaxis().SetTitle("Mean SNR") snrPerGainSlicempv[0].GetYaxis().SetTitle("Most probable SNR") residuSlice[0].GetXaxis().SetTitle("Drift Length Z (cm)") residuSlice[0].GetYaxis().SetTitle("Resolution (#mum)") residuSlice[0].GetYaxis().SetTitleOffset(1.3) residuSlicecomb[0].GetXaxis().SetTitle("Drift Length Z (cm)") residuSlicecomb[0].GetYaxis().SetTitle("Resolution (#mum)") residuSlicecomb[0].GetYaxis().SetTitleOffset(1.3) # }}} #plotting # {{{ plotting for i in range(len(files)): # {{{ digis per pad # hdigiperpad[i].SetAxisRange(0.,1.,"Y") # canvas[0].cd() # if i == 0 : # hdigiperpad[i].Draw() # if prelim==1: # drawPrelim() # else: # hdigiperpad[i].Draw("same") # legend[0].AddEntry(hdigiperpad[i],nlegend[i]) # legend[0].Draw() # }}} # {{{ samples per pad # canvas[1].cd() # if i==0: # hnsamperpad[i].Draw() # if prelim==1: # drawPrelim() # else: # hnsamperpad[i].Draw("same") # legend[1].AddEntry(hnsamperpad[i],nlegend[i]) # legend[1].Draw() # }}} # {{{ signal to noise canvas[2].cd() hsnr[i].SetLineWidth(2) hsnr[i].SetTitle("") legend[2].AddEntry(hsnr[i],nlegend[i]) if i==0: hsnr[i].GetXaxis().SetTitle("Signal to Noise Ratio") hsnr[i].GetYaxis().SetTitle("Counts (normalized)") hsnr[i].GetYaxis().SetTitleOffset(1.3) hsnr[i].Draw() if prelim==1: drawPrelim() legend[2].Draw() else: hsnr[i].Draw("same") # }}} # {{{ samples per digi # canvas[3].cd() # legend[3].AddEntry(hsamperdigi[i],nlegend[i]) # legend[3].Draw() # if i==0: # hsamperdigi[i].Draw() # else: # hsamperdigi[i].Draw("same") # }}} # {{{ max cluster amplitude canvas[4].cd() legend[4].AddEntry(hClusterMax[i],nlegend[i]) legend[4].Draw() if i==0: hClusterMax[i].Draw() else: hClusterMax[i].Draw("same") # }}} # {{{ signal to noise mean vs gain snrPerGainSlice[i].SetMarkerColor(colors[i]) snrPerGainSlice[i].SetMarkerStyle(21) snrPerGainSlice[i].GetYaxis().SetRangeUser(0.,maxmeansnr*1.2) snrPerGainSlice[i].GetHistogram().SetBins(500,0.,maxmeansnr*1.2) legend[10].AddEntry(snrPerGainSlice[i],nlegend[i],"P") csnrslice.cd() if i==0: snrPerGainSlice[i].Draw("AP") legend[10].Draw() if prelim==1: drawPrelim() else: snrPerGainSlice[i].Draw("P") # }}} # {{{ signal to noise mpv vs gain snrPerGainSlicempv[i].SetMarkerColor(colors[i]) snrPerGainSlicempv[i].SetMarkerStyle(21) snrPerGainSlicempv[i].GetYaxis().SetRangeUser(0.,maxmpvsnr*1.2) snrPerGainSlicempv[i].GetHistogram().SetBins(500,0.,maxmpvsnr*1.2) legend[11].AddEntry(snrPerGainSlicempv[i],nlegend[i],"P") csnrmpvslice.cd() if i==0: snrPerGainSlicempv[i].Draw("AP") legend[11].Draw() if prelim==1: drawPrelim() else: snrPerGainSlicempv[i].Draw("P") # }}} # {{{ residuals narrow sliced or combined if len(files) > 1: legend[13].AddEntry(residuSlice[i],nlegend[i],"P") canvas[13].SetTitle("residuals vs z, narrow") legend[14].AddEntry(residuSlicecomb[i],nlegend[i],"P") canvas[14].SetTitle("residuals vs z, mean") else: legend[13].AddEntry(residuSlice[i],"Narrow Gaussian","P") legend[13].AddEntry(residuSlicecomb[i],"Weighted Mean","P") canvas[13].SetTitle("residuals vs z") residuSlice[i].SetMarkerColor(colors[i]) residuSlice[i].SetMarkerStyle(21) residuSlice[i].GetYaxis().SetRangeUser(0.,2000.) residuSlice[i].GetXaxis().SetRangeUser(0.,80.) residuSlicecomb[i].SetMarkerColor(colors[i]) residuSlicecomb[i].SetMarkerStyle(20) residuSlicecomb[i].GetYaxis().SetRangeUser(0.,2000.) residuSlicecomb[i].GetXaxis().SetRangeUser(0.,80.) canvas[13].cd() residuSlice[i].GetHistogram().SetBins(500,0.,70.) if i==0: residuSlice[i].Draw("AP") legend[13].Draw() if len(files) > 1 : canvas[14].cd() residuSlicecomb[i].Draw("AP") legend[14].Draw() else : residuSlicecomb[i].SetMarkerColor(ROOT.kRed) residuSlicecomb[i].Draw("P") if prelim==1: drawPrelim() canvas[13].cd() diffT.Draw("same") canvas[14].cd() diffT.Draw("same") else: residuSlice[i].Draw("P") if len(files) > 1 : canvas[14].cd() residuSlicecomb[i].Draw("P") # }}} # canvas[12].cd() canvas[12].Divide(3,2) # hresidu[i].Draw() # {{{ signal to noise slices if i==0: canvas[5].Divide(3,2) canvas[5].SetTitle("signal to noise sliced") for j in range(len(zSlices[i])): canvas[5].cd(j+1) nsnr[i][j].SetLineColor(colors[i]) if j==0: legend[5].AddEntry(nsnr[i][j],nlegend[i]) if i==0: nsnr[i][j].Draw() if prelim==1: drawPrelim() legend[5].Draw() else: nsnr[i][j].Draw("same") # }}} #{{{samples per pad if i==0: canvas[6].Divide(3,2) canvas[6].SetTitle("samples per pad") for j in range(len(zSlices[i])): canvas[6].cd(j+1) if j==0: legend[6].AddEntry(nsamperpad[i][j],nlegend[i]) if i==0: nsamperpad[i][j].Draw() if prelim==1: drawPrelim() legend[6].Draw() else: nsamperpad[i][j].Draw("same") # }}} # {{{digis per pad canvas[7].cd() if i==0: canvas[7].Divide(3,2) canvas[7].SetTitle("digis per pad") for j in range(len(zSlices[i])): canvas[7].cd(j+1) if j==0: legend[7].AddEntry(ndigiperpad[i][j],nlegend[i]) if i==0: ndigiperpad[i][j].Draw() legend[7].Draw() if prelim==1: drawPrelim() else: ndigiperpad[i][j].Draw("same") # }}} # {{{ samples per cluster canvas[15].cd() if i==0: canvas[15].Divide(3,2) for j in range(len(zSlices[i])): canvas[15].cd(j+1) if j==0: legend[15].AddEntry(nsamperclust[i][j],nlegend[i]) if i==0: nsamperclust[i][j].Draw() legend[15].Draw() if prelim==1: drawPrelim() else: nsamperclust[i][j].Draw("same") # }}} # {{{samples per digi # canvas[8].cd() # if i==0: # canvas[8].Divide(3,2) # canvas[8].SetTitle("samples per digi") # for j in range(len(zSlices[i])): # canvas[8].cd(j+1) # if j==0: # legend[8].AddEntry(nsamperdigi[i][j],nlegend[i]) # if i==0: # nsamperdigi[i][j].Draw() # legend[8].Draw() # if prelim==1: # drawPrelim() # else: # nsamperdigi[i][j].Draw("same") # }}} # {{{ sampels per pad after psa # if i==0: # canvas[9].Divide(3,2) # canvas[9].SetTitle("sampels per pad after psa") # for j in range(len(zSlices[i])): # canvas[9].cd(j+1) # if j==0: # legend[9].AddEntry(nsamperpadpsa[i][j],nlegend[i]) # if i==0: # nsamperpadpsa[i][j].Draw() # if prelim==1: # drawPrelim() # legend[9].Draw() # else: # nsamperpadpsa[i][j].Draw("same") # }}} # {{{ residuals slices for j in range(len(zSlices[i])): rescanvas[i].cd(j+1) #plotting of residual slice histos banane nresidu[i][j].GetXaxis().SetRangeUser(-2000.,2000.) nresidu[i][j].Draw() drawresdata(resdata[i][j]) if prelim==1: drawPrelim() canvas[12].cd(i+1) hresidu[i].Draw() legend[12].AddEntry(hresidu[i],nlegend[i]) if prelim==1: drawPrelim() # }}} canvases=[] legends=[] canvascounter=0 fcounter=0 hcounter=0 for hd in hists: hcounter=0 for h in hd: if fcounter==0: canvases.append(ROOT.TCanvas("c"+str(canvascounter+20),h)) legends.append(ROOT.TLegend(0.55,.9-len(nlegend)*0.05,0.9,.9,"","NDC")) legends[canvascounter].SetFillColor(0) canvascounter+=1 # canvases[hcounter].cd() print "h:",h, "hcounter:",hcounter if len(hd[h])>1 : if fcounter==0: canvases[hcounter].Divide(3,2) print "len:",len(hd[h]) for i in range(len(hd[h])): print "hcounter:",hcounter,"i:",i,"fcounter:",fcounter, canvases[hcounter].GetName() canvases[hcounter].cd(i+1) hd[h][i].SetLineColor(fcounter+1) if hd[h][i].GetEntries()>0: hd[h][i].Scale(1/hd[h][i].GetEntries()) if i==0: legends[hcounter].AddEntry(hd[h][i],nlegend[fcounter]) if fcounter==0: print "first draw" hd[h][i].Draw() legends[hcounter].Draw() else: hd[h][i].Draw("same") hcounter+=1 # nix=raw_input("bla") fcounter+=1 tempc = ROOT.TCanvas("temp","temp") nresidu[0][0].Draw() drawresdata(resdata[0][0]) if prelim == 1: drawPrelim() csnr = ROOT.TCanvas("csnr","snrvsgain") csnr.cd() snrPerGain.SetMarkerStyle(21) snrPerGain.GetXaxis().SetTitle("Gain") snrPerGain.GetYaxis().SetTitle("Mean SNR") snrPerGain.Draw("AP") if prelim==1: drawPrelim() csnr.SaveAs(pdir+"/snr_mean_vs_gain.pdf") csnrmpv = ROOT.TCanvas("csnrmpv","snr mpv") csnrmpv.cd() snrPerGainmpv.SetMarkerStyle(21) snrPerGainmpv.GetXaxis().SetTitle("Gain") snrPerGainmpv.GetYaxis().SetTitle("Most probable SNR") snrPerGainmpv.Draw("AP") if prelim==1: drawPrelim() csnrmpv.SaveAs(pdir+"/snr_mpv_vs_gain.pdf") cresisgain = ROOT.TCanvas("cresislice","residuals all vs gain") resleg = ROOT.TLegend(0.75,.9-2*0.05,0.9,.9,"","NDC") resleg.SetFillColor(0) resleg.AddEntry(residupergain,"Narrow Gaussian","P") resleg.AddEntry(residupergaincomb,"Mean","P") cresisgain.cd() residupergain.SetMarkerStyle(21) residupergain.GetXaxis().SetTitle("Gain") residupergain.GetYaxis().SetTitle("Resolution") residupergain.GetYaxis().SetRangeUser(0.,1000.) residupergain.Draw("AP") residupergaincomb.SetMarkerStyle(20) residupergaincomb.Draw("P") resleg.Draw() if prelim==1: drawPrelim() cresisgain.SaveAs(pdir+"/residuals_all_vs_gain.pdf") #csnrmpvslice.cd() #snrmgmpv.Draw("A") #csnrslice.cd() #snrmg.Draw("A") for i in range(len(canvas)): canvas[i].SaveAs(pdir+"/"+cnames[i]+".pdf") for i in range(len(files)): rescanvas[i].SaveAs(pdir+"/residual_slices_gain"+str(gains[i])+"_"+str(i)+".pdf") for c in canvases: c.SaveAs(pdir+"/"+c.GetTitle()+".pdf") tempc.SaveAs(pdir+"/temp.pdf") nix = raw_input("this is the end")