import sys, os, copy pandapath = os.environ.get('PANDAPATH') sys.path.append(pandapath + '/macro/tpc/FOPI/python/argparse-1.2.1') sys.path.append(pandapath + '/macro/tpc/FOPI/mberger') import argparse, ROOT, anaFile from functions import set_palette, openTree, thisIsTheEnd, getSlices, fitres import random def plotPackage(canv,graphs,hists,slices,folder="none"): canv[0].Divide(len(graphs),2) canv[1].Divide(len(graphs)) rfile.mkdir(folder) rfile.cd(folder) for i in range(len(graphs)): canv[0].cd(i+1) hists[i].Draw() canv[0].cd(i+1+len(graphs)) graphs[i]['2d'].Draw("surf1") graphs[i]['2d'].GetXaxis().SetTitle("Z (cm)") graphs[i]['2d'].GetYaxis().SetTitle("Theta (deg)") graphs[i]['2d'].GetYaxis().SetRangeUser(20,160) graphs[i]['2d'].GetYaxis().SetRangeUser(20,160) graphs[i]['2d'].GetYaxis().SetRangeUser(20,160) graphs[i]['2d'].Write() canv[1].cd(i+1) graphs[i]['Z int'].SetMarkerStyle(20) graphs[i]['Z int'].Draw("AP") graphs[i]['Z int'].Write() rfile.mkdir(folder+"/slices") rfile.cd(folder+"/slices") for i in range(len(slices)): for j in range(len(slices[i])): slices[i][j].Write() canv[0].Update() canv[1].Update() rfile.cd("canvas") canv[0].Write() canv[1].Write() # c1.cd(1) # h.Draw() # c1.cd(4) # resXGraph['2d'].Draw("surf1") # resXGraph['2d'].GetXaxis().SetTitle("Z (cm)") # resXGraph['2d'].GetYaxis().SetTitle("Theta (deg)") # c1.cd(2) # hresY.Draw() # c1.cd(5) # resYGraph['2d'].Draw("surf1") # c1.cd(3) # hresZ.Draw() # c1.cd(6) # resZGraph['2d'].Draw("surf1") # resXGraph['2d'].GetYaxis().SetRangeUser(20,160) # resYGraph['2d'].GetYaxis().SetRangeUser(20,160) # resZGraph['2d'].GetYaxis().SetRangeUser(20,160) # c1.Update() # rfile.cd() # c1.Write() # resXGraph['2d'].Write() # resYGraph['2d'].Write() # resZGraph['2d'].Write() # hresX.Write() # hresY.Write() # hresZ.Write() # c1z.Divide(3,1) # c1z.cd(1) # resXGraph['Z int'].Draw("AP") # c1z.cd(2) # resYGraph['Z int'].Draw("AP") # c1z.cd(3) # resZGraph['Z int'].Draw("AP") # resXGraph['Z int'].Write() # resYGraph['Z int'].Write() # resZGraph['Z int'].Write() def getBinXY(h,i,j): nbins=h.GetNbinsZ() minbin=h.GetZaxis().GetXmin() maxbin=h.GetZaxis().GetXmax() retH=ROOT.TH1D(h.GetName()+"bin_{0}_{1}".format(i,j),h.GetName()+"bin_{0}_{1}".format(i,j),nbins,minbin,maxbin) for xbin in range(i-1,i+2): for ybin in range(j-1,j+2): for ibi in range(nbins): val=h.GetBinContent(xbin,ybin,ibi) if val!=0: retH.SetBinContent(ibi,retH.GetBinContent(ibi)+val) return retH def getAllBins(h,slices): rfile.mkdir(h.GetName()) rfile.cd(h.GetName()) for i in range(1,h.GetNbinsX()-1): slices.append([]) rfile.mkdir('slices2D/{1}/zbin_{0}'.format(i,h.GetName())) rfile.cd('slices2D/{1}/zbin_{0}'.format(i,h.GetName())) for j in range(1,h.GetNbinsY()-1): slices[-1].append(getBinXY(h,i,j)) slices[-1][-1].Write() rfile.cd('/') rfile.cd(h.GetName()) def fitSlices(slices,binwidthX,binwidthy,binmin=0,name="graph",fitfunc="gaus",parnum=2): graph={} graph['2d']=ROOT.TGraph2D() graph['2d'].SetName(name) graph['Z int']=ROOT.TGraphErrors() graph['Z int'].SetName(name+"ZInt") graph['Z']=[] graph['theta']=[] slicesZ=[] for xbin in range(len(slices)): graph['Z'].append(ROOT.TGraph()) slicesZ.append(slices[xbin][0].Clone(slices[xbin][0].GetName()+"_Z")) for ybin in range(len(slices[xbin])): slicesZ[xbin].Add(slices[xbin][ybin],1) if xbin==0: graph['theta'].append(ROOT.TGraph()) h=slices[xbin][ybin] if h.GetEntries()!=0: ffunc=ROOT.TF1(fitfunc,fitfunc) h.Fit(ffunc,"QN") value_of_interest=ffunc.GetParameter(parnum) graph['2d'].SetPoint(graph['2d'].GetN(),binmin+binwidthX/2.+xbin*binwidthX,binwidthy/2.+ybin*binwidthy,value_of_interest) data=fitres(slicesZ[xbin]) graph['Z int'].SetPoint(graph['Z int'].GetN(),binmin+binwidthX/2.+xbin*binwidthX,data['csig']) graph['Z int'].SetPointError(graph['Z int'].GetN()-1,data['csigerr']) for g in graph: if g=='Z' or g=='theta': for theg in graph[g]: if theg.GetN()==0: theg.SetPoint(1,-1000,-1000) else: if graph[g].GetN()==0: if g=='2d': graph[g].SetPoint(1,-1000,-1000,-1000) else: graph[g].SetPoint(1,-1000,-1000) return graph, slicesZ parser = argparse.ArgumentParser(description='plot the residual width as function of z and theta') parser.add_argument('recofile', help='the reco file', type=str,nargs='+') parser.add_argument('--list',help='the recofile is actually a list of files',action='store_true') parser.add_argument('--runs', help='if data, the runs to use (%(default)s)', type=str, default='3888') parser.add_argument('--pattern', help='if data, the file name patters ('')', type=str, default='') parser.add_argument('--events', type=int, default=-1,help="number of events to process ( all )") parser.add_argument('--rout',help='root output file name (%(default)s)',type=str,default="binSlices.root") parser.add_argument('--rangeErr',help='ranges for error hists XYZ min, XYZ max, A min, A max',nargs=4,type=float,default=[0,600,0,600]) parser.add_argument('--fopi',help='use fopi data',action='store_true') args = parser.parse_args() ROOT.gROOT.ProcessLine(".x rootlogon.C") ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') ROOT.gROOT.LoadMacro("stlPYROOT.h+") set_palette("palette",99) trees=[] files=[] if args.list: for line in open(args.recofile[0]): line=line.replace('\n','') words=line.split(';') files.append(words[0]) else: files=args.recofile nbinsz=18 nbinst=45 zmin=0 zmax=75 if args.fopi: zmin=-50 zmax= 160 c1 =ROOT.TCanvas("CresXYZ","Res XYZ" ,2000 ,1000) c1z =ROOT.TCanvas("CresXYZ_Z","Res XYZ Z" ,2000 , 500) c2 =ROOT.TCanvas("CresA","Res A" ,2000 ,1000) c2z =ROOT.TCanvas("CresA_Z","Res A Z" ,2000 , 500) c3 =ROOT.TCanvas("CresUV","Res UV" ,2000 ,1000) c3z =ROOT.TCanvas("CresUV_Z","Res UV" ,2000 , 500) c4 =ROOT.TCanvas("CpullXYZ","Pull XYZ" ,2000 ,1000) c4z =ROOT.TCanvas("CpullXYZ_Z","Pull XYZ Z" ,2000 , 500) c5 =ROOT.TCanvas("CpullA","Pull A" ,2000 ,1000) c5z =ROOT.TCanvas("CpullA_Z","Pull A Z" ,2000 , 500) c6 =ROOT.TCanvas("CpullUV","Pull UV" ,2000 ,1000) c6z =ROOT.TCanvas("CpullUV_Z","Pull UV Z" ,2000 , 500) c7 =ROOT.TCanvas("CerrXYZ","Error XYZ" ,2000 ,1000) c7z =ROOT.TCanvas("CerrXYZ_Z","Error XYZ Z" ,2000 , 500) c8 =ROOT.TCanvas("CerrA","Error A" ,2000 ,1000) c8z =ROOT.TCanvas("CerrA_Z","Error A Z" ,2000 , 500) c9 =ROOT.TCanvas("CerrUV","Error UV" ,2000 ,1000) c9z =ROOT.TCanvas("CerrUV_Z","Error UV Z" ,2000 , 500) c10 =ROOT.TCanvas("CcluSizeXYZ","Cluster Size XYZ" ,2000 ,1000) c10z=ROOT.TCanvas("CcluSizeXYZ_Z","Cluster Size XYZ Z",2000 , 500) c11 =ROOT.TCanvas("CcluSizeA","Cluster Size A" ,2000 ,1000) c11z=ROOT.TCanvas("CcluSizeA_Z","Cluster Size A Z" ,2000 , 500) if args.recofile[0].find("hist")!=-1: anafile=anaFile.anaFile(args.recofile[0]) hresX = anafile.getHist("hSurfresX") hresY = anafile.getHist("hSurfresY") hresZ = anafile.getHist("hSurfresZ") hresA1 = anafile.getHist("hSurfresA1") hresA2 = anafile.getHist("hSurfresA2") hresU = anafile.getHist("hSurfresU") hresV = anafile.getHist("hSurfresV") hpullX = anafile.getHist("hSurfpullX") hpullY = anafile.getHist("hSurfpullY") hpullZ = anafile.getHist("hSurfpullZ") hpullA1 = anafile.getHist("hSurfpullA1") hpullA2 = anafile.getHist("hSurfpullA2") hpullU = anafile.getHist("hSurfpullU") hpullV = anafile.getHist("hSurfpullV") herrX = anafile.getHist("hSurferrX") herrY = anafile.getHist("hSurferrY") herrZ = anafile.getHist("hSurferrZ") herrA1 = anafile.getHist("hSurferrA1") herrA2 = anafile.getHist("hSurferrA2") herrU = anafile.getHist("hSurferrU") herrV = anafile.getHist("hSurferrV") hclsizeX = anafile.getHist("hSurfclsizeX") hclsizeY = anafile.getHist("hSurfclsizeY") hclsizeZ = anafile.getHist("hSurfclsizeZ") hclsizeA1 = anafile.getHist("hSurfclsizeA1") hclsizeA2 = anafile.getHist("hSurfclsizeA2") else: for itree in files: trees.append(openTree(itree,args.runs,args.pattern)) trees[-1].SetBranchStatus("*", 0) if args.fopi: trees[-1].SetBranchStatus("Particles.*",1) else: trees[-1].SetBranchStatus("CutCosmics.*", 1) hresX=ROOT.TH3D("hresX","residual X;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hresY=ROOT.TH3D("hresY","residual Y;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hresZ=ROOT.TH3D("hresZ","residual Z;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hresA1=ROOT.TH3D("hresA1","residual A1;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hresA2=ROOT.TH3D("hresA2","residual A2;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hresU=ROOT.TH3D("hresU","residual U;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hresV=ROOT.TH3D("hresV","residual V;Z (cm);Theta (deg);Res (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,-1500,1500) hpullX=ROOT.TH3D("hpullX","Pull X;Z(cm); #theta (#circ);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) hpullY=ROOT.TH3D("hpullY","Pull Y;Z(cm); #theta (#circ);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) hpullZ=ROOT.TH3D("hpullZ","Pull Z;Z(cm); #theta (#circ);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) hpullA1=ROOT.TH3D("hpullA1","pull A1;Z (cm);Theta (deg);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) hpullA2=ROOT.TH3D("hpullA2","pull A2;Z (cm);Theta (deg);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) hpullU=ROOT.TH3D("hpullU","pull U;Z (cm);Theta (deg);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) hpullV=ROOT.TH3D("hpullV","pull V;Z (cm);Theta (deg);Pull",nbinsz,zmin,zmax,nbinst,0,180,100,-10,10) herrX=ROOT.TH3D("herrX","Error X;Z(cm); #theta (#circ);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[0],args.rangeErr[1]) herrY=ROOT.TH3D("herrY","Error Y;Z(cm); #theta (#circ);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[0],args.rangeErr[1]) herrZ=ROOT.TH3D("herrZ","Error Z;Z(cm); #theta (#circ);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[0],args.rangeErr[1]) herrA1=ROOT.TH3D("herrA1","Error A1;Z (cm);Theta (deg);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[2],args.rangeErr[3]) herrA2=ROOT.TH3D("herrA2","Error A2;Z (cm);Theta (deg);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[2],args.rangeErr[3]) herrU=ROOT.TH3D("herrU","Error U;Z (cm);Theta (deg);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[2],args.rangeErr[3]) herrV=ROOT.TH3D("herrV","Error V;Z (cm);Theta (deg);Error (#mu m)",nbinsz,zmin,zmax,nbinst,0,180,100,args.rangeErr[2],args.rangeErr[3]) hclsizeX=ROOT.TH3D("hclsizeX","Cluster Size X;Z(cm); #theta (#circ);Size (cm)",nbinsz,zmin,zmax,nbinst,0,180,200,0,2) hclsizeY=ROOT.TH3D("hclsizeY","Cluster Size Y;Z(cm); #theta (#circ);Size (cm)",nbinsz,zmin,zmax,nbinst,0,180,200,0,2) hclsizeZ=ROOT.TH3D("hclsizeZ","Cluster Size Z;Z(cm); #theta (#circ);Size (cm)",nbinsz,zmin,zmax,nbinst,0,180,200,0,2) hclsizeA1=ROOT.TH3D("hclsizeA1","Cluster Size A1;Z (cm);Theta (deg);Size (cm)",nbinsz,zmin,zmax,nbinst,0,180,200,0,2) hclsizeA2=ROOT.TH3D("hclsizeA2","Cluster Size A2;Z (cm);Theta (deg);Size (cm)",nbinsz,zmin,zmax,nbinst,0,180,200,0,2) iev=-1 for e in trees[-1]: iev+=1 if iev>args.events and args.events!=-1: break if iev%100==0: print iev, if iev%1000==0: print "" print 'at event:',iev if args.fopi: branch=e.Particles else: branch=e.CutCosmics for tr in branch: if not args.fopi: theta=tr.GetTheta() ncluster=tr.nhits() for icl in range(ncluster): resXYZ=tr.GetRes(icl) posXYZ=tr.GetHitPositionXYZ(icl) if args.fopi: theta=tr.getTrackMom().Theta()*ROOT.TMath.RadToDeg() hresX.Fill(posXYZ.Z(),theta,resXYZ.X()*10000) #residual in mum hresY.Fill(posXYZ.Z(),theta,resXYZ.Y()*10000) hresZ.Fill(posXYZ.Z(),theta,resXYZ.Z()*10000) resA=tr.GetResA(icl) hresA1.Fill(posXYZ.Z(),theta,resA.Y()*10000) hresA2.Fill(posXYZ.Z(),theta,resA.Z()*10000) resUV=tr.GetResUV(icl) hresU.Fill(posXYZ.Z(),theta,resUV.X()*10000) hresV.Fill(posXYZ.Z(),theta,resUV.Y()*10000) errXYZ=tr.GetSigXYZ(icl) herrX.Fill(posXYZ.Z(),theta,errXYZ.X()*10000) herrY.Fill(posXYZ.Z(),theta,errXYZ.Y()*10000) herrZ.Fill(posXYZ.Z(),theta,errXYZ.Z()*10000) errA=tr.GetSigA(icl) herrA1.Fill(posXYZ.Z(),theta,errA.Y()*10000) herrA2.Fill(posXYZ.Z(),theta,errA.Z()*10000) errUV=tr.GetSigUV(icl) herrU.Fill(posXYZ.Z(),theta,errUV.X()*10000) herrV.Fill(posXYZ.Z(),theta,errUV.Y()*10000) hpullX.Fill(posXYZ.Z(),theta,resXYZ.X()/errXYZ.X()) hpullY.Fill(posXYZ.Z(),theta,resXYZ.Y()/errXYZ.Y()) hpullZ.Fill(posXYZ.Z(),theta,resXYZ.Z()/errXYZ.Z()) hpullA1.Fill(posXYZ.Z(),theta,resA.Y()/errA.Y()) hpullA2.Fill(posXYZ.Z(),theta,resA.Z()/errA.Z()) hpullU.Fill(posXYZ.Z(),theta,resUV.X()/errUV.X()) hpullV.Fill(posXYZ.Z(),theta,resUV.Y()/errUV.Y()) clsizeXYZ=tr.GetFullClusterSize(icl) hclsizeX.Fill(posXYZ.Z(),theta,clsizeXYZ.X()) hclsizeY.Fill(posXYZ.Z(),theta,clsizeXYZ.Y()) hclsizeZ.Fill(posXYZ.Z(),theta,clsizeXYZ.Z()) clsizeA=tr.GetFullClusterSizeAxis(icl) hclsizeA1.Fill(posXYZ.Z(),theta,clsizeA.Y()) hclsizeA2.Fill(posXYZ.Z(),theta,clsizeA.Z()) binwidthZ=(zmax-zmin)/nbinsz print zmin,zmax,binwidthZ,zmax-zmin,nbinsz binwidthT=180./nbinst rfile=ROOT.TFile(args.rout,'recreate') rfile.mkdir("slices2D") rfile.mkdir("canvas") resXslices=[] resYslices=[] resZslices=[] print "getting all residual XYZ bins" getAllBins(hresX, resXslices) getAllBins(hresY, resYslices) getAllBins(hresZ, resZslices) print "Fitting all residual XYZ bins" resXGraph,ZslicesresX=fitSlices(resXslices,binwidthZ,binwidthT,zmin,"g"+hresX.GetName()) resYGraph,ZslicesresY=fitSlices(resYslices,binwidthZ,binwidthT,zmin,"g"+hresY.GetName()) resZGraph,ZslicesresZ=fitSlices(resZslices,binwidthZ,binwidthT,zmin,"g"+hresZ.GetName()) resA1slices=[] resA2slices=[] print "getting all residual A bins" getAllBins(hresA1,resA1slices) getAllBins(hresA2,resA2slices) print "Fitting all residual A bins" resA1Graph,zSlicesresA1=fitSlices(resA1slices,binwidthZ,binwidthT,zmin,"g"+hresA1.GetName()) resA2Graph,zSlicesresA2=fitSlices(resA2slices,binwidthZ,binwidthT,zmin,"g"+hresA2.GetName()) resUslices=[] resVslices=[] print "getting all residual UV bins" getAllBins(hresU,resUslices) getAllBins(hresV,resVslices) print "Fitting all residual UV bins" resUGraph,zSlicesresU=fitSlices(resUslices,binwidthZ,binwidthT,zmin,"g"+hresU.GetName()) resVGraph,zSlicesresV=fitSlices(resVslices,binwidthZ,binwidthT,zmin,"g"+hresV.GetName()) #------------------------------------------------------------------------------------------------------------------------------------------------- errXslices=[] errYslices=[] errZslices=[] print "getting all error XYZ bins" getAllBins(herrX, errXslices) getAllBins(herrY, errYslices) getAllBins(herrZ, errZslices) print "Fitting all error XYZ bins" errXGraph,ZsliceserrX=fitSlices(errXslices,binwidthZ,binwidthT,zmin,"g"+herrX.GetName(),"landau",1) errYGraph,ZsliceserrY=fitSlices(errYslices,binwidthZ,binwidthT,zmin,"g"+herrY.GetName(),"landau",1) errZGraph,ZsliceserrZ=fitSlices(errZslices,binwidthZ,binwidthT,zmin,"g"+herrZ.GetName(),"landau",1) errA1slices=[] errA2slices=[] print "getting all error A bins" getAllBins(herrA1,errA1slices) getAllBins(herrA2,errA2slices) print "Fitting all error A bins" errA1Graph,ZsliceserrA1=fitSlices(errA1slices,binwidthZ,binwidthT,zmin,"g"+herrA1.GetName(),"landau",1) errA2Graph,ZsliceserrA2=fitSlices(errA2slices,binwidthZ,binwidthT,zmin,"g"+herrA2.GetName(),"landau",1) errUslices=[] errVslices=[] print "getting all error UV bins" getAllBins(herrU,errUslices) getAllBins(herrV,errVslices) print "Fitting all error UV bins" errUGraph,ZsliceserrU=fitSlices(errUslices,binwidthZ,binwidthT,zmin,"g"+herrU.GetName(),"landau",1) errVGraph,ZsliceserrV=fitSlices(errVslices,binwidthZ,binwidthT,zmin,"g"+herrV.GetName(),"landau",1) #------------------------------------------------------------------------------------------------------------------------------------------------- pullXslices=[] pullYslices=[] pullZslices=[] print "getting all pull XYZ bins" getAllBins(hpullX, pullXslices) getAllBins(hpullY, pullYslices) getAllBins(hpullZ, pullZslices) print "Fitting all pull XYZ bins" pullXGraph,ZslicespullX=fitSlices(pullXslices,binwidthZ,binwidthT,zmin,"g"+hpullX.GetName()) pullYGraph,ZslicespullY=fitSlices(pullYslices,binwidthZ,binwidthT,zmin,"g"+hpullY.GetName()) pullZGraph,ZslicespullZ=fitSlices(pullZslices,binwidthZ,binwidthT,zmin,"g"+hpullZ.GetName()) pullA1slices=[] pullA2slices=[] print "getting all pull A bins" getAllBins(hpullA1,pullA1slices) getAllBins(hpullA2,pullA2slices) print "Fitting all pull A bins" pullA1Graph,ZslicespullA1=fitSlices(pullA1slices,binwidthZ,binwidthT,zmin,"g"+hpullA1.GetName()) pullA2Graph,ZslicespullA2=fitSlices(pullA2slices,binwidthZ,binwidthT,zmin,"g"+hpullA2.GetName()) pullUslices=[] pullVslices=[] print "getting all pull UV bins" getAllBins(hpullU,pullUslices) getAllBins(hpullV,pullVslices) print "Fitting all pull UV bins" pullUGraph,ZslicespullU=fitSlices(pullUslices,binwidthZ,binwidthT,zmin,"g"+hpullU.GetName()) pullVGraph,ZslicespullV=fitSlices(pullVslices,binwidthZ,binwidthT,zmin,"g"+hpullV.GetName()) #------------------------------------------------------------------------------------------------------------------------------------------------- clsizeXslices=[] clsizeYslices=[] clsizeZslices=[] print "getting all clsize XYZ bins" getAllBins(hclsizeX, clsizeXslices) getAllBins(hclsizeY, clsizeYslices) getAllBins(hclsizeZ, clsizeZslices) print "Fitting all clsize XYZ bins" clsizeXGraph,ZslicesClsizeX=fitSlices(clsizeXslices,binwidthZ,binwidthT,zmin,"g"+hclsizeX.GetName(),"gaus",1) clsizeYGraph,ZslicesClsizeY=fitSlices(clsizeYslices,binwidthZ,binwidthT,zmin,"g"+hclsizeY.GetName(),"gaus",1) clsizeZGraph,ZslicesClsizeZ=fitSlices(clsizeZslices,binwidthZ,binwidthT,zmin,"g"+hclsizeZ.GetName(),"gaus",1) clsizeA1slices=[] clsizeA2slices=[] print "getting all clsize A bins" getAllBins(hclsizeA1,clsizeA1slices) getAllBins(hclsizeA2,clsizeA2slices) print "Fitting all clsize A bins" clsizeA1Graph,ZslicesClsizeA1=fitSlices(clsizeA1slices,binwidthZ,binwidthT,zmin,"g"+hclsizeA1.GetName()) clsizeA2Graph,ZslicesClsizeA2=fitSlices(clsizeA2slices,binwidthZ,binwidthT,zmin,"g"+hclsizeA2.GetName()) #------------------------------------------------------------------------------------------------------------------------------------------------- #------------------------------------------------------------------------------------------------------------------------------------------------- plotPackage([c1,c1z],[resXGraph,resYGraph,resZGraph], [hresX,hresY,hresZ],[ZslicesresX, ZslicesresY,ZslicesresZ],"resXYZ") plotPackage([c2,c2z],[resA1Graph,resA2Graph],[hresA1,hresA2],[zSlicesresA1,zSlicesresA2],"resA") plotPackage([c3,c3z],[resUGraph,resVGraph],[hresU,hresV],[zSlicesresU,zSlicesresV],"resUV") plotPackage([c4,c4z],[pullXGraph,pullYGraph,pullZGraph],[hpullX,hpullY,hpullZ],[ZslicespullX,ZslicespullY,ZslicespullZ],"pullXYZ") plotPackage([c5,c5z],[pullA1Graph,pullA2Graph],[hpullA1,hpullA2],[ZslicespullA1,ZslicespullA2],"pullA") plotPackage([c6,c6z],[pullUGraph,pullVGraph],[hpullU,hpullV],[ZslicespullU,ZslicespullV],"pullUV") #-------------------------------------------------- c7.Divide(3,2) c7.cd(1) herrX.Draw() c7.cd(4) errXGraph['2d'].Draw("surf1") errXGraph['2d'].GetXaxis().SetTitle("Z (cm)") errXGraph['2d'].GetYaxis().SetTitle("Theta (deg)") c7.cd(2) herrY.Draw() c7.cd(5) errYGraph['2d'].Draw("surf1") c7.cd(3) herrZ.Draw() c7.cd(6) errZGraph['2d'].Draw("surf1") errXGraph['2d'].GetYaxis().SetRangeUser(20,160) errYGraph['2d'].GetYaxis().SetRangeUser(20,160) errZGraph['2d'].GetYaxis().SetRangeUser(20,160) c7.Update() rfile.cd() c7.Write() errXGraph['2d'].Write() errYGraph['2d'].Write() errZGraph['2d'].Write() herrX.Write() herrY.Write() herrZ.Write() c8.Divide(2,2) c8.cd(1) herrA1.Draw() c8.cd(3) errA1Graph['2d'].Draw("surf1") c8.cd(2) herrA2.Draw() c8.cd(4) errA2Graph['2d'].Draw("surf1") errA1Graph['2d'].GetYaxis().SetRangeUser(20,160) errA2Graph['2d'].GetYaxis().SetRangeUser(20,160) c8.Update() c8.Write() herrA1.Write() herrA2.Write() errA1Graph['2d'].Write() errA2Graph['2d'].Write() c9.Divide(2,2) c9.cd(1) herrU.Draw() c9.cd(3) errUGraph['2d'].Draw("surf1") c9.cd(2) herrV.Draw() c9.cd(4) errVGraph['2d'].Draw("surf1") errUGraph['2d'].GetYaxis().SetRangeUser(20,160) errVGraph['2d'].GetYaxis().SetRangeUser(20,160) c9.Update() c9.Write() herrU.Write() herrV.Write() errUGraph['2d'].Write() errVGraph['2d'].Write() #-------------------------------------------------- c10.Divide(3,2) c10.cd(1) hclsizeX.Draw() c10.cd(4) clsizeXGraph['2d'].Draw("surf1") clsizeXGraph['2d'].GetXaxis().SetTitle("Z (cm)") clsizeXGraph['2d'].GetYaxis().SetTitle("Theta (deg)") c10.cd(2) hclsizeY.Draw() c10.cd(5) clsizeYGraph['2d'].Draw("surf1") c10.cd(3) hclsizeZ.Draw() c10.cd(6) clsizeZGraph['2d'].Draw("surf1") clsizeXGraph['2d'].GetYaxis().SetRangeUser(20,160) clsizeYGraph['2d'].GetYaxis().SetRangeUser(20,160) clsizeZGraph['2d'].GetYaxis().SetRangeUser(20,160) c10.Update() rfile.cd() c10.Write() clsizeXGraph['2d'].Write() clsizeYGraph['2d'].Write() clsizeZGraph['2d'].Write() hclsizeX.Write() hclsizeY.Write() hclsizeZ.Write() c11.Divide(2,2) c11.cd(1) hclsizeA1.Draw() c11.cd(3) clsizeA1Graph['2d'].Draw("surf1") c11.cd(2) hclsizeA2.Draw() c11.cd(4) clsizeA2Graph['2d'].Draw("surf1") clsizeA1Graph['2d'].GetYaxis().SetRangeUser(20,160) clsizeA2Graph['2d'].GetYaxis().SetRangeUser(20,160) c11.Update() c11.Write() hclsizeA1.Write() hclsizeA2.Write() clsizeA1Graph['2d'].Write() clsizeA2Graph['2d'].Write() thisIsTheEnd() rfile.Close()