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 calcPOCA(thisrefpos,mom,pos): tmp1=ROOT.TVector3(mom) tmp1.SetMag(1) tmp2=ROOT.TVector3(pos-thisrefpos) fac=tmp2*mom fac/=mom.Mag() tmp1*= fac poca=ROOT.TVector3(thisrefpos+tmp1) return poca def calcAmpWeights(digipos,digiamp,mom): change=1 weights=[] sumamp=0 startpos=ROOT.TVector3(0,0,0) for i in range(len(digipos)): tmp=ROOT.TVector3(digipos[i]) #tmp*=abs(digiamp[i]) startpos+=tmp sumamp+=abs(digiamp[i]) #startpos*=1./sumamp startpos*=1./len(digipos) sumamp_unw=sumamp while (abs(change)>0.001): lastchange=change weights=[] sumamp_new=0 startpos_new=ROOT.TVector3(0,0,0) for idigi in range(len(digipos)): poca=calcPOCA(startpos,mom,digipos[idigi]) dirToPad=digipos[idigi]-poca if dirToPad.Mag()==0: weight=0 else: weight=padplane.GetPad(0).GetWeight(abs(digiamp[idigi])/sumamp_unw, dirToPad.X(), dirToPad.Y()) if weight==0: weights.append(1) else: weights.append(1./weight) tmp=ROOT.TVector3(digipos[idigi]) tmp*=weights[-1]*abs(digiamp[idigi]) startpos_new+=tmp sumamp_new+=weights[-1]*abs(digiamp[idigi]) startpos_new*=1./sumamp_new tmp=ROOT.TVector3(startpos_new-startpos) tmp*=0.5 startpos+=tmp sumamp=sumamp_new change=(startpos-startpos_new).Mag() return weights parser=argparse.ArgumentParser(description="plots digisposition projected on detector plane") parser.add_argument('recofile',help='the reco file',type=str,default='./') parser.add_argument('digifile',help='the digifile',type=str) #parser.add_argument('mcfile',help='the mc file',type=str) 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('--padf',help='the padplane file to use',type=str,default='tpc/parfiles/padPlane_FOPI.dat') parser.add_argument('--pshape',help='pad shape file',type=str,default='tpc/parfiles/PadShapePrototype.dat') parser.add_argument('--padres',help='padresponse file anme',type=str,default='tpc/FOPI/par/padresponse.txt') 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.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() tree.AddFriend('cbmsim',args.digifile) parfile=args.digifile.replace('raw.root','param.root') rparfile=ROOT.TFile(parfile,'read') digipar=rparfile.Get('TpcDigiPar') digimapper=ROOT.TpcDigiMapper.getInstance() digimapper.init(digipar) padShapePool=ROOT.TpcPadShapePool(args.pshape,args.padres) padplane=ROOT.TpcPadPlane(args.padf,padShapePool) 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) tmpAmp=ROOT.RooRealVar('tmpAmp','tmpAmp',0) landau=ROOT.RooLandau('landau','landau',tmpAmp,ROOT.RooFit.RooConst(1.79797e+03),ROOT.RooFit.RooConst(7.33996e+02)) 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') 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 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 'Event: {0}, Track: {1}, Cluster: {2}'.format(ev,trackcount,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 #calculate digi error (padsize) UVplane=ROOT.TMatrixD(3,2); for i in range(3): UVplane[i][0]=U[c][i] UVplane[i][1]=V[c][i] UVplane_inv=ROOT.TMatrixD(UVplane) UVplane_inv.Invert() 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]=0.#1./0.2598 digiCov[1][1]=0.#1./0.3 digiCov[2][2]=1./0.152867 digiCovOnPlane=digiCov.SimilarityT(UVplane); offsetU=0 if digiCovOnPlane[0][0]!=0: offsetU=1./digiCovOnPlane[0][0] offsetV=0 if digiCovOnPlane[1][1]!=0: offsetV=1./digiCovOnPlane[1][1] digiErr.Set(offsetU,offsetV); print 'the found error on plane: {0} {1}'.format(digiErr.X(),digiErr.Y()) #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) digiTuple=ROOT.TNtuple('UVdigipos','','u:v:amp') digisOnPlane=tr.GetDigisOnPlane(c) digisOnPlaneAmp=tr.GetDigisOnPlaneAmp(c) digisMCPos=tr.GetDigiMCPos(c) digisPos=tr.GetDigiPos(c) digisID=tr.GetDigiID(c) refmom=tr.GetProjectionPointsMom(c) cov=ROOT.TMatrixD(2,2) cov[0][0]=0 cov[1][1]=0 cov[1][0]=cov[0][1]=0 digiAmpSum=0 weights=calcAmpWeights(digisPos,digisOnPlaneAmp,refmom) #calculate the center of gravity cogUV=ROOT.TVector3(0,0,0) maxdigiamp=0 weighted_amp_sum=0 for idigi in range(len(digisOnPlane)): digipos=digisOnPlane[idigi] digiamp=abs(digisOnPlaneAmp[idigi]) weighted_amp_sum+=digiamp*weights[idigi] print 'digiamp:{0}, posU:{1}, posV:{2}'.format(digiamp,digipos.X(),digipos.Y()) print 'xyz pos= ({0},{1},{2})'.format(digisPos[idigi].X(),digisPos[idigi].Y(),digisPos[idigi].Z()) cogUV.SetX(cogUV.X()+(weights[idigi]*digiamp*digipos.X())) cogUV.SetY(cogUV.Y()+(weights[idigi]*digiamp*digipos.Y())) digiAmpSum+=digiamp maxdigiamp=max(maxdigiamp,abs(digiamp)) #cogUV*=1./digiAmpSum cogUV*=1./weighted_amp_sum ''' wgaus=ROOT.TF1('f','gaus',-10,10) wgaus.SetParameter(0,con) wgaus.SetParameter(1,0) wgaus.SetParameter(2,sig) ''' pdigiMC=ROOT.TPolyMarker() for idigi in range(len(digisOnPlane)): digipos=digisOnPlane[idigi] digiamp=abs(digisOnPlaneAmp[idigi]) digiweight=weights[idigi] digiID=digisID[idigi] thedigi=e.TpcDigi.At(digiID) samplemap=thedigi.getSampleMap() tmp1=ROOT.TVector3(U[c]) tmp2=ROOT.TVector3(V[c]) tmp1*=digipos.X() tmp2*=digipos.Y() digipos3D=ROOT.TVector3(tmp1+tmp2+O[c]) tmp3=ROOT.TVector3(tr.GetHitPosition(c)) dirToPad=digipos3D-tmp3 #digiweight=padplane.GetPad(0).GetWeight(digiamp/maxdigiamp,dirToPad.X(),dirToPad.Y()) digiweight=weights[idigi] print 'damp={0}, max={1}, weighted Amp={2}, weight={3}'.format(digiamp,maxdigiamp,digiamp*digiweight,digiweight) print 'U=({0},{1},{2}), V=({3},{4},{5})'.format(U[c].X(),U[c].Y(),U[c].Z(),V[c].X(),V[c].Y(),V[c].Z()) #digiamp*=digiweight for sam in samplemap: sample=e.TpcSample.At(sam.first) sampos=ROOT.TVector3(0,0,0) digimapper.map(sample,sampos) sampos.Print() samposUV=sampos-O[c] samposUV=SetXYZ(samposUV*U[c]/U[c].Mag(),samposUV*V[c]/V[c].Mag(),0) sresX=samposUV.X-cogUV.X() digiMC=digisMCPos[idigi]-O[c] pdigiMC.SetPoint(idigi,digiMC*U[c]/U[c].Mag(),digiMC*V[c]/V[c].Mag()) spreadweight1=1+ U[c].Perp() * ( (0.0343738+0.0286046*ROOT.TMath.Sqrt(O[c].Z())) -1); spreadweight2=1+ V[c].Perp() * ( (0.0343738+0.0286046*ROOT.TMath.Sqrt(O[c].Z())) -1); cm=ROOT.TMatrixD(2,1) mresX=digipos.X()-cogUV.X() mresY=digipos.Y()-cogUV.Y() print 'before: mresX={0}, mresY={1}'.format(mresX,mresY) ''' if mresX<0: mresX-=digiErr.X()/2. else: mresX+=digiErr.X()/2. if mresY<0: mresY-=digiErr.Y()/2. else: mresY+=digiErr.Y()/2. ''' cm[0][0]=mresX#*math.sqrt(spreadweight1) cm[1][0]=mresY#*math.sqrt(spreadweight2) print 'after: mresX={0}, mresY={1}'.format(mresX,mresY) cm_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,cm) acov=ROOT.TMatrixD(cm,ROOT.TMatrixD.kMult,cm_t) acov*=digiamp*digiweight #print 'digiamp',digiamp #acov.Print() #tmpAmp.setVal(digiamp) #print 'landauval={0}, amp={1}'.format(landau.getVal(),digiamp) #digiamp=4.30779e-01+6.67404e-02* math.log(digiamp) digiTuple.Fill(digipos.X(),digipos.Y(),digiamp*digiweight) #if digiweight==0: # digiweight=digiamp/maxdigiamp #digiamp/=digiweight cov+=acov if args.watchthecov: drawCov(cov,digiAmpSum,digiTuple,ctemp) #cov*=1./clAmp cov*=1./(weighted_amp_sum) cov.Print() #cov*=1./len(digisOnPlane) 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') 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) ''' umins-=digiErr.X()/2. umaxs+=digiErr.X()/2. vmins-=digiErr.Y()/2. vmaxs+=digiErr.Y()/2. ''' print 'the ranges used are: u=({0},{1}) v=({2},{3})'.format(umins,umaxs,vmins,vmaxs) 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]) varU=ROOT.RooRealVar('u','',umins,umaxs) #varU.setBins(40) varV=ROOT.RooRealVar('v','',vmins,vmaxs) #varV.setBins(40) varAmp=ROOT.RooRealVar('amp','',0,100000000) varset.add(varU) varset.add(varV) varset.add(varAmp) dataset=ROOT.RooDataSet('dataSet','digisOnPlane',digiTuple,varset,'','amp') #dataset=ROOT.RooDataSet('dataSet','digisOnPlane',digiTuple,varset,'') #datasetProj1=ROOT.RooDataSet('dataSet1','digis on Proj1',digiTuple,ROOT.RooArgSet(varU,varAmp),'','amp') #datasetProj2=ROOT.RooDataSet('dataSet2','digis on Proj2',digiTuple,ROOT.RooArgSet(varV,varAmp),'','amp') muX=ROOT.RooRealVar('muX','MeanX of 2D gaussian',0,ucogmin,ucogmax) muY=ROOT.RooRealVar('muY','MeanY of 2D gaussian',0,vcogmin,vcogmax) sigX=ROOT.RooRealVar('sigX','SigmaX of 2D gaussian',math.sqrt(cov[0][0]),0,5*math.sqrt(cov[0][0])) sigY=ROOT.RooRealVar('sigY','SigmaY of 2D gaussian',math.sqrt(cov[1][1]),0,5*math.sqrt(cov[1][1])) if math.sqrt(cov[0][0]*cov[1][1]) != 0: rho=ROOT.RooRealVar('rho','Correlation factor',cov[1][0]/math.sqrt(cov[0][0]*cov[1][1]),-1,1) else: rho=ROOT.RooRealVar('rho','Correlation factor',0,-1,1) rooGauss2D=ROOT.Roo2DimGauss('gauss2D','gauss 2D',varU,varV,rho,sigX,sigY,muX,muY) xVec=ROOT.RooArgList() xVec.add(varU) xVec.add(varV) muVec=ROOT.RooArgList() muX2=ROOT.RooRealVar('muX2','MeanX of Multivar gaussian',cogUV.X(),ucogmin,ucogmax) muY2=ROOT.RooRealVar('muY2','MeanY of Multivar gaussian',cogUV.Y(),vcogmin,vcogmax) muVec.add(muX2) muVec.add(muY2) #print '********************fitting mvg*****************' #for ghliju in digiTuple: # print 'U=',digiTuple.u,'V=',digiTuple.v,'amp=',digiTuple.amp sigma_resolU=ROOT.RooRealVar('sigma_resolU','sigma resolution U',0.1,0,0.15) sigma_resolV=ROOT.RooRealVar('sigma_resolV','sigma resolution V',0.1,0,0.15) mean_resolU=ROOT.RooRealVar('mean_resolU','mean resolution U',0) mean_resolV=ROOT.RooRealVar('mean_resolV','mean resolution V',0) s1=ROOT.RooRealVar('s1','resolution',3,0,20) s2=ROOT.RooRealVar('s2','resolution',3,0,20) mvg=ROOT.RooMultiVarGaussian("mvg","mvg",xVec,muVec,cov) gauss_relU=ROOT.RooGaussModel('gauss_relU','gauss resolutionU',varU,varU,s1,sigma_resolU) gauss_relV=ROOT.RooGaussModel('gauss_relV','gauss resolutionV',varV,varV,s2,sigma_resolV) gauss_rel=ROOT.RooProdPdf('gauss_rel','gauss resolution',ROOT.RooArgList(gauss_relU,gauss_relV)) prod_mvg=ROOT.RooProdPdf("prod_mvg","mvg*gauss_rel",ROOT.RooArgSet(mvg),ROOT.RooFit.Conditional(ROOT.RooArgSet(gauss_relU,gauss_relV),ROOT.RooArgSet(varU,varV))) start=time.time() mvgResult=mvg.fitTo(dataset, ROOT.RooFit.SumW2Error(True), ROOT.RooFit.Save(), ROOT.RooFit.Verbose(False), ROOT.RooFit.PrintLevel(-1),ROOT.RooFit.PrintEvalErrors(0), ROOT.RooFit.Warnings(False) ) #print '********************fit done********************' if (mvgResult.status()!=0): print 'RooFit mvg failed {0} {1}'.format(mvgResult.status(), mvgResult.covQual()) #continue fitCovMvg=mvgResult.covarianceMatrix() print 'mvg errors and cov' fitCovMvg.Print() print muX2.getError(),muY2.getError() print math.sqrt(fitCovMvg[0][0]),math.sqrt(fitCovMvg[1][1]) runtime=time.time()-start fittime.Fill(runtime) if not args.no2D: #print '********************fitting 2Dgauss*************' start=time.time() rooFitResGauss2D=rooGauss2D.fitTo(dataset, ROOT.RooFit.SumW2Error(True), ROOT.RooFit.Save(), ROOT.RooFit.Verbose(False), ROOT.RooFit.PrintLevel(-1), ROOT.RooFit.PrintEvalErrors(0), ROOT.RooFit.Warnings(False)) #print '********************fit done********************' if (rooFitResGauss2D.status()!=0 or rooFitResGauss2D.covQual()!=-1): print '2d gaussian fit failed' #continue status2d=rooFitResGauss2D.status() covQual2d=rooFitResGauss2D.covQual() else: status2d=0 covQual2d=0 runtime2=time.time()-start fittime2.Fill(runtime2) diffhistos['fittime'].Fill(runtime-runtime2) if (sigX.getVal()!=0 or cov[0][0]!=0): diffhistos['sigmaU'].Fill( (sigX.getVal()-math.sqrt(cov[0][0])) / ( 0.5 * (sigX.getVal()+math.sqrt(cov[0][0])) ) ) if (sigY.getVal()!=0 or cov[1][1]!=0): diffhistos['sigmaV'].Fill( (sigY.getVal()-math.sqrt(cov[1][1])) / ( 0.5 * (sigY.getVal()+math.sqrt(cov[1][1])) ) ) if math.sqrt(cov[0][0]*cov[1][1]) !=0: meanrho=0.5*(rho.getVal()+cov[1][0]/math.sqrt(cov[0][0]*cov[1][1])) if meanrho!=0: diffhistos['rho'].Fill((rho.getVal()-cov[1][0]/math.sqrt(cov[0][0]*cov[1][1]))/meanrho) else: meanrho=rho.getVal() dmeanU=9999 dmeanV=9999 dmeanErrU=9999 dmeanErrV=9999 if (muX.getVal()!=0 or muX2.getVal()!=0): dmeanU = ( muX.getVal()-muX2.getVal() ) dmeanU/= ( 0.5*(muX.getVal()+muX2.getVal()) ) diffhistos['meanU'].Fill( dmeanU ) if muY.getVal()!=0 or muY2.getVal()!=0: dmeanV = (muY.getVal()-muY2.getVal()) dmeanV/= ( 0.5*(muY.getVal()+muY2.getVal()) ) 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(muX.getError()) errhistos['meanErrV2D'].Fill(muY.getError()) errhistos['meanErrUMulti'].Fill(muX2.getError()) errhistos['meanErrVMulti'].Fill(muY2.getError()) sigmahistos['U2D'].Fill(sigX.getVal()) sigmahistos['V2D'].Fill(sigY.getVal()) 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() 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, muX.getError(), muY.getError(), muX2.getError(), muY2.getError(), sigX.getVal(), sigY.getVal(), math.sqrt(cov[0][0]), math.sqrt(cov[1][1]), status2d, mvgResult.status(), covQual2d, mvgResult.covQual(), muX.getVal(), muY.getVal(), muX2.getVal(), muY2.getVal(), rho.getVal(), rho2D, mcPosUV.X(), mcPosUV.Y(), runtime2, runtime, cogUV.X(), cogUV.Y()] #tofill resultTuple.Fill(array.array("f",tofill)) if not args.batch: hd=varU.createHistogram('digis on plane',varV,'') dataset.fillHistogram(hd,ROOT.RooArgList(varU,varV)) canv.cd(1) hd.SetStats(0) hd.Draw('colz') 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(cogUV.X(),cogUV.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=ROOT.TPolyMarker() ppos.SetPoint(0,clPosUV.X(),clPosUV.Y()) ppos.SetMarkerStyle(20) ppos.Draw('same') ''' cogpos=ROOT.TPolyMarker() cogpos.SetPoint(0,cogUV.X(),cogUV.Y()) cogpos.SetMarkerStyle(22) cogpos.Draw('same') mcpos=ROOT.TPolyMarker() mcpos.SetPoint(0,mcPosUV.X(),mcPosUV.Y()) mcpos.SetMarkerStyle(29) mcpos.SetMarkerColor(ROOT.kMagenta) mcpos.Draw('same') pdigiMC.SetMarkerStyle(4) pdigiMC.SetMarkerColor(1) pdigiMC.Draw('same') pnewpos=ROOT.TPolyMarker() print 'setting fit marker:',muX2.getVal(),muY2.getVal() pnewpos.SetPoint(0,muX2.getVal(),muY2.getVal()) pnewpos.SetMarkerStyle(21) pnewpos.SetMarkerColor(ROOT.kPink) pnewpos.Draw() canv.cd(2) hgmult=varU.createHistogram('digis on plane Mult',varV,'') hg2D=varU.createHistogram('digis on plane 2D',varV,'') rooGauss2D.fillHistogram(hg2D,ROOT.RooArgList(varU,varV)) mvg.fillHistogram(hgmult,ROOT.RooArgList(varU,varV)) hgmult.Draw('surf') hd.Draw('lego same') canv.cd(3) fittime.Draw() canv.cd(4) canv.cd(4).Clear() #hg2D.Draw('colz') cogRes= math.sqrt( (cogUV.X()- mcPosUV.X())**2 + (cogUV.Y()- mcPosUV.Y())**2 ) fitRes= math.sqrt( (muX2.getVal()- mcPosUV.X())**2 + (muY2.getVal()- mcPosUV.Y())**2 ) fiteigenproblem=ROOT.TMatrixDEigen(ROOT.TMatrixD(fitCovMvg)) fiteigenwert=eigenproblem.GetEigenValues() fiteigenvector=eigenproblem.GetEigenVectors() tfiteigenvec1=ROOT.TVector3(fiteigenvector[0][0],fiteigenvector[1][0],0) tfiteigenvec2=ROOT.TVector3(fiteigenvector[0][1],fiteigenvector[1][1],0) resUV=ROOT.TVector3((muX2.getVal()- mcPosUV.X()),(muY2.getVal()- mcPosUV.Y()),0) resA=ROOT.TVector3(resUV*tfiteigenvec1,resUV*tfiteigenvec2,0) text=[] text.append('Cluster pos: {0:6.3f} {1:6.3f}'.format(clPosUV.X(),clPosUV.Y())) text.append('Cog pos: {0:6.3f} {1:6.3f}'.format(cogUV.X(),cogUV.Y())) text.append('New Pos: {0:6.3f} {1:6.3f}'.format(muX2.getVal(),(muY2.getVal()))) text.append('Difference%: {0:8.4f} {1:8.4f}'.format( 100*abs(muX2.getVal()-clPosUV.X())/clPosUV.X() ,100*abs(muY2.getVal()-clPosUV.Y())/clPosUV.Y() )) text.append('MC res: fit={0:6.5f} cog={1:6.5f}'.format(fitRes ,cogRes)) text.append('fit-cog/cog={0:8.6f}'.format((fitRes-cogRes)/cogRes)) #text.append('error: U={0:6.3f} V={1:6.3f}'.format(muX2.getError(),muY2.getError())) text.append('error: A1={0:6.3f} A2={1:6.3f}'.format(math.sqrt(fiteigenwert[0][0]), math.sqrt(fiteigenwert[1][1]))) text.append('variance (cov): U={0:6.3f} V={1:6.3f}'.format(cov[0][0],cov[1][1])) #text.append('pullU: {0:6.3f} pullV: {1:6.3f}'.format( (muX2.getVal()- mcPosUV.X())/muX2.getError() , (muY2.getVal()- mcPosUV.Y())/muY2.getError())) text.append('pullA1:{0:6.3f} pullA2: {1:6.3f}'.format(resA.X()/math.sqrt(fiteigenwert[0][0]), resA.Y()/math.sqrt(fiteigenwert[1][1]))) text.append('fittime: {0:6.3f}'.format(runtime)) 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 del varU del varV del muX del muY del muX2 del muY2 del muVec del xVec del varset del mvgResult del mvg del dataset #del datasetProj1 #del datasetProj2 del digiTuple del varAmp del sigX del sigY del rho 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?')