import glob,math,sys,os,copy pandapath=os.environ.get('PANDAPATH') sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1') import argparse sys.path.append('$PANDAPATH/macro/tpc/FOPI/mberger') from anaFile import anaFile from functions import * def DrawPulls(c,h,g,opt): c.Divide(3,3) for i in range(3): c.cd(i+1) h[i].Draw('colz') c.cd(i+4) g[i]['Mean'].Draw(opt) c.cd(i+7) g[i]['RMS'].Draw(opt) def setGraph(g,title,col=1): colors=[] for i in range(0,10): colors.append(ROOT.kViolet-i) for i in range(1,11): colors.append(ROOT.kViolet+i) g.SetTitle(title) g.SetMarkerStyle(20) g.SetMarkerSize(1) g.SetMarkerColor(colors[col]) return def fillGraph(hist,graph,fnum=0,mini={},maxi={}): nbins=hist.GetNbinsX() projs=[] graph.append({}) graph[-1]['RMS']=ROOT.TGraphErrors() setGraph(graph[-1]['RMS'],'RMS',fnum+1) graph[-1]['Mean']=ROOT.TGraphErrors() setGraph(graph[-1]['Mean'],'Mean',fnum+1) pcounter=0 binwidth=args.binwidth i=binwidth while 1: projs.append(hist.ProjectionY(hist.GetName()+str(i),i-binwidth,i)) graph[-1]['RMS'].SetPoint(pcounter,i-binwidth+1,projs[-1].GetRMS()) graph[-1]['RMS'].SetPointError(pcounter,float(1./2.),projs[-1].GetRMSError()) graph[-1]['Mean'].SetPoint(pcounter,i-binwidth+1,projs[-1].GetMean()) graph[-1]['Mean'].SetPointError(pcounter,float(1/2.),projs[-1].GetMeanError()) graph[-1]['RMS'].SetName(hist.GetName()+'_RMS') graph[-1]['Mean'].SetName(hist.GetName()+'_Mean') mini['RMS']=min(mini.get('RMS',1000),projs[-1].GetRMS()) maxi['RMS']=max(maxi.get('RMS',0),projs[-1].GetRMS()) mini['Mean']=min(mini.get('Mean',1000),projs[-1].GetMean()) maxi['Mean']=max(maxi.get('Mean',0),projs[-1].GetMean()) pcounter+=1 i+=1 if (i>=nbins): break return projs def fitSlices(hist,graph,func=0,fnum=0,mini={},maxi={}): print 'starting fitslices for {0}, with func={1}'.format(hist.GetName(),func) if func==1: ffunc=ROOT.TF1("landau","landau",0,hist.GetRMS(2)*3) ffunc.SetParLimits(1,0,10000) nbins=hist.GetNbinsX() projs=[] parname=['Sigma ','Mean ','Const '] if func==1: parname=['Sigma ','MPV ','Const '] graph.append({}) graph[-1]['RMS']=ROOT.TGraphErrors() setGraph(graph[-1]['RMS'],parname[0]+' '+hist.GetTitle(),fnum+1) graph[-1]['Mean']=ROOT.TGraphErrors() setGraph(graph[-1]['Mean'],parname[1]+' '+hist.GetTitle(),fnum+1) graph[-1]['Const']=ROOT.TGraphErrors() setGraph(graph[-1]['Const'],parname[2]+' '+hist.GetTitle(),fnum+1) pcounter=0 binwidth=args.binwidth i=binwidth #allData[hist.GetTitle()]=[] while 1: projs.append(hist.ProjectionY(hist.GetName()+str(i),i-binwidth,i)) projs[-1].SetTitle(hist.GetTitle()+'_'+str(i)) #print 'func=',func if func==1: projs[-1].Fit(ffunc,"QREM+") projs[-1].SetLineColor(2) #print 'fitting landau' data={} data['nsig']=ffunc.GetParameter(2) data['nsigerr']=ffunc.GetParError(2) data['nmean']=ffunc.GetParameter(1) data['nmeanerr']=ffunc.GetParError(1) data['nconst']=ffunc.GetParameter(0) data['nconsterr']=ffunc.GetParError(0) else: data=fitres(projs[-1],debug) allData[projs[-1].GetTitle()]=data graph[-1]['RMS'].SetPoint(pcounter,i-int(binwidth/2.),data['nsig']) graph[-1]['RMS'].SetPointError(pcounter,float(1./2.),data['nsigerr']) graph[-1]['Mean'].SetPoint(pcounter,i-int(binwidth/2.),data['nmean']) graph[-1]['Mean'].SetPointError(pcounter,float(1/2.),data['nmeanerr']) graph[-1]['Const'].SetPoint(pcounter,i-int(binwidth/2.),data['nconst']) graph[-1]['Const'].SetPointError(pcounter,float(1/2.),data['nconsterr']) mini['RMS']=min(mini.get('RMS',1000),data['nsig']) maxi['RMS']=max(maxi.get('RMS',0),data['nsig']) mini['Mean']=min(mini.get('Mean',1000),data['nmean']) maxi['Mean']=max(maxi.get('Mean',0),data['nmean']) pcounter+=1 i+=1 if (i>=nbins): break graph[-1]['RMS'].SetName(hist.GetName()+'_Sigma') graph[-1]['Mean'].SetName(hist.GetName()+'_Mean') graph[-1]['Const'].SetName(hist.GetName()+'_Const') return projs def getSlices(h,projs,opt=0,func=0): graph=[] if opt==0: tmp=fitSlices(h,graph,func) elif opt==1: tmp=fillGraph(h,graph) for p in tmp: projs.append(copy.deepcopy(p)) allProj[h.GetTitle()]=(copy.deepcopy(projs)) allGraph.append(copy.deepcopy(graph[-1])) return graph[-1] def calcproj(res,ax,amp): #calculate projection if ax.Mag()>0: mag=ax.Mag() else: mag=1 proj=(ax*res)/mag # mag=math.sqrt(mag) #for test val=proj/mag #for calibration #val*=math.sqrt(amp) return val def calcproj2(ax,amp): val=ax.Mag() return val def calcproj3(res,ax,scale=1): if ax.Mag()>0: mag=ax.Mag() else: mag=1 proj=(ax*res)/mag #scale=math.sqrt(scale) val=proj/(mag*scale) return val def scalefunc(z,i,ax): val=1. if z>0: #val=calcxyzproj(Dl,Dl,Dt,ax)*math.sqrt(z) #sz=math.sqrt(z) #val=calcxyzproj(Dl*sz,Dl*sz,Dt*sz,ax,True) #val=calcxyzproj(Dl*sz,0,0,ax,True) #val+=calcxyzproj(0,Dl*sz,0,ax,True) #val+=calcxyzproj(0,0,Dt*sz,ax,True) geomMat=ROOT.TMatrixDSym(3) tax=ROOT.TVectorD(3) if ax.Mag()!=0: tax[0]=ax[0]/ax.Mag() tax[1]=ax[1]/ax.Mag() tax[2]=ax[2]/ax.Mag() else: return 1 geomMat[0][0]=3**2 geomMat[1][1]=2.598**2 #geomMat[2][2]=1.52867**2 geomMat[2][2]=1 geomMat[0][1]=geomMat[0][2]=geomMat[1][0]=geomMat[1][2]=geomMat[2][0]=geomMat[2][1]=0 geomMat=geomMat.InvertFast() geomMat*=10 val=geomMat.Similarity(tax) if val!=0: val=1./val else: val=1 #else: # print 'scalefunc z<0: i=',i,'z=',z ''' if i==1: val*=0.253/math.sqrt(z) elif i==2: val*=0.683/math.sqrt(z) ''' return abs(val) def calcxyzproj(x,y,z,ax,dodiff=False): val=1. if ax.Mag()>0: Diff=ROOT.TVector3(x,y,z) if dodiff: val=(ax*Diff)/ax.Mag() else: val=(ax*Diff) #else: #print ax.Mag(),ax #ax.Print() return val def calcAxProj(vec,ax): val=0 if ax.Mag()>0: val=(ax*vec)/ax.Mag() return val def calcprojT(res,phi): alpha=res.Phi()-phi longiR=res.Perp()*ROOT.TMath.Cos(alpha) perp1R=res.Perp()*ROOT.TMath.Sin(alpha) perp2R=res.Z() val=ROOT.TVector3(longiR,perp1R,perp2R) return val def getStartPars(hist): pars=[] f1=ROOT.TF1('f1','[0]',0,73) f2=ROOT.TF1('f2','[0]+[1]*sqrt(x)',10,73) f3=ROOT.TF1('f3','[0]+[1]*exp(-x^[2])',0,20) opt='QRN' if args.getParV: opt='R' hist.Fit(f1,opt) f2.SetParameter(0,f1.GetParameter(0)) hist.Fit(f2,opt) f3.SetParameter(0,f2.GetParameter(0)) f3.SetParLimits(1,0,100) f3.SetParLimits(2,0,10) hist.Fit(f3,opt) if args.getParV: print 'getting startpars for',hist.GetName() ctemp=ROOT.TCanvas() ctemp.cd() hist.Draw(graphopt) f1.SetLineColor(2) f2.SetLineColor(3) f3.SetLineColor(4) f1.SetLineWidth(1) f2.SetLineWidth(1) f3.SetLineWidth(1) f1.Draw('same') f2.Draw('same') f3.Draw('same') u=raw_input('done looking at start pars?') offset=(f2.GetParameter(0)+f3.GetParameter(0))/2. pars=[offset,f2.GetParameter(1),f3.GetParameter(1),f3.GetParameter(2)] return pars def calcCovScale(cov,res,amp): cov2=ROOT.TMatrixDSym(3,cov.GetMatrixArray()) cov_i=ROOT.TMatrixDSym(cov2) #print 'inverting cov' #cov_i.Print() cov_i.InvertFast() #cov_i.Print() #print 'inverted cov' #cov_i=ROOT.TMatrixDSym(ROOT.TMatrixDSym.kInverted,cov2) # if cov[0][0]==0 or cov[1][1]==0 or cov[2][2]==0: # cov.Print() residual=ROOT.TMatrixD(3,1) if res.Mag()==0: return -1,-1 for i in range(3): residual[i][0]=res[i]/res.Mag() residual_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,residual) if args.cutCov: temp2=ROOT.TMatrixD(cov_i,ROOT.TMatrixD.kMult,residual) #print 'using inverted' else: temp2=ROOT.TMatrixD(cov,ROOT.TMatrixD.kMult,residual) lam=ROOT.TMatrixD(residual_t,ROOT.TMatrixD.kMult,temp2) #print'orig lam:',lam.Print() if args.cutCov: lam=1./math.sqrt(abs(lam[0][0])) #print 'inverting' else: lam=math.sqrt(abs(lam[0][0])) num=res.Mag()/abs(lam) #print num, abs(lam) return num,abs(lam) #def Draw(c,h): # c.Divide(len(h),3) parser=argparse.ArgumentParser(description='calculate a parametrization for the cluster error') parser.add_argument('recofile',help='the reco file') parser.add_argument('--events',help='number of events to analyze',type=int, default=-1) parser.add_argument('--new',help='create the histo file anew',action='store_true') parser.add_argument('--basepath',help='the path to store the output files',type=str,default='./') parser.add_argument('--outfile',help='name for the temp hist file (default=clerrcal.root)',type=str,default='clerrcal.root') parser.add_argument('--phi',help='the phi angle used for this simulation in degree',type=float,default=-1) parser.add_argument('--fitfile',help='the file to write the fit information into',type=str,default="None") parser.add_argument('--batch',help='start in batchmode',action='store_true') parser.add_argument('--plot',help='plot pictures in the folder where the recofile is',action='store_true') parser.add_argument('--pfile',help='the file to plot all the stuff into',type=str,default='None') parser.add_argument('--genfit',help='use cluster to genfit track residual',action='store_true') parser.add_argument('--track',help='use mcpos to genfit track residual',action='store_true') parser.add_argument('--binwidth',help='width of bins',type=int,default=2) parser.add_argument('--getParV',help='be verbose with getPars',action='store_true') parser.add_argument('--sliceopt',help='0 to fit gauss, 1 to use RMS of slices',type=int,default=0) parser.add_argument('--showRMS',help='show RMS distribution for i',type=int,default=-1) parser.add_argument('--cutCov',help='use cut instead of proj for scaling calculation',action='store_true') parser.add_argument('--autorms',help='try to find rms automaticly',action='store_true') parser.add_argument('--rmsfile',help='the file with rms values (only if not autorms)',type=str,default='./rms.txt') parser.add_argument('--ax3Zax',help='cut on angle of axis 3 and z axis (inside window)',type=float,nargs=2,default=[0,180]) parser.add_argument('--clsize2D',help='cut on 2D clustersize (inside window)',type=int,nargs=2,default=[0,99999999]) parser.add_argument('--clsize3D',help='cut on 3D clustersize (inside window)',type=int,nargs=2,default=[0,99999999]) parser.add_argument('--title',help='the title',type=str,default='') parser.add_argument('--stderr',help='disable cuts the std err cant provide',action='store_true') parser.add_argument('--sort',help='dont sort the axes according to their orientation',action='store_true') parser.add_argument('--hlp',help='print help',action='store_true') args=parser.parse_args() global debug debug=0 if args.batch: sys.argv.append( '-b-' ) import ROOT if args.batch: ROOT.gROOT.SetBatch(True) ROOT.gROOT.ProcessLine(".x rootlogon.C") ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') ROOT.gStyle.SetOptFit(1111) set_palette() if args.hlp: parser.print_help() exit(0) graphopt='AP' #if args.sliceopt==1: # graphopt='AP' c2=ROOT.TCanvas('c2','FuncX') cfuncy=ROOT.TCanvas('cfuncy','FuncY') cfuncz=ROOT.TCanvas('cfuncZ','FuncZ') cerr=ROOT.TCanvas('cerr','Error',1200,1000) cerr.Divide(3,3) cerr2=ROOT.TCanvas('cerr2','Error 2',1200,1000) cerr2.Divide(3,3) cerrT=ROOT.TCanvas('cerrT','Error T',1200,1000) cerrT2=ROOT.TCanvas('cerrT2','Error T2',1200,1000) cerrA=ROOT.TCanvas('cerrA','Error A',1200,1000) cerrA2=ROOT.TCanvas('cerrA2','Error A2',1200,1000) cerrA3=ROOT.TCanvas('cerrA3','Error A3',1200,1000) cerrNum=ROOT.TCanvas('cerrNum','Scaling factor',1200,1000) cScale=ROOT.TCanvas('cScale','Diff Scaling',1200,1000) cErrUV=ROOT.TCanvas('cErrUV','Error UV',1200,1000) cPullUV=ROOT.TCanvas('cPullUV','Pull UV',1200,1000) cprojA1=ROOT.TCanvas('cprojA1','Proj A1',1200,1000) cprojA2=ROOT.TCanvas('cprojA2','Proj A2',1200,1000) cprojA3=ROOT.TCanvas('cprojA3','Proj A3',1200,1000) cprojA1_3=ROOT.TCanvas('cprojA1_3','Proj A1 3',1200,1000) cprojA2_3=ROOT.TCanvas('cprojA2_3','Proj A2 3',1200,1000) cprojA3_3=ROOT.TCanvas('cprojA3_3','Proj A3 3',1200,1000) cprojNum=ROOT.TCanvas('cprojNum','Proj Num',1200,1000) cprojPullU=ROOT.TCanvas('cprojPullU','Proj Pull U',1200,1000) cprojPullV=ROOT.TCanvas('cprojPullV','Proj Pull V',1200,1000) cerrC=ROOT.TCanvas('cerrC','Err XYZ Cov cut',1200,1000) cPullC=ROOT.TCanvas('cPullC','Pull XYZ Cov cut',1200,1000) cangle=ROOT.TCanvas('cangle','Angles',1200,1000) cangle2=ROOT.TCanvas('cangle2','Angles 2',1200,1000) cres=ROOT.TCanvas('cres','residuals',1200,1000) cerr3=ROOT.TCanvas('cerr3','errors',1200,1000) if args.pfile !='None': cproj=ROOT.TCanvas('cproj','Projection') dx=0.3 dy=0.2598 dz=0.152867 Dl=0.0259143 Dt=0.0197879 Xaxis=ROOT.TVector3(1,0,0) Yaxis=ROOT.TVector3(0,1,0) Zaxis=ROOT.TVector3(0,0,1) htheta=[] hphi=[] hdangle=[] hdanglex=[] hdangley=[] hdanglez=[] if args.new: rfile=ROOT.TFile(args.recofile,'read') print args.recofile tree=rfile.Get('cbmsim') if tree==None: print 'no tree found' exit() tree.SetBranchStatus('*',0) #tree.SetBranchStatus('TrackFitStat_0.*',1) tree.SetBranchStatus('CutCosmics.*',1) #tree.SetBranchStatus('TpcSPHit.*',1) #tree.SetBranchStatus('TpcCluster.*',1) i=-1 rmsX=[] rmsY=[] rmsZ=[] for line in open(args.rmsfile,'r'): i+=1 words=line.split() rmsX.append(float(words[0])) rmsY.append(float(words[1])) rmsZ.append(float(words[2])) print 'using rms:',i,'X=',rmsX[-1],'Y=',rmsY[-1],'Z=',rmsZ[-1] hErrXAmpVsZ=ROOT.TH2D('hErrXrAmpVsZ','mcresX/ferrX vs Z '+args.title,75,0,75,150,-rmsX[0],rmsX[0]) hErrYAmpVsZ=ROOT.TH2D('hErrYrAmpVsZ','mcresY/ferrY vs Z '+args.title,75,0,75,150,-rmsY[0],rmsY[0]) hErrZAmpVsZ=ROOT.TH2D('hErrZrAmpVsZ','mcresZ/ferrZ vs Z '+args.title,75,0,75,150,-rmsZ[0],rmsZ[0]) hErrXAmpVsZ2=ROOT.TH2D('hErrXAmpVsZ2','ferrX vs Z '+args.title,75,0,75,150,0,rmsX[1]) hErrYAmpVsZ2=ROOT.TH2D('hErrYAmpVsZ2','ferrY vs Z '+args.title,75,0,75,150,0,rmsY[1]) hErrZAmpVsZ2=ROOT.TH2D('hErrZAmpVsZ2','ferrZ vs Z '+args.title,75,0,75,150,0,rmsZ[1]) hErrLAmpVsZ2 =ROOT.TH2D('hErrLAmpVsZ2' ,'ferrL vs Z '+args.title ,75,0,75,100,0,rmsX[2]) hErrP1AmpVsZ2=ROOT.TH2D('hErrP1AmpVsZ2','ferrP1 vs Z '+args.title,75,0,75,100,0,rmsY[2]) hErrP2AmpVsZ2=ROOT.TH2D('hErrP2AmpVsZ2','ferrP2 vs Z '+args.title,75,0,75,100,0,rmsZ[2]) hErrLAmpVsZ =ROOT.TH2D('hErrLrAmpVsZ' ,'mcresT/ferrL vs Z '+args.title ,75,0,75,100,-rmsX[3],rmsX[3]) hErrP1AmpVsZ=ROOT.TH2D('hErrP1rAmpVsZ','mcresT/ferrP1 vs Z '+args.title,75,0,75,100,-rmsY[3],rmsY[3]) hErrP2AmpVsZ=ROOT.TH2D('hErrP2rAmpVsZ','mcresT/ferrP2 vs Z '+args.title,75,0,75,100,-rmsZ[3],rmsZ[3]) hErrA1AmpVsZ=ROOT.TH2D('hErrA1rAmpVsZ','mcresProj/ferrA1 vs Z '+args.title,75,0,75,150,-rmsX[4],rmsX[4]) hErrA2AmpVsZ=ROOT.TH2D('hErrA2rAmpVsZ','mcresProj/ferrA2 vs Z '+args.title,75,0,75,150,-rmsY[4],rmsY[4]) hErrA3AmpVsZ=ROOT.TH2D('hErrA3rAmpVsZ','mcresProj/ferrA3 vs Z '+args.title,75,0,75,150,-rmsZ[4],rmsZ[4]) hErrA1AmpVsZ2=ROOT.TH2D('hErrA1AmpVsZ2','ferrA1 vs Z '+args.title,75,0,75,100,0,rmsX[5]) hErrA2AmpVsZ2=ROOT.TH2D('hErrA2AmpVsZ2','ferrA2 vs Z '+args.title,75,0,75,100,0,rmsY[5]) hErrA3AmpVsZ2=ROOT.TH2D('hErrA3AmpVsZ2','ferrA3 vs Z '+args.title,75,0,75,100,0,rmsZ[5]) hErrSA1AmpVsZ2=ROOT.TH2D('hErrSA1AmpVsZ2','(fit,geom cut) ferrA1 vs Z '+args.title,75,0,75,100,0,rmsX[5]) hErrSA2AmpVsZ2=ROOT.TH2D('hErrSA2AmpVsZ2','(fit,geom cut) ferrA2 vs Z '+args.title,75,0,75,100,0,rmsY[5]) hErrSA3AmpVsZ2=ROOT.TH2D('hErrSA3AmpVsZ2','(fit,geom cut) ferrA3 vs Z '+args.title,75,0,75,100,0,rmsZ[5]) hResA1AmpVsZ2=ROOT.TH2D('hResA1AmpVsZ2','Res A1 vs Z '+args.title,75,0,75,100,-rmsX[14],rmsX[14]) hResA2AmpVsZ2=ROOT.TH2D('hResA2AmpVsZ2','Res A2 vs Z '+args.title,75,0,75,100,-rmsY[14],rmsY[14]) hResA3AmpVsZ2=ROOT.TH2D('hResA3AmpVsZ2','Res A3 vs Z '+args.title,75,0,75,100,-rmsZ[14],rmsZ[14]) hErrA1AmpVsZ3=ROOT.TH2D('hErrA1AmpVsZ3','mcresProj/ferrA1 (fit,geom cut) vs Z '+args.title,75,0,75,150,-rmsX[6],rmsX[6]) hErrA2AmpVsZ3=ROOT.TH2D('hErrA2AmpVsZ3','mcresProj/ferrA2 (fit,geom cut) vs Z '+args.title,75,0,75,150,-rmsY[6],rmsY[6]) hErrA3AmpVsZ3=ROOT.TH2D('hErrA3AmpVsZ3','mcresProj/ferrA3 (fit,geom cut) vs Z '+args.title,75,0,75,150,-rmsZ[6],rmsZ[6]) hErrNumVsZ=ROOT.TH2D('hErrNumVsZ','COV mcres/ferrProj '+args.title,75,0,75,125,0,rmsX[7]) hErrLamVsZ=ROOT.TH2D('hErrLamVsZ','COV ferrProj '+args.title,75,0,75,125,0,rmsY[7]) hResMagVsZ=ROOT.TH2D('hResMagVsZ','MC Res Mag '+args.title,75,0,75,125,0,rmsZ[7]) hScaleA1VsZ=ROOT.TH2D('hScaleA1VsZ','Diff A1 scale '+args.title,75,0,75,1000,0,5*scalefunc(73,0,ROOT.TVector3(1,0,0))) hScaleA2VsZ=ROOT.TH2D('hScaleA2VsZ','Diff A2 scale '+args.title,75,0,75,1000,0,5*scalefunc(73,0,ROOT.TVector3(0,1,0))) hScaleA3VsZ=ROOT.TH2D('hScaleA3VsZ','Diff A3 scale '+args.title,75,0,75,1000,0,5*scalefunc(73,0,ROOT.TVector3(0,0,1))) hErrUvsZ=ROOT.TH2D('hErrUvsZ','Err U '+args.title,75,0,75,100,0,rmsX[8]) hErrVvsZ=ROOT.TH2D('hErrVvsZ','Err V '+args.title,75,0,75,100,0,rmsY[8]) hErrUVvsZ=ROOT.TH2D('hErrUVvsZ','Err UV '+args.title,75,0,75,100,0,rmsZ[8]) hResUvsZ=ROOT.TH2D('hResUvsZ','Res U '+args.title,75,0,75,100,-rmsX[9],rmsX[9]) hResVvsZ=ROOT.TH2D('hResVvsZ','Res V '+args.title,75,0,75,100,-rmsY[9],rmsY[9]) hResUVvsZ=ROOT.TH2D('hResUVvsZ','Res UV '+args.title,75,0,75,100,-rmsZ[9],rmsZ[9]) hPullUvsZ=ROOT.TH2D('hPullUvsZ','Pull U '+args.title,75,0,75,100,-rmsX[10],rmsX[10]) hPullVvsZ=ROOT.TH2D('hPullVvsZ','Pull V '+args.title,75,0,75,100,-rmsY[10],rmsY[10]) hPullUVvsZ=ROOT.TH2D('hPullUVvsZ','Pull UV '+args.title,75,0,75,100,-rmsZ[10],rmsZ[10]) hErrXCov=ROOT.TH2D('hErrXC','Err X Cov cut '+args.title,75,0,75,100,0,rmsX[11]) hErrYCov=ROOT.TH2D('hErrYC','Err Y Cov cut '+args.title,75,0,75,100,0,rmsY[11]) hErrZCov=ROOT.TH2D('hErrZC','Err Z Cov cut '+args.title,75,0,75,100,0,rmsZ[11]) hPullXCov=ROOT.TH2D('hPullXCov','Pull X Cov cut '+args.title,75,0,75,100,-rmsX[12],rmsX[12]) hPullYCov=ROOT.TH2D('hPullYCov','Pull Y Cov cut '+args.title,75,0,75,100,-rmsY[12],rmsY[12]) hPullZCov=ROOT.TH2D('hPullZCov','Pull Z Cov cut '+args.title,75,0,75,100,-rmsZ[12],rmsZ[12]) hResX=ROOT.TH2D('hResX','Res X '+args.title,75,0,75,100,-rmsX[13],rmsX[13]) hResY=ROOT.TH2D('hResY','Res Y '+args.title,75,0,75,100,-rmsY[13],rmsY[13]) hResZ=ROOT.TH2D('hResZ','Res Z '+args.title,75,0,75,100,-rmsZ[13],rmsZ[13]) for i in range(3): htheta.append(ROOT.TH2D('htheta{0}'.format(i),'theta axis {0} {1}'.format(i,args.title),75,0,75,180,0,180)) hphi.append(ROOT.TH2D('hphi{0}'.format(i),'phi axis {0} {1}'.format(i,args.title),75,0,75,360,-180,180)) hdangle.append(ROOT.TH2D('hdangle{0}'.format(i),'angle to mom {0} {1}'.format(i,args.title),75,0,75,360,-180,180)) hdanglex.append(ROOT.TH2D('hdanglex{0}'.format(i),'angle to X {0} {1}'.format(i,args.title),75,0,75,360,-180,180)) hdangley.append(ROOT.TH2D('hdangley{0}'.format(i),'angle to Y {0} {1}'.format(i,args.title),75,0,75,360,-180,180)) hdanglez.append(ROOT.TH2D('hdanglez{0}'.format(i),'angle to Z {0} {1}'.format(i,args.title),75,0,75,360,-180,180)) ev=0 for e in tree: ev+=1 if args.events!=-1 and ev>args.events: break if ev%100==0: print 'Events:',ev # if e.CutCosmics.GetEntries()>1: # continue for tr in e.CutCosmics: phi=tr.GetPhi() theta=tr.GetTheta() amp=tr.GetAmps() zpos=tr.GetHitPositionsZ() ids=tr.GetHitIDs() mcresX=tr.GetResMcX() mcresY=tr.GetResMcY() mcresZ=tr.GetResMcZ() mcresU=tr.GetResMCU() mcresV=tr.GetResMCV() errU=tr.GetErrU() errV=tr.GetErrV() errMCRes=tr.GetSigResMC() trresX=tr.GetResX() trresY=tr.GetResY() trresZ=tr.GetResZ() trresU=tr.GetResU() trresV=tr.GetResU() trmcres=tr.GetResMCTrack() trmcresUV=tr.GetResMCTrackUV() errX=tr.GetSigX() errY=tr.GetSigY() errZ=tr.GetSigZ() errT=tr.GetSigT() errC=tr.GetSigC() axes=[tr.GetAxis1(),tr.GetAxis2(),tr.GetAxis3()] moms=tr.GetProjectionPointsMom() fitgood=tr.GetClusterFitGood() isGeomErr=tr.GetIsGeomErr() sizes=tr.GetClSizes() sizes2D=tr.Get2DClSizes() for cl in range(len(amp)): iax1=-1 iax2=-1 iax3=-1 minMOMangle=99 minYangle=99 minZangle=99 minXangle=99 min3angle=99 Xaxis=ROOT.TVector3(1,0,0) Yaxis=ROOT.TVector3(0,1,0) Zaxis=ROOT.TVector3(0,0,1) compax=Zaxis if args.sort: axc=-1 for ax in axes: axc+=1 if abs(ax[cl].Angle(moms[cl]))45 and abs(mtheta)<135: compax=Zaxis else: compax=Yaxis #angle if abs(ax[cl].Angle(compax)) < min3angle and axc!=iax1: min3angle=abs(ax[cl].Angle(compax)) iax3=axc #axc=-1 #for ax in axes: # axc+=1 # if ax[cl].Angle(Xaxis)90: dangleAx3Z=180-dangleAx3Z if dangleAx3Zargs.ax3Zax[1]: continue if sizes2D[cl]args.clsize2D[1]: continue if sizes[cl]args.clsize3D[1]: continue cov=tr.GetClusterCov(cl) scale,lam=calcCovScale(cov,res,amp[cl]) resUV=(resU+resV)/2. errUV=math.sqrt(errU[cl]**2+errV[cl]**2) if (fitgood[cl] and isGeomErr[cl].Mag2()==1) or args.stderr: if errU[cl]!=0 and resU!=0: hPullUvsZ.Fill(zpos[cl],resU/errU[cl]) if errV[cl]!=0 and resV!=0: hPullVvsZ.Fill(zpos[cl],resV/errV[cl]) if errUV!=0 and resUV!=0: hPullUVvsZ.Fill(zpos[cl],resUV/errUV) hErrUvsZ.Fill(zpos[cl],errU[cl]) hErrVvsZ.Fill(zpos[cl],errV[cl]) hErrUVvsZ.Fill(zpos[cl],errUV) hResUvsZ.Fill(zpos[cl],resU) hResVvsZ.Fill(zpos[cl],resV) hResUVvsZ.Fill(zpos[cl],resUV) hErrNumVsZ.Fill(zpos[cl],scale) hErrLamVsZ.Fill(zpos[cl],lam) hResMagVsZ.Fill(zpos[cl],res.Mag()) hScaleA1VsZ.Fill(zpos[cl],scalefunc(zpos[cl],0,axes[iax1][cl])) hScaleA2VsZ.Fill(zpos[cl],scalefunc(zpos[cl],1,axes[iax2][cl])) hScaleA3VsZ.Fill(zpos[cl],scalefunc(zpos[cl],2,axes[iax3][cl])) hErrXCov.Fill(zpos[cl],errC[cl][0]) hErrYCov.Fill(zpos[cl],errC[cl][1]) hErrZCov.Fill(zpos[cl],errC[cl][2]) if errC[cl][0]>0: hPullXCov.Fill(zpos[cl],res.X()/errC[cl][0]) if errC[cl][1]>0: hPullYCov.Fill(zpos[cl],res.Y()/errC[cl][1]) if errC[cl][2]>0: hPullZCov.Fill(zpos[cl],res.Z()/errC[cl][2]) hResX.Fill(zpos[cl],res.X()) hResY.Fill(zpos[cl],res.Y()) hResZ.Fill(zpos[cl],res.Z()) if 1: resT=calcprojT(res,moms[cl].Phi()) if errX[cl]!=0: hErrXAmpVsZ.Fill(zpos[cl],resX/errX[cl]) hErrXAmpVsZ2.Fill(zpos[cl],errX[cl]) if errY[cl]>0: hErrYAmpVsZ.Fill(zpos[cl],resY/errY[cl]) hErrYAmpVsZ2.Fill(zpos[cl],errY[cl]) #else: #print 'errY<0:',amp[cl],errY[cl] if errZ[cl]>0: hErrZAmpVsZ.Fill(zpos[cl],resZ/errZ[cl]) hErrZAmpVsZ2.Fill(zpos[cl],errZ[cl]) #print 'axes:' #axes[0][cl].Print() #axes[1][cl].Print() #axes[2][cl].Print() #print '' #print 'proj1' hErrA1AmpVsZ.Fill(zpos[cl],calcproj(res,axes[iax1][cl],amp[cl])) #print 'proj2' hErrA2AmpVsZ.Fill(zpos[cl],calcproj(res,axes[iax2][cl],amp[cl])) #print 'proj3' hErrA3AmpVsZ.Fill(zpos[cl],calcproj(res,axes[iax3][cl],amp[cl])) #print 'err axis:',calcproj2(axes[0][cl],amp[cl]),calcproj2(axes[1][cl],amp[cl]),calcproj2(axes[2][cl],amp[cl]) #print 'err xyz:',errX[cl],errY[cl],errZ[cl] #print '' hErrA1AmpVsZ2.Fill(zpos[cl],calcproj2(axes[iax1][cl],amp[cl])) hErrA2AmpVsZ2.Fill(zpos[cl],calcproj2(axes[iax2][cl],amp[cl])) hErrA3AmpVsZ2.Fill(zpos[cl],calcproj2(axes[iax3][cl],amp[cl])) hResA1AmpVsZ2.Fill(zpos[cl],calcAxProj(res,axes[iax1][cl])) hResA2AmpVsZ2.Fill(zpos[cl],calcAxProj(res,axes[iax2][cl])) hResA3AmpVsZ2.Fill(zpos[cl],calcAxProj(res,axes[iax3][cl])) if (fitgood[cl] and isGeomErr[cl].Mag2()==1) or args.stderr : hErrA1AmpVsZ3.Fill(zpos[cl],calcproj3(res,axes[iax1][cl])) hErrA2AmpVsZ3.Fill(zpos[cl],calcproj3(res,axes[iax2][cl])) hErrA3AmpVsZ3.Fill(zpos[cl],calcproj3(res,axes[iax3][cl])) hErrSA1AmpVsZ2.Fill(zpos[cl],calcproj2(axes[iax1][cl],amp[cl])) hErrSA2AmpVsZ2.Fill(zpos[cl],calcproj2(axes[iax2][cl],amp[cl])) hErrSA3AmpVsZ2.Fill(zpos[cl],calcproj2(axes[iax3][cl],amp[cl])) if errT[cl].X()!=0: hErrLAmpVsZ.Fill(zpos[cl],resT.X()/errT[cl].X()) if errT[cl].Y()!=0: hErrP1AmpVsZ.Fill(zpos[cl],resT.Y()/errT[cl].Y()) if errT[cl].Z()!=0: hErrP2AmpVsZ.Fill(zpos[cl],resT.Z()/errT[cl].Z()) hErrLAmpVsZ2.Fill(zpos[cl],errT[cl].X()/amp[cl]) hErrP1AmpVsZ2.Fill(zpos[cl],errT[cl].Y()/amp[cl]) hErrP2AmpVsZ2.Fill(zpos[cl],errT[cl].Z()/amp[cl]) rfile2=ROOT.TFile(args.basepath+'/'+args.outfile,'recreate') hErrXAmpVsZ.Write() hErrYAmpVsZ.Write() hErrZAmpVsZ.Write() hErrXAmpVsZ2.Write() hErrYAmpVsZ2.Write() hErrZAmpVsZ2.Write() hErrLAmpVsZ2.Write() hErrP1AmpVsZ2.Write() hErrP2AmpVsZ2.Write() hErrA1AmpVsZ.Write() hErrA2AmpVsZ.Write() hErrA3AmpVsZ.Write() hErrA1AmpVsZ2.Write() hErrA2AmpVsZ2.Write() hErrA3AmpVsZ2.Write() hErrSA1AmpVsZ2.Write() hErrSA2AmpVsZ2.Write() hErrSA3AmpVsZ2.Write() hResA1AmpVsZ2.Write() hResA2AmpVsZ2.Write() hResA3AmpVsZ2.Write() hErrA1AmpVsZ3.Write() hErrA2AmpVsZ3.Write() hErrA3AmpVsZ3.Write() hErrLAmpVsZ.Write() hErrP1AmpVsZ.Write() hErrP2AmpVsZ.Write() hErrNumVsZ.Write() hErrLamVsZ.Write() hResMagVsZ.Write() hScaleA1VsZ.Write() hScaleA2VsZ.Write() hScaleA3VsZ.Write() hErrUvsZ.Write() hErrVvsZ.Write() hErrUVvsZ.Write() hResUvsZ.Write() hResVvsZ.Write() hResUVvsZ.Write() hPullUvsZ.Write() hPullVvsZ.Write() hPullUVvsZ.Write() hErrXCov.Write() hErrYCov.Write() hErrZCov.Write() hPullXCov.Write() hPullYCov.Write() hPullZCov.Write() hResX.Write() hResY.Write() hResZ.Write() for i in range(3): htheta[i].Write() hphi[i].Write() hdangle[i].Write() hdanglex[i].Write() hdangley[i].Write() hdanglez[i].Write() if args.autorms: hrmsNum.Write() for i in range(len(hrmsX)): hrmsX[i].Write() hrmsY[i].Write() hrmsZ[i].Write() print 'histos written to:',args.basepath+'/'+args.outfile rfile2.Close() else: rfile=ROOT.TFile(args.recofile,'r') hErrXAmpVsZ=rfile.Get('hErrXrAmpVsZ') hErrYAmpVsZ=rfile.Get('hErrYrAmpVsZ') hErrZAmpVsZ=rfile.Get('hErrZrAmpVsZ') hErrXAmpVsZ2=rfile.Get('hErrXAmpVsZ2') hErrYAmpVsZ2=rfile.Get('hErrYAmpVsZ2') hErrZAmpVsZ2=rfile.Get('hErrZAmpVsZ2') hErrLAmpVsZ2=rfile.Get('hErrLAmpVsZ2') hErrP1AmpVsZ2=rfile.Get('hErrP1AmpVsZ2') hErrP2AmpVsZ2=rfile.Get('hErrP2AmpVsZ2') hErrLAmpVsZ=rfile.Get('hErrLrAmpVsZ') hErrP1AmpVsZ=rfile.Get('hErrP1rAmpVsZ') hErrP2AmpVsZ=rfile.Get('hErrP2rAmpVsZ') hErrA1AmpVsZ=rfile.Get('hErrA1rAmpVsZ') hErrA2AmpVsZ=rfile.Get('hErrA2rAmpVsZ') hErrA3AmpVsZ=rfile.Get('hErrA3rAmpVsZ') hErrA1AmpVsZ2=rfile.Get('hErrA1AmpVsZ2') hErrA2AmpVsZ2=rfile.Get('hErrA2AmpVsZ2') hErrA3AmpVsZ2=rfile.Get('hErrA3AmpVsZ2') hErrSA1AmpVsZ2=rfile.Get('hErrSA1AmpVsZ2') hErrSA2AmpVsZ2=rfile.Get('hErrSA2AmpVsZ2') hErrSA3AmpVsZ2=rfile.Get('hErrSA3AmpVsZ2') hErrA1AmpVsZ3=rfile.Get('hErrA1AmpVsZ3') hErrA2AmpVsZ3=rfile.Get('hErrA2AmpVsZ3') hErrA3AmpVsZ3=rfile.Get('hErrA3AmpVsZ3') hErrNumVsZ=rfile.Get('hErrNumVsZ') hErrLamVsZ=rfile.Get('hErrLamVsZ') hResMagVsZ=rfile.Get('hResMagVsZ') hScaleA1VsZ=rfile.Get('hScaleA1VsZ') hScaleA2VsZ=rfile.Get('hScaleA2VsZ') hScaleA3VsZ=rfile.Get('hScaleA3VsZ') hErrUvsZ=rfile.Get('hErrUvsZ') hErrVvsZ=rfile.Get('hErrVvsZ') hErrUVvsZ=rfile.Get('hErrUVvsZ') hResUvsZ=rfile.Get('hResUvsZ') hResVvsZ=rfile.Get('hResVvsZ') hResUVvsZ=rfile.Get('hResUVvsZ') hPullUvsZ=rfile.Get('hPullUvsZ') hPullVvsZ=rfile.Get('hPullVvsZ') hPullUVvsZ=rfile.Get('hPullUVvsZ') hResA1AmpVsZ2=rfile.Get('hResA1AmpVsZ2') hResA2AmpVsZ2=rfile.Get('hResA2AmpVsZ2') hResA3AmpVsZ2=rfile.Get('hResA3AmpVsZ2') hErrXCov=rfile.Get('hErrXC') hErrYCov=rfile.Get('hErrYC') hErrZCov=rfile.Get('hErrZC') hPullXCov=rfile.Get('hPullXCov') hPullYCov=rfile.Get('hPullYCov') hPullZCov=rfile.Get('hPullZCov') hResX=rfile.Get('hResX') hResY=rfile.Get('hResY') hResZ=rfile.Get('hResZ') for i in range(3): htheta.append(rfile.Get('htheta{0}'.format(i))) hphi.append(rfile.Get('hphi{0}'.format(i))) hdangle.append(rfile.Get('hdangle{0}'.format(i))) hdanglex.append(rfile.Get('hdanglex{0}'.format(i))) hdangley.append(rfile.Get('hdangley{0}'.format(i))) hdanglez.append(rfile.Get('hdanglez{0}'.format(i))) global allGraph global allProj global allData allProj={} allData={} allGraph=[] doXYZ=1 doAxis=1 doUV=1 doXYZCov=1 if doXYZ: funcstr='' funcstr+='1+[0]' funcstr+='+[1]*sqrt(x)*'+str(Dl) funcstr+='+[2]*exp(-1*x^[3]*'+str(Dl)+')' if args.phi!=-1: funcstr+='+'+str(args.phi)+'*[4]' func=ROOT.TF1('func',funcstr,0,73) func.SetLineWidth(1) func.SetTitle('X') print 'function for X:' print funcstr func.SetParameter(0,11) func.SetParameter(1,12) func.SetParameter(2,94) func.SetParameter(3,2) funcstr='' funcstr+='1+[0]' funcstr+='+[1]*sqrt(x)*'+str(Dl) funcstr+='+[2]*exp(-1*x^[3]*'+str(Dl)+')' if args.phi!=-1: funcstr+='+'+str(args.phi)+'*[4]' funcY=ROOT.TF1('funcY',funcstr,0,73) funcY.SetLineWidth(1) funcY.SetTitle('Y') print 'function for Y:' print funcstr funcY.SetParameter(0,11) funcY.SetParameter(1,12) funcY.SetParameter(2,94) funcY.SetParameter(3,2) funcstr='' if args.phi!=-1: funcstr+='(' funcstr+='[0]' funcstr+='+[1]*x' funcstr+='+[2]*'+str(dz)+'/sqrt(12)*exp(-1*x^[3]*'+str(Dt)+')' if args.phi!=-1: funcstr+=')*[4]' funcZ=ROOT.TF1('funcZ',funcstr,0,73) funcZ.SetLineWidth(1) funcZ.SetTitle('Z') print 'function for Z:' print funcstr funcZ.SetParameter(3,1) funcZ.SetParameter(4,1) f1=ROOT.TF1('f1',"[0]",0,73) f2=ROOT.TF1('f2','[0]*sqrt(x)*'+str(Dl),0,73) f3=ROOT.TF1('f3','[0]*'+str(dx)+"/sqrt(12)*exp(-1*x^[1]*"+str(Dl)+")",0,73) cerr.cd(1) projsX=[] gErrXAmpVsZ=getSlices(hErrXAmpVsZ,projsX,args.sliceopt) hErrXAmpVsZ.Draw('colz') cerr.cd(4) gErrXAmpVsZ['Mean'].Draw(graphopt) cerr.cd(7) hsig=gErrXAmpVsZ['RMS'] hsig.Draw(graphopt) hsig.Fit(func,'Q') print 'errX: chi^2=',func.GetChisquare() if func.GetNDF()!=0: chiNDF=func.GetChisquare()/func.GetNDF() else: chiNDF=-999 print 'errX: chi^2/NDF=',chiNDF c2.cd() f1.SetParameter(0,func.GetParameter(0)) f2.SetParameter(0,func.GetParameter(1)) f3.SetParameter(0,func.GetParameter(2)) f3.SetParameter(1,func.GetParameter(3)) f1.SetLineColor(1) f2.SetLineColor(2) f3.SetLineColor(3) f1.Draw() f2.Draw("same") f3.Draw("same") c2.Update() cerr.cd(2) projsY=[] gErrYAmpVsZ=getSlices(hErrYAmpVsZ,projsY,args.sliceopt) hErrYAmpVsZ.Draw('colz') cerr.cd(5) gErrYAmpVsZ['Mean'].Draw(graphopt) cerr.cd(8) hsigY=gErrYAmpVsZ['RMS'] hsigY.Draw(graphopt) hsigY.Fit(funcY,'Q') print 'errY: chi^2=',funcY.GetChisquare() if funcY.GetNDF()!=0: chiNDF=funcY.GetChisquare()/funcY.GetNDF() else: chiNDF=-999 print 'errY: chi^2/NDF=',chiNDF cerr.Update() cerr.cd(3) projsZ=[] gErrZAmpVsZ=getSlices(hErrZAmpVsZ,projsZ,args.sliceopt) hErrZAmpVsZ.Draw('colz') cerr.cd(6) gErrZAmpVsZ['Mean'].Draw(graphopt) cerr.cd(9) hsigZ=gErrZAmpVsZ['RMS'] hsigZ.Draw(graphopt) hsigZ.Fit(funcZ,'Q') print 'errZ: chi^2=',funcZ.GetChisquare() if funcZ.GetNDF()!=0: chiNDF=funcZ.GetChisquare()/funcZ.GetNDF() else: chiNDF=-999 print 'errZ: chi^2/NDF=',chiNDF cerr.Update() #error 2 *************** cerr2.cd(1) hErrXAmpVsZ2.Draw('colz') cerr2.cd(2) hErrYAmpVsZ2.Draw('colz') cerr2.cd(3) hErrZAmpVsZ2.Draw('colz') gErrXAmpVsZ2=[] fillGraph(hErrXAmpVsZ2,gErrXAmpVsZ2) cerr2.cd(4) gErrXAmpVsZ2[0]['Mean'].Draw("AP") cerr2.cd(7) gErrXAmpVsZ2[0]['RMS'].Draw("AP") gErrYAmpVsZ2=[] fillGraph(hErrYAmpVsZ2,gErrYAmpVsZ2) cerr2.cd(5) gErrYAmpVsZ2[0]['Mean'].Draw("AP") cerr2.cd(8) gErrYAmpVsZ2[0]['RMS'].Draw("AP") gErrZAmpVsZ2=[] fillGraph(hErrZAmpVsZ2,gErrZAmpVsZ2) cerr2.cd(6) gErrZAmpVsZ2[0]['Mean'].Draw("AP") cerr2.cd(9) gErrZAmpVsZ2[0]['RMS'].Draw("AP") cerr2.Update() # end error2 ******************************** ''' #error T2 ****************************** cerrT2.Divide(3,3) cerrT2.cd(1) hErrLAmpVsZ2.Draw('colz') cerrT2.cd(2) hErrP1AmpVsZ2.Draw('colz') cerrT2.cd(3) hErrP2AmpVsZ2.Draw('colz') cerrT2.Update() gErrLAmpVsZ2=[] fillGraph(hErrLAmpVsZ2,gErrLAmpVsZ2) cerrT2.cd(4) gErrLAmpVsZ2[0]['Mean'].Draw("AP") cerrT2.cd(7) gErrLAmpVsZ2[0]['RMS'].Draw("AP") gErrP1AmpVsZ2=[] fillGraph(hErrP1AmpVsZ2,gErrP1AmpVsZ2) cerrT2.cd(5) gErrP1AmpVsZ2[0]['Mean'].Draw("AP") cerrT2.cd(8) gErrP1AmpVsZ2[0]['RMS'].Draw("AP") gErrP2AmpVsZ2=[] fillGraph(hErrP2AmpVsZ2,gErrP2AmpVsZ2) cerrT2.cd(6) gErrP2AmpVsZ2[0]['Mean'].Draw("AP") cerrT2.cd(9) gErrP2AmpVsZ2[0]['RMS'].Draw("AP") cerrT2.Update() #end error T2 *************************** #error T ****************************** cerrT.Divide(3,3) cerrT.cd(1) projsL=[] gErrLAmpVsZ=getSlices(hErrLAmpVsZ,projsL,args.sliceopt) hErrLAmpVsZ.Draw('colz') cerrT.cd(4) gErrLAmpVsZ['Mean'].Draw(graphopt) funcstr='' funcstr+='[0]' funcstr+='+[1]*sqrt(x)'#*'+str(Dl) funcstr+='+[2]*exp(-x^[3])'#+str(Dl)+')' funcL=ROOT.TF1('funcL',funcstr,0,73) funcL.SetLineWidth(1) funcL.SetTitle('L') print 'function for L:' print funcstr hrmsL=gErrLAmpVsZ['RMS'] parsL=getStartPars(hrmsL) for i in range(len(parsL)): funcL.SetParameter(i,parsL[i]) cerrT.cd(7) hrmsL.Fit(funcL,'Q') hrmsL.Draw(graphopt) cerrT.cd(2) projsP1=[] gErrP1AmpVsZ=getSlices(hErrP1AmpVsZ,projsP1,args.sliceopt) hErrP1AmpVsZ.Draw('colz') cerrT.cd(5) gErrP1AmpVsZ['Mean'].Draw(graphopt) funcP1=ROOT.TF1('funcP1',funcstr,0,73) funcP1.SetLineWidth(1) funcP1.SetTitle('P1') print 'function used for P1:' print funcstr hrmsP1=gErrP1AmpVsZ['RMS'] parsP1=getStartPars(hrmsP1) for i in range(len(parsP1)): funcP1.SetParameter(i,parsP1[i]) cerrT.cd(8) hrmsP1.Fit(funcP1,'Q') hrmsP1.Draw(graphopt) cerrT.cd(3) projsP2=[] gErrP2AmpVsZ=getSlices(hErrP2AmpVsZ,projsP2,args.sliceopt) hErrP2AmpVsZ.Draw('colz') cerrT.cd(6) gErrP2AmpVsZ['Mean'].Draw(graphopt) funcstr='' funcstr+='[0]' funcstr+='+[1]*sqrt(x)*'+str(Dt) funcstr+='+[2]*exp(-x^[3]*'+str(Dt)+')' funcP2=ROOT.TF1('funcP2',funcstr,0,73) funcP2.SetLineWidth(1) funcP2.SetTitle('P2') print 'function used for P2:' print funcstr hrmsP2=gErrP2AmpVsZ['RMS'] parsP2=getStartPars(hrmsP2) for i in range(len(parsP2)): funcP2.SetParameter(i,parsP2[i]) cerrT.cd(9) hrmsP2.Fit(funcP2,'Q') hrmsP2.Draw(graphopt) cerrT.Update() #end error T *************************** ''' if doAxis: # error A ********************** cerrA.Divide(3,3) cerrA.cd(1) projsA1=[] gErrA1AmpVsZ=getSlices(hErrA1AmpVsZ,projsA1,args.sliceopt) hErrA1AmpVsZ.Draw('colz') cerrA.cd(4) gErrA1AmpVsZ['Mean'].Draw(graphopt) funcstr='' funcstr+='[0]' funcstr+='+[1]*sqrt(x)' funcstr+='+[2]*exp(-(x^[3]))'#'+str(Dl)+')' funcA1=ROOT.TF1('funcA1',funcstr,0,73) funcA1.SetLineWidth(1) funcA1.SetTitle('A1') print 'function for A1:' print funcstr hrmsA1=gErrA1AmpVsZ['RMS'] parsA1=getStartPars(hrmsA1) for i in range(len(parsA1)): funcA1.SetParameter(i,parsA1[i]) if args.sliceopt==-1: funcA1.SetParLimits(3,0,10) if args.sliceopt==1: funcA1.SetParameter(3,0.4) funcA1.SetParameter(2,0.5) funcA1.SetParLimits(3,0.,5.) funcA1.SetParLimits(2,0.,10.) cerrA.cd(7) hrmsA1.Fit(funcA1,'') hrmsA1.Draw(graphopt) cprojA1.Divide(int(math.ceil(math.sqrt(len(projsA1)))), int(math.ceil(math.sqrt(len(projsA1))))) for ip in range(len(projsA1)): cprojA1.cd(ip+1) projsA1[ip].Draw() cprojA1.Update() debug=0 cerrA.cd(2) projsA2=[] gErrA2AmpVsZ=getSlices(hErrA2AmpVsZ,projsA2,args.sliceopt) hErrA2AmpVsZ.Draw('colz') cerrA.cd(5) gErrA2AmpVsZ['Mean'].Draw(graphopt) funcA2=ROOT.TF1('funcA2',funcstr,0,73) funcA2.SetLineWidth(1) funcA2.SetTitle('A2') print 'function used for A2:' print funcstr hrmsA2=gErrA2AmpVsZ['RMS'] parsA2=getStartPars(hrmsA2) for i in range(len(parsA2)): funcA2.SetParameter(i,parsA2[i]) if args.sliceopt==-1: funcA2.SetParLimits(3,0,10) #funcA2.FixParameter(1,0) cerrA.cd(8) hrmsA2.Fit(funcA2,'') hrmsA2.Draw(graphopt) cprojA2.Divide(int(math.ceil(math.sqrt(len(projsA2)))), int(math.ceil(math.sqrt(len(projsA2))))) for ip in range(len(projsA2)): cprojA2.cd(ip+1) projsA2[ip].Draw() cprojA2.Update() cerrA.cd(3) projsA3=[] gErrA3AmpVsZ=getSlices(hErrA3AmpVsZ,projsA3,args.sliceopt) hErrA3AmpVsZ.Draw('colz') cerrA.cd(6) gErrA3AmpVsZ['Mean'].Draw(graphopt) funcstr='' funcstr+='[0]' funcstr+='+[1]*sqrt(x)'#*'+str(Dt) funcstr+='+[2]*exp(-(x^[3]))'#*'+str(Dt)+')' funcA3=ROOT.TF1('funcA3',funcstr,0,73) funcA3.SetLineWidth(1) funcA3.SetTitle('A3') print 'function used for A3:' print funcstr hrmsA3=gErrA3AmpVsZ['RMS'] parsA3=getStartPars(hrmsA3) for i in range(len(parsA3)): funcA3.SetParameter(i,parsA3[i]) cerrA.cd(9) hrmsA3.Fit(funcA3,'') hrmsA3.Draw(graphopt) cprojA3.Divide(int(math.ceil(math.sqrt(len(projsA3)))), int(math.ceil(math.sqrt(len(projsA3))))) for ip in range(len(projsA3)): cprojA3.cd(ip+1) projsA3[ip].Draw() cprojA3.Update() cerrA.Update() debug=0 #end Error A ********************** #error A2 ********************** cerrA2.Divide(3,3) cerrA2.cd(1) hErrA1AmpVsZ2.Draw('colz') cerrA2.cd(2) hErrA2AmpVsZ2.Draw('colz') cerrA2.cd(3) hErrA3AmpVsZ2.Draw('colz') cerrA2.cd(4) hErrSA1AmpVsZ2.Draw('colz') cerrA2.cd(5) hErrSA2AmpVsZ2.Draw('colz') cerrA2.cd(6) hErrSA3AmpVsZ2.Draw('colz') cerrA2.cd(7) hResA1AmpVsZ2.Draw('colz') cerrA2.cd(8) hResA2AmpVsZ2.Draw('colz') cerrA2.cd(9) hResA3AmpVsZ2.Draw('colz') cerrA2.Update() #end error A2 ********************** #error A3 ************************** cerrA3.Divide(3,3) cerrA3.cd(1) projsA1_3=[] gErrA1AmpVsZ3=getSlices(hErrA1AmpVsZ3,projsA1_3,args.sliceopt) hErrA1AmpVsZ3.Draw('colz') cerrA3.cd(4) gErrA1AmpVsZ3['Mean'].Draw(graphopt) cerrA3.cd(7) funcA1_3=ROOT.TF1('funcA1_3',funcstr,0,73) funcA1_3.SetLineWidth(1) funcA1_3.SetTitle('A1_3') print 'function for proj3 A1:',funcstr hrmsA1_3=gErrA1AmpVsZ3['RMS'] hrmsA1_3.Fit(funcA1_3,'') hrmsA1_3.Draw(graphopt) cprojA1_3.Divide(int(math.ceil(math.sqrt(len(projsA1_3)))), int(math.ceil(math.sqrt(len(projsA1_3))))) for ip in range(len(projsA1_3)): cprojA1_3.cd(ip+1) projsA1_3[ip].Draw() cprojA1_3.Update() cerrA3.cd(2) projsA2_3=[] debug=0 gErrA2AmpVsZ3=getSlices(hErrA2AmpVsZ3,projsA2_3,args.sliceopt) hErrA2AmpVsZ3.Draw('colz') cerrA3.cd(5) gErrA2AmpVsZ3['Mean'].Draw(graphopt) cerrA3.cd(8) funcA2_3=ROOT.TF1('funcA2_3',funcstr,0,73) funcA2_3.SetLineWidth(1) funcA2_3.SetTitle('A2_3') hrmsA2_3=gErrA2AmpVsZ3['RMS'] hrmsA2_3.Fit(funcA2_3,'') hrmsA2_3.Draw(graphopt) cprojA2_3.Divide(int(math.ceil(math.sqrt(len(projsA2_3)))), int(math.ceil(math.sqrt(len(projsA2_3))))) for ip in range(len(projsA2_3)): cprojA2_3.cd(ip+1) projsA2_3[ip].Draw() cprojA2_3.Update() projsA3_3=[] gErrA3AmpVsZ3=getSlices(hErrA3AmpVsZ3,projsA3_3,args.sliceopt) cerrA3.cd(3) hErrA3AmpVsZ3.Draw('colz') cerrA3.cd(6) gErrA3AmpVsZ3['Mean'].Draw(graphopt) cerrA3.cd(9) funcA3_3=ROOT.TF1('funcA3_3',funcstr,0,73) funcA3_3.SetLineWidth(1) funcA3_3.SetTitle('A3_3') hrmsA3_3=gErrA3AmpVsZ3['RMS'] hrmsA3_3.Fit(funcA3_3,'') hrmsA3_3.Draw(graphopt) cprojA3_3.Divide(int(math.ceil(math.sqrt(len(projsA3_3)))), int(math.ceil(math.sqrt(len(projsA3_3))))) for ip in range(len(projsA3_3)): cprojA3_3.cd(ip+1) projsA3_3[ip].Draw() cprojA3_3.Update() cerrA3.Update() debug=0 #end error A3 ********************** #scaling factor ******************** print'now doing scaling factor' cerrNum.Divide(3,3) projsNum=[] gErrNumVsZ=getSlices(hErrNumVsZ,projsNum,args.sliceopt,1) cerrNum.cd(1) hErrLamVsZ.Draw('colz') cerrNum.cd(2) hResMagVsZ.Draw('colz') cerrNum.cd(3) hErrNumVsZ.Draw('colz') cerrNum.cd(6) gErrNumVsZ['Mean'].Draw(graphopt) cerrNum.cd(9) gErrNumVsZ['RMS'].Draw(graphopt) cerrNum.Update() cprojNum.Divide(int(math.ceil(math.sqrt(len(projsNum)))), int(math.ceil(math.sqrt(len(projsNum))))) for ip in range(len(projsNum)): cprojNum.cd(ip+1) projsNum[ip].Draw() cprojNum.Update() #end scaling factor **************** #**diff based scaling print 'now doing diff based scaling' cScale.Divide(3,3) cScale.cd(1) hScaleA1VsZ.Draw("colz") cScale.cd(2) hScaleA2VsZ.Draw("colz") cScale.cd(3) hScaleA3VsZ.Draw("colz") cScale.Update() #end diff based scaling if doUV: #Error and Residual UV print 'now doing UV error' cErrUV.Divide(3,3) cErrUV.cd(1) hErrUvsZ.Draw("colz") cErrUV.cd(2) hErrVvsZ.Draw("colz") cErrUV.cd(3) hErrUVvsZ.Draw("colz") cErrUV.cd(4) hResUvsZ.Draw("colz") cErrUV.cd(5) hResVvsZ.Draw("colz") cErrUV.cd(6) hResUVvsZ.Draw("colz") cErrUV.cd(7) cErrUV.cd(8) cErrUV.cd(9) cErrUV.Update() #End error and residual UV #Pull UV projsPullU=[] gPullUvsZ=getSlices(hPullUvsZ,projsPullU,args.sliceopt) projsPullV=[] gPullVvsZ=getSlices(hPullVvsZ,projsPullV,args.sliceopt) projsPullUV=[] gPullUVvsZ=getSlices(hPullUVvsZ,projsPullUV,args.sliceopt) cPullUV.Divide(3,3) cPullUV.cd(1) hPullUvsZ.Draw("colz") cPullUV.cd(2) hPullVvsZ.Draw("colz") cPullUV.cd(3) hPullUVvsZ.Draw("colz") cPullUV.cd(4) gPullUvsZ['Mean'].Draw(graphopt) cPullUV.cd(5) gPullVvsZ['Mean'].Draw(graphopt) cPullUV.cd(6) gPullUVvsZ['Mean'].Draw(graphopt) cPullUV.cd(7) gPullUvsZ['RMS'].Draw(graphopt) cPullUV.cd(8) gPullVvsZ['RMS'].Draw(graphopt) cPullUV.cd(9) gPullUVvsZ['RMS'].Draw(graphopt) cPullUV.Update() cprojPullU.Divide(int(math.ceil(math.sqrt(len(projsPullUV)))), int(math.ceil(math.sqrt(len(projsPullUV))))) for ip in range(len(projsPullU)): cprojPullU.cd(ip+1) projsPullU[ip].Draw() cprojPullU.Update() cprojPullV.Divide(int(math.ceil(math.sqrt(len(projsPullUV)))), int(math.ceil(math.sqrt(len(projsPullUV))))) for ip in range(len(projsPullV)): cprojPullV.cd(ip+1) projsPullV[ip].Draw() cprojPullV.Update() #end Pull UV if doXYZCov: #err xyzCov cerrC.Divide(3,3) cerrC.cd(1) hErrXCov.Draw('colz') cerrC.cd(2) hErrYCov.Draw('colz') cerrC.cd(3) hErrZCov.Draw('colz') cerrC.cd(4) hResX.Draw('colz') cerrC.cd(5) hResY.Draw('colz') cerrC.cd(6) hResZ.Draw('colz') cerrC.cd(7) hErrXAmpVsZ2.Draw('colz') cerrC.cd(8) hErrYAmpVsZ2.Draw('colz') cerrC.cd(9) hErrZAmpVsZ2.Draw('colz') #end errxyzCov #pull xyzCov projsXC=[] projsYC=[] projsZC=[] gPullXC=getSlices(hPullXCov,projsXC,args.sliceopt) gPullYC=getSlices(hPullYCov,projsYC,args.sliceopt) gPullZC=getSlices(hPullZCov,projsZC,args.sliceopt) DrawPulls(cPullC,[hPullXCov,hPullYCov,hPullZCov],[gPullXC,gPullYC,gPullZC],graphopt) cPullC.Update() #end pullxyzCov #********************residuals cres.Divide(3,3) cres.cd(1) hResA1AmpVsZ2.Draw('colz') cres.cd(2) hResA2AmpVsZ2.Draw('colz') cres.cd(3) hResA3AmpVsZ2.Draw('colz') resproj1=[] gresmean1=getSlices(hResA1AmpVsZ2,resproj1,args.sliceopt) resproj2=[] gresmean2=getSlices(hResA2AmpVsZ2,resproj2,args.sliceopt) resproj3=[] gresmean3=getSlices(hResA3AmpVsZ2,resproj3,args.sliceopt) cres.cd(4) gresmean1['Mean'].Draw(graphopt) cres.cd(5) gresmean2['Mean'].Draw(graphopt) cres.cd(6) gresmean3['Mean'].Draw(graphopt) cres.cd(7) gresmean1['RMS'].Draw(graphopt) cres.cd(8) gresmean2['RMS'].Draw(graphopt) cres.cd(9) gresmean3['RMS'].Draw(graphopt) #******************end residuals #********************errors cerr3.Divide(3,3) cerr3.cd(1) hErrSA1AmpVsZ2.Draw('colz') cerr3.cd(2) hErrSA2AmpVsZ2.Draw('colz') cerr3.cd(3) hErrSA3AmpVsZ2.Draw('colz') errproj1=[] gresmean1=getSlices(hErrSA1AmpVsZ2,errproj1,1) errproj2=[] gerrmean2=getSlices(hErrSA2AmpVsZ2,errproj2,1) errproj3=[] gerrmean3=getSlices(hErrSA3AmpVsZ2,errproj3,1) cerr3.cd(4) gresmean1['Mean'].Draw(graphopt) cerr3.cd(5) gerrmean2['Mean'].Draw(graphopt) cerr3.cd(6) gerrmean3['Mean'].Draw(graphopt) cerr3.cd(7) gresmean1['RMS'].Draw(graphopt) cerr3.cd(8) gerrmean2['RMS'].Draw(graphopt) cerr3.cd(9) gerrmean3['RMS'].Draw(graphopt) #******************end errors cangle.Divide(3,3) cangle2.Divide(3,3) for i in range(3): cangle.cd( i+1 ) htheta[i].Draw('colz') cangle.cd( i+4 ) hphi[i].Draw('colz') cangle.cd( i+7 ) hdangle[i].Draw('colz') cangle2.cd( i+1 ) if hdanglex[i]!=None: hdanglex[i].Draw('colz') cangle2.cd( i+4 ) if hdangley[i]!=None: hdangley[i].Draw('colz') cangle2.cd( i+7 ) if hdanglez[i]!=None: hdanglez[i].Draw('colz') rfile2=ROOT.TFile(args.basepath+'/'+args.outfile,'update') cerrT2.Write() cerrT.Write() cerr.Write() cerr2.Write() cerrA.Write() cerrA2.Write() cerrA3.Write() cerrNum.Write() cScale.Write() cErrUV.Write() cPullUV.Write() cprojA1.Write() cprojA2.Write() cprojA3.Write() cprojA1_3.Write() cprojA2_3.Write() cprojA3_3.Write() cprojNum.Write() cprojPullU.Write() cprojPullV.Write() cerrC.Write() cPullC.Write() cangle.Write() cangle2.Write() for arr in allProj: rfile2.mkdir("{0}".format(arr)) rfile2.cd("{0}".format(arr)) print 'entering directory',arr for h in allProj[arr]: #print 'writing Proj histo ',h.GetName() h.Write() rfile2.cd() for arr in allGraph: for g in arr: #print 'writing graph',g arr[g].Write() carray=[cerrT2,cerrT,cerr,cerr2,cerrA,cerrA2,cerrA3,cres,cerr3,cerrNum,cScale,cErrUV,cPullUV,cerrC,cPullC,cangle,cangle2 #,cprojA1,cprojA2,cprojA3,cprojA1_3,cprojA2_3,cprojA3_3,cprojNum,cprojPullU,cprojPullV ] #cerry.Write() #c1.Write() rfile2.Close() print 'histos and canvas written to:',args.basepath+'/'+args.outfile if args.fitfile!="None": outfile=open(args.basepath+'/'+args.fitfile,'a') funcs=[func,funcY,funcZ,funcA1,funcA2,funcA3,funcA1_3,funcA2_3,funcA3_3] for f in funcs: line='' line+=f.GetTitle()+';' line+=str(f.GetParameter(0))+';' line+=str(f.GetParameter(1))+';' line+=str(f.GetParameter(2))+';' line+=str(f.GetParameter(3))+';' if args.phi!=-1: line+=str(f.GetParameter(4))+';' line+='\n' outfile.write(line) print 'fitinfo written to:',args.basepath+'/'+args.fitfile outfile.close ''' if args.plot and args.pfile=='None': path=args.recofile[0:args.recofile.rfind('/')+1] if not os.path.exists(path): os.system('mkdir '+path) for c in carray: if args.genfit: c.SaveAs(path+c.GetTitle()+'_tr.eps') else: c.SaveAs(path+c.GetTitle()+'_mc.eps') ''' if args.pfile !='None': carray[0].SaveAs(args.pfile+'[') for c in carray: c.SaveAs(args.pfile) for arr in allProj: for h in allProj[arr]: cproj.cd() h.Draw() if allData.get(h.GetTitle(),None)!=None: drawGaussians(allData[h.GetTitle()]) cproj.Update() cproj.SaveAs(args.pfile) carray[-1].SaveAs(args.pfile+']') if not args.batch: u='n' while u!='q': u=raw_input("done?")