import ROOT, argparse, operator, numpy, math, sys sys.path.append('$PANDAPATH/macro/tpc/FOPI/mberger') from functions import * parser=argparse.ArgumentParser(description='Plot stuff along a track') parser.add_argument('filename',help='filename of the reco file',type=str) parser.add_argument('--rootfilename', help='name of the root file', type=str, default="not") parser.add_argument('--cut', help='Real Data?', type=str, default="not") #parser.add_argument('--doallmode',help='if you want to do more', type=str) parser.add_argument('--ev',help='which event to go to',type=int) parser.add_argument('--zcut',help='zcut',type=int,default=-80) parser.add_argument('--data',help='analyse real data',action='store_true') parser.add_argument('--pdir',help='the picture print directory',type=str,default="not") args=parser.parse_args() ROOT.gROOT.ProcessLine(".x rootlogon.C") ROOT.gROOT.ProcessLine('gSystem->Load("libPhysics")') ROOT.gROOT.ProcessLine('gStyle->SetPalette(1)') ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') ROOT.gROOT.ProcessLine('gSystem->Load("libGeom")') ROOT.gROOT.LoadMacro("stlPYROOT.h+") rfile=ROOT.TFile.Open(args.filename,"read") if rfile.IsOpen()==ROOT.kTRUE: tree = rfile.Get("cbmsim") if tree==None: print "No cbmsim found" if type(Rfile)==ROOT.TFile: Rfile.Close() exit() Slices = [ -48.0, -36.0, -24.0, -12.0, 0.0, 20.0] resbins = 200 #Create Hisogramms--------------------------------------------------------- hdpcz = ROOT.TH2D("hdpcz", "Digis per cm over z", 90, -70, 20, 50, 0, 5) hresz = ROOT.TH2D("hresz", "ResX over z", 90, -70, 20, 2000, -0.5, 0.5) htcoc = ROOT.TH1D("htcoc", "count of the cluster - track", 20,0,20) hacoc = ROOT.TH1D("hacoc", "count of the cluster - Array", 70,0,70) hdcoc = ROOT.TH1D("hdcoc", "difference - cout of the cluster", 10,0,10) hcs2D = ROOT.TH2D("hcs2D", "ClusterSize 2D over z", 90,-70,20,40,0,40) hcs3D = ROOT.TH2D("hcs3D", "ClusterSize 3D over z", 90,-70,20,40,0,40) hdsoc = ROOT.TH2D("hdsoc", "difference ClusterSize", 90,-70,20,20,0,20) hslices = ROOT.TH1D("hslices", "hslices", 90, -70, 20) htheta = ROOT.TH1D("htheta", "htheta", 180, 0, 180) hphi = ROOT.TH1D("hphi", "hphi", 180, 0, 180) hchiz = ROOT.TH2D("hchiz", "Chi^2 over z", 90, -70, 20, 200, 0, 20) hchi = ROOT.TH1D("hchi", "Chi^2", 200, 0, 20 ) HXres1 = ROOT.TH1D("StatsResX Z:"+str(int(Slices[0])), "ResX" ,resbins,-10000, 10000) HXres2 = ROOT.TH1D("StatsResX Z:"+str(int(Slices[1])), "ResX", resbins,-10000, 10000) HXres3 = ROOT.TH1D("StatsResX Z:"+str(int(Slices[2])), "ResX",resbins,-10000, 10000) HXres4 = ROOT.TH1D("StatsResX Z:"+str(int(Slices[3])), "ResX" ,resbins,-10000, 10000) HXres5 = ROOT.TH1D("StatsResX Z:"+str(int(Slices[4])), "ResX" ,resbins,-10000, 10000) HXres6 = ROOT.TH1D("StatsResX Z:"+str(int(Slices[5])), "ResX" ,resbins,-10000, 10000) #hzpos1 = ROOT.TH1D("hzpos1", "z pos", 90, -70, 20) #hzpos2 = ROOT.TH1D("hzpos1", "z pos", 90, -70, 20) #hzpos3 = ROOT.TH1D("hzpos1", "z pos", 90, -70, 20) #hzpos4 = ROOT.TH1D("hzpos1", "z pos", 90, -70, 20) #hzpos5 = ROOT.TH1D("hzpos1", "z pos", 90, -70, 20) #hzpos6 = ROOT.TH1D("hzpos1", "z pos", 90, -70, 20) c1 = ROOT.TCanvas("c1", "Plot 1", 700,1000) c2 = ROOT.TCanvas("c2", "Plot 2", 700,1000) c3 = ROOT.TCanvas("c3", "Plot 3", 700,1000) c4 = ROOT.TCanvas("c4", "Plot 4", 700,1000) c5 = ROOT.TCanvas("c5", "Plot 5", 700,1000) c6 = ROOT.TCanvas("c6", "Plot 6", 700,1000) c7 = ROOT.TCanvas("c7", "Plot 6", 700,1000) c1.Divide(1,2) c2.Divide(1,2) c3.Divide(1,2) c4.Divide(1,2) c5.Divide(1,2) c6.Divide(1,2) c7.Divide(1,2) #---------------------------------------------------------------------- oldpos = ROOT.TVector3(0,0,0) evnum = -1 dist = 0 summe_z = 0 digisize = 0 anzahl = 0 entries = 0 cuts = 0 print "events:", tree.GetEntries() for e in tree: evnum += 1 trackindex = -1 print "tracks:", e.TrackPostFit.GetEntries() if e.TrackPostFit.GetEntries() == 1: entries = 1 if evnum < args.ev: continue unterbrechen = 0 for tr in e.TrackPostFit: trackindex += 1 cand = tr.getCand() tfs = e.TrackFitStat_0.At(trackindex) trackrep=tr.getTrackRep(0) fitflag=trackrep.getStatusFlag() resX=tfs.GetResX() resY=tfs.GetResY() resZ=tfs.GetResZ() resU=tfs.GetResU() resV=tfs.GetResV() theta=tr.getMom().Theta()*ROOT.TMath.RadToDeg() phi = tr.getMom().Phi()*ROOT.TMath.RadToDeg() if args.data: hitIds = cand.getHitIDs(8) if args.cut!= "not": hitIds = cand.getHitIDs(11) else: hitIds = cand.getHitIDs(11) hitindex=0 if resX.size()==0: continue #Cut------------------------------------------ # print "theta:", theta, " ", "phi:", phi, " fitflag:", fitflag htheta.Fill(theta) hphi.Fill(abs(phi)) if args.cut!="not": if theta<40 or theta>140 or fitflag ==1: cuts+=1 print "*************CUT" break #Cut------------------------------------------ dist = 0 laenge=0 cluster = 0 dist1 = 0 oldpos = ROOT.TVector3(0,0,0) #--------------------------------------------------------------------------- rep = tr.getTrackRep(0) chi2 = rep.getChiSqu() print "Chi2:", chi2 NDF = rep.getNDF() print "NDF:", NDF if NDF == 0: quod = 5000 else: quod = chi2/NDF print "QUOD:", quod #---------------------------------------------------------------------------- cluster_posz = [] clsumamp = 0 clsum = 0 for ihit in hitIds: cluster +=1 hitindex +=1 if hitindex>=resX.size(): continue hit=e.TpcSPHit.At(ihit) clid=hit.getClusterID() cl=e.TpcCluster.At(clid) pos=cl.pos() size = cl.size() size2D = cl.get2DSize() cluster_posz.append(cl.pos().Z()) #print cl.pos().Z() #print cluster_posz #-------Cut---------------------- if args.cut!="not": amp = cl.amp() clsumamp += amp clsum += 1 #--------------------------------- if pos.Z() 4: dist +=0 else: dist+=(pos-oldpos).Mag() summe_z += cl.pos().Z() digisize += size anzahl = hitindex oldpos=pos diff_clustersize = math.fabs(size-size2D) # cluster_posz.append(cl.pos().Z()) hcs2D.Fill(cl.pos().Z(), size2D) hcs3D.Fill(cl.pos().Z(), size) hdsoc.Fill(cl.pos().Z(), diff_clustersize) if args.cut!="not": if clsum != 0: clmean = clsumamp/clsum else: clmean = -1 if args.cut!="not" and clsum < 400: break # else: # continue # print "Track Cluster", cluster # print "Array Cluster", cluster_array #1111 if laenge = dist zpos = summe_z summe_z = 0 gesdigi = digisize digisize = 0 if anzahl != 0: mw_pos = zpos/anzahl else: mw_pos = 100 if laenge != 0: cpc = anzahl/laenge else: cpc = -1 if laenge != 0: dpc = gesdigi/laenge else: dpc = 0 #--------------------------------------------------------------------- if args.cut == "not": for resCol in e.MC_Cluster_Residuals: for iRes in range(0, resCol.getSize()): res = resCol.getResidual(iRes); z = res.getRefPos().Z()# - zGem; cl = e.TpcCluster.At(res.getHitIndex()) clAmp = cl.amp() clSize2D = cl.get2DSize() clSize3D = cl.size(); resX = res.getResXYZ().X() errX = cl.sig().X(); hresz.Fill(z, resX) hslices.Fill(z) if z < -48.0: HXres1.Fill(resX*10000) if -48.0 < z < -36.0: HXres2.Fill(resX*10000) if -36.0 < z < -24.0: HXres3.Fill(resX*10000) if -24.0 < z < -12.0: HXres4.Fill(resX*10000) if -12.0 < z < 0.0: HXres5.Fill(resX*10000) if 0.0 < z < 20.0 : HXres6.Fill(resX*10000) # for i in resX: # print i # print "****" # for j in cluster_posz: # print # if args.cut != "not": # for i in range(0, len(resX)): # hresz.Fill(cluster_posz[i], resX[i]) # if cluster_posz[i] < -48.0: # HXres1.Fill(resX[i]*10000) # if -48.0 < cluster_posz[i] < -36.0: # HXres2.Fill(resX[i]*10000) # if -36.0 < cluster_posz[i] < -24.0: # HXres3.Fill(resX[i]*10000) #if -24.0 < cluster_posz[i] < -12.0: # HXres4.Fill(resX[i]*10000) #if -12.0 < cluster_posz[i] < 0.0: # HXres5.Fill(resX[i]*10000) #if 0.0 < cluster_posz[i]< 20.0 : # HXres6.Fill(resX[i]*10000) cluster_array = 0 for cl in e.TpcCluster: cluster_array += 1 # print "Track Cluster", cluster # print "Array Cluster", cluster_array difference_count_cluster = math.fabs(cluster - cluster_array) #Fill----------------------------------------------------------------- hdpcz.Fill(mw_pos, dpc) htcoc.Fill(cluster) hacoc.Fill(cluster_array) hdcoc.Fill(difference_count_cluster) hchiz.Fill(mw_pos, quod) hchi.Fill(quod) #Mean----------------------------------------------------------------- #lable Histogramms------------------------------------------------- hresz.GetXaxis().SetTitle("z-pos [cm]") hresz.GetYaxis().SetTitle("ResX [cm]") htcoc.GetXaxis().SetTitle("count of the cluster") htcoc.GetYaxis().SetTitle("count") hacoc.GetXaxis().SetTitle("count of the cluster") hacoc.GetYaxis().SetTitle("count") hdcoc.GetXaxis().SetTitle("difference count of the cluster") hdcoc.GetYaxis().SetTitle("count") hdcoc.GetXaxis().SetTitle("residuals (x-pos)") hdcoc.GetYaxis().SetTitle("count") hchi.GetXaxis().SetTitle("chi^2/ndf") hchi.GetYaxis().SetTitle("count") hchiz.GetXaxis().SetTitle("z-pos [cm]") hchiz.GetYaxis().SetTitle("chi^2/ndf") #Output-------------------------------------------------------- print "*******Events:", evnum print "*******Cuts:", cuts c1.cd(1) hdpcz.Draw("colz") c1.cd(2) hdpcz.FitSlicesY(0,10,-1,0,"QNRL") hdpcz2=ROOT.gDirectory.Get("hdpcz_1") hdpcz2.Draw() c1.Update() c2.cd(1) hresz.Draw("colz") c2.cd(2) hresz.FitSlicesY(0,10,-1,0,"QNRL") hresz2=ROOT.gDirectory.Get("hresz_2") #print "histname:",hresz2.GetName() hresz2.Draw() c2.Update() c3.cd(1) htcoc.Draw() c3.cd(2) hacoc.Draw() c3.Update() c4.cd(1) hdcoc.Draw() c4.cd(2) hdsoc.Draw("colz") c4.Update() c5.cd(1) hcs2D.Draw("colz") c5.cd(2) hcs3D.Draw("colz") c5.Update() c6.cd(1) hchi.Draw() c6.cd(2) hchiz.Draw("colz") c6.Update() c7.cd(1) htheta.Draw() c7.cd(2) hphi.Draw() c7.Update() liste = [HXres1,HXres2,HXres3,HXres4,HXres5,HXres6, hslices,hdpcz,hresz,htcoc,hresz2,hacoc,hdcoc,hcs2D,hcs3D,hdsoc,hdpcz2] #Save-Pictures --------------------------------------------------------- if args.rootfilename!="not": residualPlot = ROOT.TFile("pics/root/res_plot_"+args.rootfilename+".root", "create") for i in liste: i.Write() residualPlot.Close() else: residualPlot = ROOT.TFile("pics/root/residualPlot.root", "create") for i in liste: i. Write() residualPlot.Close() if args.pdir!="not": c1.SaveAs(args.pdir+"Digis per cm over z.eps") c2.SaveAs(args.pdir+"resX over z.eps") c3.SaveAs(args.pdir+"count of cluster.eps") c4.SaveAs(args.pdir+"difference.eps") c5.SaveAs(args.pdir+"ClusterSize.eps") c6.SaveAs(args.pdir+"chi2_ndf.eps") c7.SaveAs(args.pdir+"phi_theta.eps") #------------------------------------------------------------ u = raw_input("Exit with q") if u == 'q': exit()