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 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.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) 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 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 for idigi in range(len(digisOnPlane)): digipos=digisOnPlane[idigi] digiamp=digisOnPlaneAmp[idigi] #digiTuple.Fill(digipos.X(),digipos.Y(),digiamp*0.05) graph1.SetPoint(idigi,digipos.X(),digiamp) graph2.SetPoint(idigi,digipos.Y(),digiamp) graph2D.SetPoint(idigi,digipos.X(),digipos.Y(),digiamp) 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() teigenvector=ROOT.TVector3(eigenvector[0][0],eigenvector[1][0],0) teigenvector2=ROOT.TVector3(eigenvector[0][1],eigenvector[1][1],0) if cov[0][0]==0 and cov[1][1]==0: skipTrack=True print 'skipping track due to cov!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' continue #reclaculate digipositions in coordinate system of cov for idigi in range(len(digisOnPlane)): digipos=digisOnPlane[idigi] digiamp=digisOnPlaneAmp[idigi] eigenDigiPos=ROOT.TVector3(0,0,0) eigenDigiPos.SetX(digipos*teigenvector) eigenDigiPos.SetY(digipos*teigenvector2) digiTuple.Fill(eigenDigiPos.X(),eigenDigiPos.Y(),digiamp*0.05) x=cogUV*teigenvector y=cogUV*teigenvector2 cogUV.SetX(x) cogUV.SetY(y) x=mcPosUV*teigenvector y=mcPosUV*teigenvector2 mcPosUV.SetX(x) mcPosUV.SetY(y) 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) ucogmin=cogUV.X()-0.5*math.sqrt(cov[0][0]) ucogmax=cogUV.X()+0.5*math.sqrt(cov[0][0]) vcogmin=cogUV.Y()-0.5*math.sqrt(cov[1][1]) vcogmax=cogUV.Y()+0.5*math.sqrt(cov[1][1]) varU=ROOT.RooRealVar('u','',umins,umaxs) varV=ROOT.RooRealVar('v','',vmins,vmaxs) varAmp=ROOT.RooRealVar('amp','',0,1.5*digiTuple.GetMaximum('amp')) varset.add(varU) varset.add(varV) varset.add(varAmp) dataset=ROOT.RooDataSet('dataSet','digisOnPlane',digiTuple,varset,'') muX=ROOT.RooRealVar('mux','mean of gauss in U',0,ucogmin,ucogmax) muY=ROOT.RooRealVar('muy','mean of gauss in V',0,vcogmin,vcogmax) if not args.batch: print 'drawing' hd=varU.createHistogram('digis on plane',varV,'') dataset.fillHistogram(hd,ROOT.RooArgList(varU,varV)) canv.cd(1) hd.SetStats(0) hd.Draw('colz') 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) pcov.SetFillStyle(4000) pcov.Draw("same") #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') pnewpos=ROOT.TPolyMarker() pnewpos.SetPoint(0,muX.getVal(),muY.getVal()) pnewpos.SetMarkerStyle(21) pnewpos.SetMarkerColor(ROOT.kPink) pnewpos.Draw('same') print 'drawing gaus' canv.cd(2) h3D=ROOT.RooAbsData.createHistogram(dataset,'h3D',varU,ROOT.RooFit.YVar(varV),ROOT.RooFit.ZVar(varAmp)) h3D.Draw() #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=gausvar2damp.createHistogram('varU,varV',100,100,100) #hgmult.Draw('surf') 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( (muX.getVal()- mcPosUV.X())**2 + (muY.getVal()- mcPosUV.Y())**2 ) 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(muX.getVal(),(muY.getVal()))) text.append('Difference%: {0:8.4f} {1:8.4f}'.format( 100*abs(muX.getVal()-clPosUV.X())/clPosUV.X() ,100*abs(muY.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} correlation={2:6.3f}'.format(muX.getError(),muY.getError(),1)) #text.append('pullU: {0:6.3f} pullV: {1:6.3f}'.format( (muX.getVal()- mcPosUV.X())/muX.getError() , (muY.getVal()- mcPosUV.Y())/muY.getError())) #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 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?')