import sys,os,time,math,array pandapath=os.environ.get('PANDAPATH') sys.path.append(pandapath+'/macro/tpc/FOPI/mberger') sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1') import glob,math,argparse,ROOT from functions import set_palette, openTree 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)): weights.append(1) return weights 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 def getDiffusionWidth(zpos,direction=1): retval=0 if direction==1 or direction==2: retval=0.0197879*math.sqrt(zpos) if direction==3: retval=0.0259143*math.sqrt(zpos) return (retval**2) def getParamCov(zpos,amp,planeMatrix): cov=ROOT.TMatrixDSym(3) cov[0][0]=20*getDiffusionWidth(zpos,1) cov[1][1]=20*getDiffusionWidth(zpos,2) cov[2][2]=20*getDiffusionWidth(zpos,3) print 'generating param cov' cov.Print() print 'making covOnPlane zpos=',zpos planeMatrix.Print() print 'cov' cov.Print() print 'cov inverted' cov.Invert() cov.Print() covOnPlane=cov.SimilarityT(planeMatrix) print 'cov On plane' covOnPlane.Print() covOnPlane.Invert() print 'cov on plane inverted' covOnPlane.Print() cov.Print() return covOnPlane parser=argparse.ArgumentParser(description="plots digisposition projected on detector plane") parser.add_argument('recofile',help='the reco file',type=str,default='./') parser.add_argument('--goto',help='go to event',type=int,default=-1) parser.add_argument('--gotoCl',help='go to Cluster',type=int,default=-1) parser.add_argument('--events',help='number of events to run',type=int,default=-1) 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('--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('--batch',help='batch mode',action='store_true') parser.add_argument('--incrRange',help='increase the range of the fit',action='store_true') parser.add_argument('--extraPlot',help='plot some extra stuff',action='store_true') 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) tree=openTree(args.recofile) rfile=ROOT.TFile(args.outfile,'recreate') 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.pfile!='None': canv.SaveAs(args.pfile+'[') if os.path.isdir(args.recofile): tree=openTree(args.recofile) else: RFile=ROOT.TFile(args.recofile,'read') tree=RFile.Get('cbmsim') fittime=ROOT.TH1D('fittime','fittime multi gauss',1000,0,2) ev=-1 u='' 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.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() W=tr.GetMidPlaneW() if len(U)==0: U=tr.GetPlaneU() V=tr.GetPlaneV() O=tr.GetPlaneO() W=tr.GetPlaneW() if len(W)==0: W=[] for i in range(len(U)): W.append(U[i].Cross(V[i])) trackPoints=tr.GetProjectionPoints() ids=tr.GetHitIDs() first=True print tr.nhits(),'cluster in track' for c in range(tr.nhits()): if args.gotoCl>c and args.gotoCl!=-1: continue 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) if cluster.size2D()<3: continue print len(O),len(U) clPos=ROOT.TVector3(tr.GetHitPosition(c)-O[c]) #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() #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) graph2D=ROOT.TGraph2D() graph2D.SetTitle('Digis On Plane') digiTuple=ROOT.TNtuple('UVdigipos','','u:v:amp') digisOnPlane=tr.GetDigisOnPlane(c) digisOnPlaneAmp=tr.GetDigisOnPlaneAmp(c) digisMCPos=tr.GetDigiMCPos(c) digisPos=tr.GetDigiPos(c) refmom=tr.GetProjectionPointsMom(c) digisLength=tr.GetDigiLen(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 sumlength=0 for idigi in range(len(digisOnPlane)): weight=weights[idigi] digipos=digisOnPlane[idigi] digiamp=abs(digisOnPlaneAmp[idigi]) weighted_amp_sum+=digiamp*weight cogUV.SetX(cogUV.X()+(weight*digiamp*digipos.X())) cogUV.SetY(cogUV.Y()+(weight*digiamp*digipos.Y())) digiAmpSum+=digiamp maxdigiamp=max(maxdigiamp,abs(digiamp)) sumlength+=digisLength[idigi] #cogUV*=1./digiAmpSum cogUV*=1./weighted_amp_sum pdigiMC=ROOT.TPolyMarker() for idigi in range(len(digisOnPlane)): digiposUV=digisOnPlane[idigi] digiamp=abs(digisOnPlaneAmp[idigi]) digiweight=weights[idigi] digiamp*=digiweight tmp1=ROOT.TVector3(U[c]) tmp2=ROOT.TVector3(V[c]) tmp1*=digiposUV.X() tmp2*=digiposUV.Y() digipos3D=ROOT.TVector3(tmp1+tmp2+O[c]) tmp3=ROOT.TVector3(tr.GetHitPosition(c)) dirToPad=digipos3D-tmp3 graph2D.SetPoint(idigi,digiposUV.X(),digiposUV.Y(),digiamp) 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=digiposUV.X()-cogUV.X() mresY=digiposUV.Y()-cogUV.Y() ''' if mresX<0: mresX-=math.sqrt( (digisLength[idigi]/(2*sumlength)) *abs(zjitterOnPlane .X())) else: mresX+=math.sqrt( (digisLength[idigi]/(2*sumlength)) *abs(zjitterOnPlane .X())) if mresY<0: mresY-=math.sqrt( (digisLength[idigi]/(2*sumlength)) *abs(zjitterOnPlane .Y())) else: mresY+=math.sqrt( (digisLength[idigi]/(2*sumlength)) *abs(zjitterOnPlane .Y())) ''' cm[0][0]=mresX*math.sqrt(spreadweight1) cm[1][0]=mresY*math.sqrt(spreadweight2) cm_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,cm) acov=ROOT.TMatrixD(cm,ROOT.TMatrixD.kMult,cm_t) acov*=digiamp*digiweight digiTuple.Fill(digiposUV.X(),digiposUV.Y(),3771*digiamp/weighted_amp_sum)#*digiweight) cov+=acov cov*=1./(weighted_amp_sum) cov=getParamCov(O[c].Z(), weighted_amp_sum, UVplane) #cov*=1./len(digisOnPlane) cov=ROOT.TMatrixDSym(2,cov.GetMatrixArray()) cov.Print() #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=cogUV.x()-3*math.sqrt(cov[0][0]) umaxs=cogUV.X()+3*math.sqrt(cov[0][0]) vmins=cogUV.Y()-3*math.sqrt(cov[1][1]) vmaxs=cogUV.Y()+3*math.sqrt(cov[1][1]) ''' umins-=digiErr.X()/2. umaxs+=digiErr.X()/2. vmins-=digiErr.Y()/2. vmaxs+=digiErr.Y()/2. ''' 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) varV=ROOT.RooRealVar('v','',vmins,vmaxs) varAmp=ROOT.RooRealVar('amp','',0,digiTuple.GetMaximum('amp')+500) varset.add(varU) varset.add(varV) varset.add(varAmp) dataset=ROOT.RooDataSet('dataSet','digisOnPlane',digiTuple,varset,'','amp') #dataset=ROOT.RooDataSet('dataSet','digisOnPlane',digiTuple,varset) 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) mvg=ROOT.RooMultiVarGaussian("mvg","mvg",xVec,muVec,cov) start=time.time() mvgResult=mvg.fitTo(dataset, ROOT.RooFit.Save(), ROOT.RooFit.SumW2Error(False)#, #ROOT.RooFit.Verbose(False), #ROOT.RooFit.PrintLevel(-1),ROOT.RooFit.PrintEvalErrors(0), #ROOT.RooFit.Warnings(False) ) if (mvgResult.status()!=0): print 'RooFit mvg failed {0} {1}'.format(mvgResult.status(), mvgResult.covQual()) #continue fitCovMvg=mvgResult.covarianceMatrix() print 'the fitcov: ' fitCovMvg.Print() print 'the hesse:' hesse=ROOT.TMatrixDSym(2,cov.GetMatrixArray()) hesse*=1./weighted_amp_sum#len(digisOnPlane) hesse.Print() runtime=time.time()-start fittime.Fill(runtime) setGraph(graph2D) 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') hd.GetXaxis().SetTitle("X (cm)") hd.GetYaxis().SetTitle("Y (cm)") hd.GetZaxis().SetTitle("Amp") 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") cogpos=ROOT.TPolyMarker() cogpos.SetPoint(0,cogUV.X(),cogUV.Y()) cogpos.SetMarkerStyle(22) cogpos.SetMarkerSize(2) 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() 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,'') mvg.fillHistogram(hgmult,ROOT.RooArgList(varU,varV)) hgmult.GetXaxis().SetTitle("X (cm)") hgmult.GetYaxis().SetTitle("Y (cm)") hgmult.GetZaxis().SetTitle("Amp") hgmult.Draw('surf') hd.Draw('lego same') canv.cd(3) fittime.Draw() canv.cd(4) canv.cd(4).Clear() 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=fiteigenproblem.GetEigenValues() fiteigenvector=fiteigenproblem.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()))) if(clPosUV.X()!=0 and clPosUV.Y()!=0): 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: 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('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.extraPlot: set_palette('pal',99) exCanv=ROOT.TCanvas("extra","extra",1000,1000) exCanv.SetLeftMargin(0.2) hd.SetTitle("") hd.GetXaxis().SetTitleOffset(2) hd.GetYaxis().SetTitleOffset(2) hd.GetZaxis().SetTitleOffset(2.5) hd.Draw('lego2 0') exCanv.Update() u=raw_input('next') hgmult.SetTitle("") hgmult.GetXaxis().SetTitleOffset(2) hgmult.GetYaxis().SetTitleOffset(2) hgmult.GetZaxis().SetTitleOffset(2.5) hgmult.SetStats(0) hgmult.Draw('surf2') hgmult.Draw('cont1 same') exCanv.Update() u=raw_input('next') 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 muX2 del muY2 del muVec del xVec del varset del mvgResult del mvg del dataset del digiTuple del varAmp 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() while u!='q' and not args.batch: u=raw_input('done?')