import sys,os,time,math, array sys.path.append('/home/mberger/fopiroot/fopiroot_dev/macro/tpc/FOPI/mberger') pandapath=os.environ.get('PANDAPATH') sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1') import glob,math,argparse,ROOT from functions import set_palette def drawHist(c,h,fopt,u,v): latex=ROOT.TLatex() text1=[] hproj1=h.ProjectionY() hproj2=h.ProjectionX() hproj1.SetTitle('Projection 1') hproj2.SetTitle('Projection 2') c.cd(1) h.SetStats(0) h.Draw('colz') c.cd(2) fitfunc.SetParameter(0,hproj1.GetMaximum()) fitfunc.SetParameter(1,hproj1.GetMean()) fitfunc.SetParameter(2,hproj1.GetRMS()) hproj1.Fit(fitfunc,fopt) text1.append('Projection 1 Fit:') text1.append('Mean: ' +str(fitfunc.GetParameter(1))+' #pm '+str(fitfunc.GetParError(1)) ) hproj1.Draw() c.cd(3) hproj2.Draw() fitfunc.SetParameter(0,hproj2.GetMaximum()) fitfunc.SetParameter(1,hproj2.GetMean()) fitfunc.SetParameter(2,hproj2.GetRMS()) hproj2.Fit(fitfunc,fopt) text1.append('Projection 2 Fit:') text1.append('Mean: ' +str(fitfunc.GetParameter(1))+' #pm '+str(fitfunc.GetParError(1))) c.cd(4) text1.append('U: theta='+str(u.Theta()*ROOT.TMath.RadToDeg()) +', Phi='+str(u.Phi()*ROOT.TMath.RadToDeg())) text1.append('V: theta='+str(v.Theta()*ROOT.TMath.RadToDeg()) +', Phi='+str(v.Phi()*ROOT.TMath.RadToDeg())) x=0 y=.9 for t in text1: latex.DrawLatex(x,y,t) y-=0.1 c.Update() def drawText(text): latex=ROOT.TLatex() x=0 y=.9 for t in text: latex.DrawLatex(x,y,t) y-=0.1 def multIf(var,m1,m2): if var>0: return m1 else: return m2 def drawGraph(c,g1,g2,g3,fopt,u,v): latex=ROOT.TLatex() text1=[] c.cd(1) setGraph(g1) g1.Draw('pcol') c.cd(2) setGraph(g2) g2.Draw('AP') c.cd(3) setGraph(g3) g3.Draw('AP') def setGraph(g,msize=1,mstyle=20,mcol=1): g.SetMarkerSize(msize) g.SetMarkerStyle(mstyle) g.SetMarkerColor(mcol) def drawCov(cov,amp,tup,c): c.cd() tcov=ROOT.TMatrixDSym(2,cov.GetMatrixArray()) tcov*=1./amp eigenproblem=ROOT.TMatrixDEigen(ROOT.TMatrixD(tcov)) eigenwert=eigenproblem.GetEigenValues() eigenvector=eigenproblem.GetEigenVectors() teigenvector=ROOT.TVector2(eigenvector[0][0],eigenvector[1][0]) teigenvector2=ROOT.TVector2(eigenvector[0][1],eigenvector[1][1]) eigenvector_i=ROOT.TMatrixD(ROOT.TMatrixD.kInverted,eigenvector) pcov=ROOT.TEllipse(clPosUV.X(),clPosUV.Y(),math.sqrt(abs(eigenwert[0][0])),math.sqrt(abs(eigenwert[1][1])),0,360,ROOT.TMath.RadToDeg()*teigenvector.Phi()) pcov.SetFillStyle(4000) h=ROOT.TH2D('h','h',100,tup.GetMinimum('u')-1,tup.GetMaximum('u')+1,100,tup.GetMinimum('v')-1,tup.GetMaximum('v')+1) for e in tup: h.Fill(tup.u,tup.v,tup.amp) h.Draw('colz') ppos.Draw('same') pcov.Draw('same') c.Update() u=raw_input('next!') #def vectorOnPlane(): def fillAxisProjHist(dtup,ev,cog,hc,ht,hg): ew1=ROOT.TVector2(eigenvector[0][0],eigenvector[1][0]) ew2=ROOT.TVector2(eigenvector[0][1],eigenvector[1][1]) ew1.SetMagPhi(1,ew1.Phi()) ew2.SetMagPhi(1,ew2.Phi()) retval=[] for dig in digiTuple: dpos=ROOT.TVector2(digiTuple.u-cog.X(),digiTuple.v-cog.Y()) dpose=ROOT.TVector2(dpos*ew1,dpos*ew2) retval.append(dpose) hc[0].Fill(dpose.X(),digiTuple.amp) hc[1].Fill(dpose.Y(),digiTuple.amp) hc[2].Fill(dpose.X()) hc[3].Fill(dpose.Y()) ht[0].Fill(dpose.X(),digiTuple.amp) ht[1].Fill(dpose.Y(),digiTuple.amp) ht[2].Fill(dpose.X()) ht[3].Fill(dpose.Y()) hg[0].Fill(dpose.X(),digiTuple.amp) hg[1].Fill(dpose.Y(),digiTuple.amp) hg[2].Fill(dpose.X()) hg[3].Fill(dpose.Y()) return retval def makeAxisProjHist(pos,hname,htit): posU=[] posV=[] for p in pos: posU.append(p.X()) posV.append(p.Y()) posU.sort() posV.sort() meanDU=0 meanDV=0 for ip in range(len(posU)-1): meanDU+=abs(posU[ip]-posU[ip+1]) meanDV+=abs(posV[ip]-posV[ip+1]) meanDU/=len(posU)-1 meanDV/=len(posV)-1 def gauss2d(x,par): #par=[sigX,sigY,rho,meanX,meanY] factor=2*ROOT.TMath.Pi()*par[0]*par[1] factor*=ROOT.TMath.Sqrt(1-par[2]*par[2]) if factor!=0: factor=1./factor val1=(x[0]-par[3])*(x[0]-par[3]) if par[0]!=0: val1/=par[0]*par[0] val2=(x[1]-par[4])*(x[1]-par[4]) if par[1]!=0: val2/=par[1]*par[1] val3=2*par[2]*(x[0]-par[3])*(x[1]-par[4]) if par[0]!=0 and par[1]!=0: val3/=par[0]*par[1] exponent=-1./(2*(1-par[2]*par[2])) exponent*=(val1+val2-val3) retval=factor*ROOT.TMath.Exp(exponent) return retval parser=argparse.ArgumentParser(description="plots digisposition projected on detector plane") parser.add_argument('recofile',help='the reco file',type=str,default='./') parser.add_argument('--outfile',help='the root file to write into',type=str,default='./fitcompare.root') parser.add_argument('--pfile',help='the pdf to write into',type=str,default='None') parser.add_argument('--pfile2',help='the pdf to write into',type=str,default='None') parser.add_argument('--events',help='number of events to run',type=int,default=-1) parser.add_argument('--watchthecov',help='draw cov at each step',action='store_true') parser.add_argument('--goto',help='go to event',type=int,default=-1) parser.add_argument('--clusterwise',help='stop at each cluster',action='store_true') parser.add_argument('--batch',help='batch mode',action='store_true') parser.add_argument('--incrRange',help='increase the range of the fit',action='store_true') parser.add_argument('--no2D',help='skip the 2D gauss fit',action='store_true') parser.add_argument('--zcut',help='only cluter inside z-window are taken into account',type=float,default=[-1,80],nargs=2) parser.add_argument args=parser.parse_args() ROOT.gROOT.ProcessLine(".x rootlogon.C") ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') ROOT.gROOT.LoadMacro("stlPYROOT.h+") ROOT.gROOT.LoadMacro("2dgauss.h+") ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.ERROR) ROOT.RooMsgService.instance().setSilentMode(True) set_palette() if args.recofile=='./': print 'no recofile given' exit() if args.batch: ROOT.gROOT.SetBatch(True) RFile=ROOT.TFile(args.recofile,'read') rfile=ROOT.TFile(args.outfile,'recreate') tree=RFile.Get('cbmsim') if tree==None: print 'no tree found' RFile.Close() exit() canv=ROOT.TCanvas('canv','digis on plane',1000,1000) canv.Divide(2,2) if args.watchthecov: ctemp=ROOT.TCanvas('ctemp','ctemp') fitfunc=ROOT.TF1('gaus','gaus',-2,2) fitfunc.SetLineWidth(1) fitfunc.SetLineColor(2) if args.pfile!='None': canv.SaveAs(args.pfile+'[') ev=-1 fittime=ROOT.TH1D('fittime','fittime multi gauss',1000,0,2) fittime2=ROOT.TH1D('fittime2','fittime 2dgauss',1000,0,2) diffhistos={} diffhistos['fittime']=ROOT.TH1D('dfittime','Fittime Diff',1000,-5,5) diffhistos['sigmaU']=ROOT.TH1D('dsigmaU','SigmaU Diff',1000,-0.5,0.5) diffhistos['sigmaV']=ROOT.TH1D('dsigmaV','SigmaV Diff',1000,-0.5,0.5) diffhistos['rho']=ROOT.TH1D('drho','Rho Diff',1000,-0.2,0.2) diffhistos['meanU']=ROOT.TH1D('dmeanU','meanU Diff',1000,-0.02,0.02) diffhistos['meanV']=ROOT.TH1D('dmeanV','meanV Diff',1000,-0.02,0.02) diffhistos['meanErrU']=ROOT.TH1D('dmeanErrU','mean Err U Diff',1000,-10,10) diffhistos['meanErrV']=ROOT.TH1D('dmeanErrV','mean Err V Diff',1000,-10,10) diffcanvas=ROOT.TCanvas('Differences','Differences',1000,1000) diffcanvas.Divide(3,3) errhistos={} errhistos['meanErrU2D']=ROOT.TH1D('MeanErrU1','Mean Err U 2D',1000,-0.1,0.1) errhistos['meanErrV2D']=ROOT.TH1D('MeanErrV1','Mean Err V 2D',1000,-0.1,0.1) errhistos['meanErrUMulti']=ROOT.TH1D('MeanErrU2','Mean Err U Multi',1000,-0.1,0.1) errhistos['meanErrVMulti']=ROOT.TH1D('MeanErrV2','Mean Err V Multi',1000,-0.1,0.1) errcanvas=ROOT.TCanvas('Errors','Errors',1000,1000) errcanvas.Divide(2,2) sigmahistos={} sigmahistos['U2D']=ROOT.TH1D('sigmaU2D','sigma U 2D',1000,0,1) sigmahistos['V2D']=ROOT.TH1D('sigmaV2D','sigma V 2D',1000,0,1) sigmahistos['UMulti']=ROOT.TH1D('sigmaUMulti','sigma U Multi',1000,0,1) sigmahistos['VMulti']=ROOT.TH1D('sigmaVMulti','sigma V Multi',1000,0,1) sigmacanvas=ROOT.TCanvas('sigmas','sigmas',1000,1000) sigmacanvas.Divide(2,2) hdproj_clu=[] dproj_clu=[] hdproj_clu.append(ROOT.TH1D('hdproj_clu1','Digis Projected Axis1 Amp (cluster)',50,-2,2)) hdproj_clu.append(ROOT.TH1D('hdproj_clu2','Digis Projected Axis2 Amp (cluster)',50,-2,2)) hdproj_clu.append(ROOT.TH1D('hdproj_clu1_n','Digis Projected Axis1 (cluster)',50,-2,2)) hdproj_clu.append(ROOT.TH1D('hdproj_clu2_n','Digis Projected Axis2 (cluster)',50,-2,2)) hdproj_clu.append(ROOT.TH1D('hdproj_clu1_norm','Digis Projected Axis1 Normalized(cluster)',50,-2,2)) hdproj_clu.append(ROOT.TH1D('hdproj_clu2_norm','Digis Projected Axis2 Normalized(cluster)',50,-2,2)) hdproj_track=[] dproj_track=[] hdproj_track.append(ROOT.TH1D('hdproj_tr1','Digis Projected Axis1 Amp(track)',100,-2,2)) hdproj_track.append(ROOT.TH1D('hdproj_tr2','Digis Projected Axis2 Amp(track)',100,-2,2)) hdproj_track.append(ROOT.TH1D('hdproj_tr1_n','Digis Projected Axis1 (track)',100,-2,2)) hdproj_track.append(ROOT.TH1D('hdproj_tr2_n','Digis Projected Axis2 (track)',100,-2,2)) hdproj_track.append(ROOT.TH1D('hdproj_tr1_norm','Digis Projected Axis1 Normalized(track)',100,-2,2)) hdproj_track.append(ROOT.TH1D('hdproj_tr2_norm','Digis Projected Axis2 Normalized(track)',100,-2,2)) hdproj_glob=[] dproj_glob=[] hdproj_glob.append(ROOT.TH1D('hdproj_gl1','Digis Projected Axis1 Amp(global)',1000,-2,2)) hdproj_glob.append(ROOT.TH1D('hdproj_gl2','Digis Projected Axis2 Amp(global)',1000,-2,2)) hdproj_glob.append(ROOT.TH1D('hdproj_gl1_n','Digis Projected Axis1 (global)',1000,-2,2)) hdproj_glob.append(ROOT.TH1D('hdproj_gl2_n','Digis Projected Axis2 (global)',1000,-2,2)) hdproj_glob.append(ROOT.TH1D('hdproj_gl1_norm','Digis Projected Axis1 Normalized(global)',1000,-2,2)) hdproj_glob.append(ROOT.TH1D('hdproj_gl2_norm','Digis Projected Axis2 Normalized(global)',1000,-2,2)) digiProjCanvas=ROOT.TCanvas('digiproj','Digis Projected on Axis',1500,1000) digiProjCanvas.Divide(2,3) canv3d=ROOT.TCanvas('canv3d','3D',1000,1000) u='' #create result ntuple resultTuple=ROOT.TNtuple("Result",'','dmeanU:dmeanV:dmeanErrU:dmeanErrV:meanErrU2D:meanErrV2D:meanErrUMult:meanErrVMult:sigU2D:sigV2D:sigUMult:sigVMult:status2D:statusMult:qual2D:qualMult:meanU2D:meanV2D:meanUMult:meanVMult:rho2D:rhoMult:mcU:mcV:ftime2D:ftimeMult:cogU:cogV') fit2dgauss=ROOT.TF2("fitGauss",ROOT.gauss2d,-10,10,-10,10,6) for e in tree: ev+=1 if args.events!=-1 and ev==args.events: if args.pfile!='None': canv.SaveAs(args.pfile+']') break if args.goto!=-1 and args.goto>ev: continue print "EVENT:",ev trackcount=-1 #for tr in e.TrackFitStat_0: for tr in e.CutCosmics: trackcount+=1 skipTrack=False print 'track',trackcount Umean=ROOT.TVector3(0,0,0) Vmean=ROOT.TVector3(0,0,0) U=tr.GetMidPlaneU() V=tr.GetMidPlaneV() O=tr.GetMidPlaneO() trackPoints=tr.GetProjectionPoints() ids=tr.GetHitIDs() first=True print tr.nhits(),'cluster in track' for h in hdproj_track: h.Reset() for c in range(tr.nhits()): print "CLUSTER:",c hit=e.TpcSPHit.At(ids[c]) cluster=e.TpcCluster.At(hit.getClusterID()) mcPos=tr.GetMCRefPos(c) #print 'cluster has size: 2D=',cluster.size2D(),'3D=',cluster.size() if cluster.size2D()<3: continue clPos=ROOT.TVector3(tr.GetHitPosition(c)-O[c]) if clPos.Z()args.zcut[1]: continue #clPos=clPos-O[c] mcPos=mcPos-O[c] clPosUV=ROOT.TVector3(clPos*U[c]/U[c].Mag(),clPos*V[c]/V[c].Mag(),0) mcPosUV=ROOT.TVector3(mcPos*U[c]/U[c].Mag(),mcPos*V[c]/V[c].Mag(),0) trPos=trackPoints[c] clAmp=tr.GetAmp(c) ppos=ROOT.TPolyMarker() ppos.SetPoint(0,clPosUV.X(),clPosUV.Y()) ppos.SetMarkerStyle(20) graph2D=ROOT.TGraph2D() graph2D.SetTitle('Digis On Plane') graph1=ROOT.TGraph() graph1.SetTitle('Projection1') graph2=ROOT.TGraph() graph2.SetTitle('Projection2') digiTuple=ROOT.TNtuple('UVdigipos','','u:v:amp') digisOnPlane=tr.GetDigisOnPlane(c) digisOnPlaneAmp=tr.GetDigisOnPlaneAmp(c) cov=ROOT.TMatrixD(2,2) cov[0][0]=0 cov[1][1]=0 cov[1][0]=cov[0][1]=0 digiAmpSum=0 cogUV=ROOT.TVector3(0,0,0) for idigi in range(len(digisOnPlane)): digipos=digisOnPlane[idigi] digiamp=digisOnPlaneAmp[idigi] cogUV.SetX(cogUV.X()+(digiamp*digipos.X())) cogUV.SetY(cogUV.Y()+(digiamp*digipos.Y())) digiAmpSum+=digiamp cogUV*=1./digiAmpSum pdigipos=ROOT.TPolyMarker() for idigi in range(len(digisOnPlane)): digipos=digisOnPlane[idigi] digiamp=digisOnPlaneAmp[idigi] digiTuple.Fill(digipos.X(),digipos.Y(),digiamp) graph1.SetPoint(idigi,digipos.X(),digiamp) graph2.SetPoint(idigi,digipos.Y(),digiamp) graph2D.SetPoint(idigi,digipos.X(),digipos.Y(),digiamp) pdigipos.SetPoint(idigi,digipos.X(),digipos.Y()) cm=ROOT.TMatrixD(2,1) cm[0][0]=digipos.X()-cogUV.X() cm[1][0]=digipos.Y()-cogUV.Y() cm_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,cm) acov=ROOT.TMatrixD(cm,ROOT.TMatrixD.kMult,cm_t) acov*=digiamp cov+=acov if args.watchthecov: drawCov(cov,digiAmpSum,digiTuple,ctemp) cov*=1./clAmp cov=ROOT.TMatrixDSym(2,cov.GetMatrixArray()) #get eigenvectors and eigenvalues of cov eigenproblem=ROOT.TMatrixDEigen(ROOT.TMatrixD(cov)) eigenwert=eigenproblem.GetEigenValues() eigenvector=eigenproblem.GetEigenVectors() if cov[0][0]==0 and cov[1][1]==0: skipTrack=True print 'skipping track due to cov!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' continue varset=ROOT.RooArgSet('digicoords') ucogmin=umins=umin=digiTuple.GetMinimum('u') ucogmax=umaxs=umax=digiTuple.GetMaximum('u') vcogmin=vmins=vmin=digiTuple.GetMinimum('v') vcogmax=vmaxs=vmax=digiTuple.GetMaximum('v') maxamp=digiTuple.GetMaximum('amp') if args.incrRange: umins*=multIf(umin,0.5,1.5) umaxs*=multIf(umax,1.5,0.5) vmins*=multIf(vmin,0.5,1.5) vmaxs*=multIf(vmax,1.5,0.5) ucogmin=cogUV.X()-math.sqrt(cov[0][0]) ucogmax=cogUV.X()+math.sqrt(cov[0][0]) vcogmin=cogUV.Y()-math.sqrt(cov[1][1]) vcogmax=cogUV.Y()+math.sqrt(cov[1][1]) #calculate digi error UVplane=ROOT.TMatrixD(3,2); for i in range(3): UVplane[i][0]=U[c][i]; UVplane[i][1]=V[c][i]; digiErr=ROOT.TVector2(0,0) digiCov=ROOT.TMatrixDSym(3) digiCov[0][1]=digiCov[0][2]=digiCov[1][0]=digiCov[2][0]=digiCov[1][2]=digiCov[2][1]=0; digiCov[0][0]=1./(0.2598/ROOT.TMath.Sqrt(12)); digiCov[1][1]=1./(0.3/ROOT.TMath.Sqrt(12)); digiCov[2][2]=1./(0.152867/ROOT.TMath.Sqrt(12)); #digiCov=digiCov.Invert() digiCovOnPlane=ROOT.TMatrixDSym(2) digiCovOnPlane=digiCov.SimilarityT(UVplane) digiErr.Set(1./digiCovOnPlane[0][0],1./digiCovOnPlane[1][1]) print 'error of digis' digiErr.Print() #fill tgraph graph=ROOT.TGraph2DErrors(digiTuple.GetEntries()) graph.Clear() graph.SetMarkerStyle(20) pcounter=-1 for tupcl in digiTuple: pcounter+=1 graph.SetPoint(pcounter,tupcl.u,tupcl.v,tupcl.amp) graph.SetPointError(pcounter,digiErr.X(),digiErr.Y(),clAmp*0.025) #graph.SetPointError(pcounter,0,0,tupcl.amp*0.05) print pcounter,tupcl.u,tupcl.v,tupcl.amp #graph=ROOT.TGraph2DErrors(len(ex),px,py,pz,ex,ey,ez) #graph.SetMarkerStyle(20) fit2dgauss.SetParName(0,"sigU") fit2dgauss.SetParameter(0,math.sqrt(cov[0][0])) fit2dgauss.FixParameter(0,math.sqrt(cov[0][0])) fit2dgauss.SetParLimits(0,0.5*math.sqrt(cov[0][0]),1.2*math.sqrt(cov[0][0])) fit2dgauss.SetParName(1,"sigV") fit2dgauss.SetParameter(1,math.sqrt(cov[1][1])) fit2dgauss.FixParameter(1,math.sqrt(cov[1][1])) fit2dgauss.SetParLimits(1,0.5*math.sqrt(cov[1][1]),1.2*math.sqrt(cov[1][1])) fit2dgauss.SetParName(2,"rho") if math.sqrt(cov[0][0]*cov[1][1]) != 0: fit2dgauss.SetParameter(2,cov[1][0]/math.sqrt(cov[0][0]*cov[1][1])) fit2dgauss.FixParameter(2,cov[1][0]/math.sqrt(cov[0][0]*cov[1][1])) # fit2dgauss.SetParLimits(2, #0.9*cov[1][0]/math.sqrt(cov[0][0]*cov[1][1]), #1.1*cov[1][0]/math.sqrt(cov[0][0]*cov[1][1])) else: fit2dgauss.SetParameter(2,0) fit2dgauss.SetParLimits(2,-1,1) fit2dgauss.SetRange(umins,vmins,umaxs,vmaxs) fit2dgauss.SetParName(3,"MeanU") fit2dgauss.SetParameter(3,cogUV.X()) fit2dgauss.SetParLimits(3,cogUV.X()-1*math.sqrt(cov[0][0]),cogUV.X()+1*math.sqrt(cov[0][0])) fit2dgauss.SetParName(4,"MeanV") fit2dgauss.SetParameter(4,cogUV.Y()) fit2dgauss.SetParLimits(4,cogUV.Y()-1*math.sqrt(cov[1][1]),cogUV.Y()+1*math.sqrt(cov[1][1])) fit2dgauss.SetParName(5,"Amp") fit2dgauss.SetParameter(5,maxamp) #fit2dgauss.SetParLimits(5,0.5*maxamp,1.2*maxamp) x,y,z=ROOT.Double(0),ROOT.Double(0),ROOT.Double(0) xarr=graph.GetX() yarr=graph.GetY() zarr=graph.GetZ() print '' for iii in range(graph.GetN()): print iii, xarr[iii],yarr[iii],zarr[iii] result=graph.Fit(fit2dgauss,'S E M') #result=digiTuple.UnbinnedFit("fitGauss",'u:v:amp') resl_result=ROOT.TFitResultPtr.Get(result) real_result=ROOT.TFitResultPtr(result) fitcov=ROOT.TMatrixDSym(6) fitcov=real_result.GetCovarianceMatrix() # print 'the cov i got' fitcov.Print() # print 'the cov from print' # real_result.PrintCovMatrix() fitcorr=real_result.GetCorrelationMatrix() fitcorr.Print() fiteigenproblem=ROOT.TMatrixDSymEigen(fitcov) fiteigenvector=fiteigenproblem.GetEigenVectors() fiteigenwert=fiteigenproblem.GetEigenValues() fitter = ROOT.TVirtualFitter.GetFitter() fitcov2=ROOT.TMatrixD(6,6,fitter.GetCovarianceMatrix()) fitcov2.Print() #fiteigenwert.Print() print 'errors:',real_result.NormalizedErrors() for iii in range(6): print 'sqrt(cov):', math.sqrt(abs(fitcov[iii][iii])),'sqrt(cov2)',math.sqrt(abs(fitcov2[iii][iii])) print 'parerror:',result.ParError(iii),math.sqrt(real_result.CovMatrix(iii,iii)) fit2dgauss.SetNpy(100) fit2dgauss.SetNpx(100) canv3d.Clear() canv3d.cd() fit2dgauss.Draw('surf4') graph.Draw('same pcol err') print '' for iii in range(graph.GetN()): print iii, xarr[iii],yarr[iii],zarr[iii] canv3d.Update() #u=raw_input('pause') start=time.time() runtime2=time.time()-start fittime2.Fill(runtime2) sigX=fit2dgauss.GetParameter(0) sigY=fit2dgauss.GetParameter(1) rho=fit2dgauss.GetParameter(2) muX=fit2dgauss.GetParameter(3) muY=fit2dgauss.GetParameter(4) muXError=fit2dgauss.GetParError(3) muYError=fit2dgauss.GetParError(4) #diffhistos['fittime'].Fill(runtime-runtime2) if (sigX!=0 or cov[0][0]!=0): diffhistos['sigmaU'].Fill( (sigX-math.sqrt(cov[0][0])) / ( 0.5 * (sigX+math.sqrt(cov[0][0])) ) ) if (sigY!=0 or cov[1][1]!=0): diffhistos['sigmaV'].Fill( (sigY-math.sqrt(cov[1][1])) / ( 0.5 * (sigY+math.sqrt(cov[1][1])) ) ) if math.sqrt(cov[0][0]*cov[1][1]) !=0: meanrho=0.5*(rho+cov[1][0]/math.sqrt(cov[0][0]*cov[1][1])) if meanrho!=0: diffhistos['rho'].Fill((rho-cov[1][0]/math.sqrt(cov[0][0]*cov[1][1]))/meanrho) else: meanrho=rho dmeanU=9999 dmeanV=9999 dmeanErrU=9999 dmeanErrV=9999 ''' if (muX!=0 or muX2!=0): dmeanU = ( muX-muX2 ) dmeanU/= ( 0.5*(muX+muX2) ) diffhistos['meanU'].Fill( dmeanU ) if muY!=0 or muY2!=0: dmeanV = (muY-muY2) dmeanV/= ( 0.5*(muY+muY2) ) diffhistos['meanV'].Fill( dmeanV ) if (0.5*(muX.getError()+muX2.getError())!=0): dmeanErrU =(muX.getError()-muX2.getError()) dmeanErrU/=(0.5*(muX.getError()+muX2.getError())) diffhistos['meanErrU'].Fill( dmeanErrU ) if (0.5*(muY.getError()+muY2.getError())!=0): dmeanErrV =(muY.getError()-muY2.getError()) dmeanErrV/=(0.5*(muY.getError()+muY2.getError())) diffhistos['meanErrV'].Fill( dmeanErrV ) ''' errhistos['meanErrU2D'].Fill(muXError) errhistos['meanErrV2D'].Fill(muYError) #errhistos['meanErrUMulti'].Fill(muX2.getError()) #errhistos['meanErrVMulti'].Fill(muY2.getError()) sigmahistos['U2D'].Fill(sigX) sigmahistos['V2D'].Fill(sigY) sigmahistos['UMulti'].Fill(math.sqrt(cov[0][0])) sigmahistos['VMulti'].Fill(math.sqrt(cov[1][1])) for h in hdproj_clu: h.Reset() dpos_tmp=fillAxisProjHist(digiTuple,eigenvector,cogUV,hdproj_clu,hdproj_track,hdproj_glob) ''' for dp in dpos_tmp: dproj_clu.append(ROOT.TVector2(dp)) dproj_track.append(ROOT.TVector2(dp)) dproj_glob.append(ROOT.TVector2(dp)) ''' digiProjCanvas.cd(1) hdproj_clu[4].Divide(hdproj_clu[0], hdproj_clu[2]) hdproj_clu[0].Draw() digiProjCanvas.cd(2) hdproj_clu[5].Divide(hdproj_clu[1], hdproj_clu[3]) hdproj_clu[1].Draw() digiProjCanvas.cd(3) hdproj_track[0].Draw() digiProjCanvas.cd(4) hdproj_track[1].Draw() digiProjCanvas.cd(5) hdproj_glob[0].Draw() digiProjCanvas.cd(6) hdproj_glob[1].Draw() digiProjCanvas.Update() setGraph(graph2D) rho2D=-100 if cov[0][0]*cov[1][1] !=0: rho2D=cov[1][0]/math.sqrt(cov[0][0]*cov[1][1]) tofill=[dmeanU, dmeanV, dmeanErrU, dmeanErrV, muXError, muYError, #muX2.getError(), #muY2.getError(), sigX, sigY, math.sqrt(cov[0][0]), math.sqrt(cov[1][1]), #status2d, #mvgResult.status(), #covQual2d, #mvgResult.covQual(), muX, muY, #muX2, #muY2, rho, rho2D, mcPosUV.X(), mcPosUV.Y(), runtime2, #runtime, cogUV.X(), cogUV.Y()] #tofill #resultTuple.Fill(array.array("f",tofill)) if not args.batch: digihist=ROOT.TH2D('digihist','digis',100,umins,umaxs,100,vmins,vmaxs) canv.cd(1) digihist.Draw() pdigipos.SetMarkerStyle(20) pdigipos.SetMarkerColor(ROOT.kAzure) pdigipos.Draw('same') teigenvector=ROOT.TVector2(eigenvector[0][0],eigenvector[1][0]) teigenvector2=ROOT.TVector2(eigenvector[0][1],eigenvector[1][1]) eigenvector_i=ROOT.TMatrixD(ROOT.TMatrixD.kInverted,eigenvector) pcov=ROOT.TEllipse(clPosUV.X(),clPosUV.Y(),math.sqrt(abs(eigenwert[0][0])),math.sqrt(abs(eigenwert[1][1])),0,360,ROOT.TMath.RadToDeg()*teigenvector.Phi()) pcov.SetFillStyle(4000) pcov.Draw("same") ppos.Draw('same') mcpos=ROOT.TPolyMarker() mcpos.SetPoint(0,mcPosUV.X(),mcPosUV.Y()) mcpos.SetMarkerStyle(29) mcpos.SetMarkerColor(ROOT.kMagenta) mcpos.Draw('same') pnewpos=ROOT.TPolyMarker() pnewpos.SetPoint(0,muX,muY) pnewpos.SetMarkerStyle(21) pnewpos.SetMarkerColor(ROOT.kPink) pnewpos.Draw('same') canv.cd(2) canv.cd(3) fittime.Draw() canv.cd(4) canv.cd(4).Clear() #hg2D.Draw('colz') text=[] text.append('Cluster pos:') text.append(str(clPosUV.X())+' '+str(clPosUV.Y())) text.append('New Pos:') text.append(str(muX)+' '+str(muY)) text.append('Difference%: {0:8.4f} {1:8.4f}'.format( 100*abs(muX-clPosUV.X())/clPosUV.X() ,100*abs(muY-clPosUV.Y())/clPosUV.Y() )) text.append('MC res: {0:6.3f}'.format(math.sqrt( (muX- mcPosUV.X())**2 + (muY- mcPosUV.Y())**2 ) )) text.append('error: U={0:6.3f} V={1:6.3f} correlation={2:6.3f}'.format(muXError,muYError,rho2D)) text.append('pullU: {0:6.3f} pullV: {1:6.3f}'.format( (muX-mcPosUV.X())/muXError , (muY- mcPosUV.Y())/muYError)) drawText(text) canv.Update() if args.pfile!='None': canv.SaveAs(args.pfile) tr.GetHitPosition(c).Print() newposU=ROOT.TVector3(U[c]) newposU*=clPosUV.X() newposV=ROOT.TVector3(V[c]) newposV*=clPosUV.Y() newpos=newposU+newposV+O[c] newpos.Print() if args.events==-1: u=raw_input('next cluster?') if u=='q': if args.pfile!='None': canv.SaveAs(args.pfile+']') break if u=='n': if args.pfile!='None': canv.SaveAs(args.pfile+']') break if args.pfile!='None': canv.SaveAs(args.pfile) if args.events==-1 and not skipTrack: if u=='q': if args.pfile!='None': canv.SaveAs(args.pfile+']') break u=raw_input('next track?') if u=='n': if args.pfile!='None': canv.SaveAs(args.pfile+']') break if u=='q' or u=='n': break canv.cd(4).Clear() #rfile=ROOT.TFile(args.outfile,'recreate') rfile.cd() diffcanvas.cd(1) diffhistos['sigmaU'].Draw() diffcanvas.cd(2) diffhistos['sigmaV'].Draw() diffcanvas.cd(3) diffhistos['fittime'].Draw() diffcanvas.cd(4) diffhistos['meanU'].Draw() diffcanvas.cd(5) diffhistos['meanV'].Draw() diffcanvas.cd(6) diffhistos['rho'].Draw() diffcanvas.cd(7) diffhistos['meanErrU'].Draw() diffcanvas.cd(8) diffhistos['meanErrV'].Draw() for h in diffhistos: diffhistos[h].Write() errcanvas.cd(1) errhistos['meanErrU2D'].Draw() errcanvas.cd(2) errhistos['meanErrV2D'].Draw() errcanvas.cd(3) errhistos['meanErrUMulti'].Draw() errcanvas.cd(4) errhistos['meanErrVMulti'].Draw() for h in errhistos: errhistos[h].Write() sigmacanvas.cd(1) sigmahistos['U2D'].Draw() sigmacanvas.cd(2) sigmahistos['V2D'].Draw() sigmacanvas.cd(3) sigmahistos['UMulti'].Draw() sigmacanvas.cd(4) sigmahistos['VMulti'].Draw() for h in sigmahistos: sigmahistos[h].Write() ''' counter=-1 for h in diffhistos: counter+=1 diffcanvas.cd(counter+1) diffhistos[h].Draw() counter=-1 for h in errhistos: counter+=1 errcanvas.cd(counter+1) errhistos[h].Draw() ''' diffcanvas.Update() errcanvas.Update() sigmacanvas.Update() if not args.batch: diffcanvas.Write() errcanvas.Write() sigmacanvas.Write() if args.pfile2!='None': diffcanvas.SaveAs(args.pfile2+'(') errcanvas.SaveAs(args.pfile2) sigmacanvas.SaveAs(args.pfile2+')') resultTuple.Write() rfile.Close() while u!='q' and not args.batch: u=raw_input('done?')