import ROOT, argparse, operator, numpy, math, sys 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('--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() #Create Hisogramms--------------------------------------------------------- hcpt = ROOT.TH1D("hcpt", "Cluster per track",45, 0, 45 ) hcpc = ROOT.TH1D("hcpc", "Cluster per cm", 200, 0,2) hldt = ROOT.TH1D("hldt", "Length of the tracks", 40, 0, 40) hnof = ROOT.TH1D("hnof", "number of tracks", 2,0,2) hdbc = ROOT.TH1D("hdbc", "distance between cluster", 150, 0, 15) #hdpcz = ROOT.TH2D("hdpcz", "Digis per cm over z", 90, -70, 20, 300, 0, 30) #hresz = ROOT.TH2D("hresz", "ResX over z", 90, -70, 20, 2000, -0.5, 0.5) hdpc1 = ROOT.TH1D("hdpc", "Digis per Cluster, 1st sector", 20, 0, 20) hdpc2 = ROOT.TH1D("hdpc", "Digis per Cluster, 2nd sector", 25, 0, 25) hdpc3 = ROOT.TH1D("hdpc", "Digis per Cluster, 3rd sector", 30, 0, 30) hdpc4 = ROOT.TH1D("hdpc", "Digis per Cluster, 4th sector", 35, 0, 35) hdpc5 = ROOT.TH1D("hdpc", "Digis per Cluster, 5th sector", 40, 0, 40) hdpc6 = ROOT.TH1D("hdpc", "Digis per Cluster, 6th sector", 45, 0, 45) hcpc1 = ROOT.TH1D("hcpc1", "Cluster per cm, 1st sector", 200, 0, 2) hcpc2 = ROOT.TH1D("hcpc2", "Cluster per cm, 2nd sector", 200, 0, 2) hcpc3 = ROOT.TH1D("hcpc3", "Cluster per cm, 3rd sector", 200, 0, 2) hcpc4 = ROOT.TH1D("hcpc4", "Cluster per cm, 4th sector", 200, 0, 2) hcpc5 = ROOT.TH1D("hcpc5", "Cluster per cm, 5th sector", 200, 0, 2) hcpc6 = ROOT.TH1D("hcpc6", "Cluster per cm, 6th sector", 200, 0, 2) hdpcm1 = ROOT.TH1D("hdpc", "Digis per cm, 1st sector", 30, 0, 30) hdpcm2 = ROOT.TH1D("hdpc", "Digis per cm, 2nd sector", 30, 0, 30) hdpcm3 = ROOT.TH1D("hdpc", "Digis per cm, 3rd sector", 30, 0, 30) hdpcm4 = ROOT.TH1D("hdpc", "Digis per cm, 4th sector", 30, 0, 30) hdpcm5 = ROOT.TH1D("hdpc", "Digis per cm, 5th sector", 30, 0, 30) hdpcm6 = ROOT.TH1D("hdpc", "Digis per cm, 6th sector", 30, 0, 30) c1 = ROOT.TCanvas("c1", "Cluster", 700,1000) c2 = ROOT.TCanvas("c2", "Cluster", 700,1000) c3 = ROOT.TCanvas("c3", "Cluster", 1000,700) c4 = ROOT.TCanvas("c4", "Cluster", 1000,700) c5 = ROOT.TCanvas("c5", "Cluster", 1000,700) c6 = ROOT.TCanvas("c6", "Cluster", 700,1000) c7 = ROOT.TCanvas("c7", "Cluster", 700,1000) c8 = ROOT.TCanvas("c8", "Cluster", 700,1000) c1.Divide(1,2) c2.Divide(1,2) c3.Divide(3,2) c4.Divide(3,2) c5.Divide(3,2) c6.Divide(1,2) c7.Divide(1,2) c8.Divide(1,2) #------------------------------------------------------------------------- #------------------------------------------------------------------------- oldpos = ROOT.TVector3(0,0,0) evnum = -1 dist = 0 summe_z = 0 digisize = 0 anzahl = 0 entries = 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 # hnof.Fill(entries) if evnum < args.ev: continue 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------------------------------------------ if args.cut!="not": if theta<40 or theta>140 or fitflag==1: print "*************CUT" break #Cut------------------------------------------ dist = 0 laenge=0 cluster = 0 dist1 = 0 oldpos = ROOT.TVector3(0,0,0) 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() # print "Y-Pos:", pos.Y() #Cut------------------------------------------- if args.cut!="not": amp = cl.amp() clsumamp += amp clsum += 1 #Cut--------------------------------------------- if pos.Z() 4: dist +=0 else: dist+=(pos-oldpos).Mag() print dist summe_z += cl.pos().Z() digisize += size anzahl = hitindex oldpos=pos if args.cut!="not": if clsum != 0: clmean = clsumamp/clsum else: clmean = -1 if args.cut!="not" and clsum < 400: break #Digis per Cluster - Fill Histogramm if pos.Z() <= -49.85: hdpc1.Fill(size) elif -49.85 < pos.Z() <= -37.7: hdpc2.Fill(size) elif -37.7 < pos.Z() <= -25.55: hdpc3.Fill(size) elif -25.5 < pos.Z() <= -13.4: hdpc4.Fill(size) elif -13.4 < pos.Z() <= -1.25: hdpc5.Fill(size) elif -1.25 < pos.Z() <= 10.9: hdpc6.Fill(size) laenge = dist zpos = summe_z summe_z = 0 gesdigi = digisize digisize = 0 #print "******Event:", evnum #print "Laenge:", laenge #print "Anzahl der Cluster:", anzahl #if laenge > 50: #print "*********zu lang:", evnum 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 # print "mw:", mw_pos # print "cpc:", cpc # print "dcp:", dpc #Fill Histogramms---------------------------------------------- hcpt.Fill(anzahl) hcpc.Fill(cpc) hldt.Fill(laenge) # hnof.Fill(entries) if mw_pos <= -49.85: hcpc1.Fill(cpc) hdpcm1.Fill(dpc) elif -49.85 < mw_pos <= -37.7: hcpc2.Fill(cpc) hdpcm2.Fill(dpc) elif -37.7 < mw_pos <= -25.55: hcpc3.Fill(cpc) hdpcm3.Fill(dpc) elif -25.5 < mw_pos <= -13.4: hcpc4.Fill(cpc) hdpcm4.Fill(dpc) elif -13.4 < mw_pos <= -1.25: hcpc5.Fill(cpc) hdpcm5.Fill(dpc) elif -1.25 < mw_pos <= 10.9: hcpc6.Fill(cpc) hdpcm6.Fill(dpc) if e.TrackPostFit.GetEntries() == 1: entries = 1 hnof.Fill(entries) # hdpcz.Fill(mw_pos, dpc) #----------------------------------------------------------------- #lable Histogramms------------------------------------------------- hcpt.GetXaxis().SetTitle("count of cluster per track") hcpt.GetYaxis().SetTitle("count") hcpc.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc.GetYaxis().SetTitle("count") hldt.GetXaxis().SetTitle("length of the track [cm]") hldt.GetYaxis().SetTitle("count") hnof.GetYaxis().SetTitle("count") hcpc1.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc1.GetYaxis().SetTitle("count") hcpc2.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc2.GetYaxis().SetTitle("count") hcpc3.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc3.GetYaxis().SetTitle("count") hcpc4.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc4.GetYaxis().SetTitle("count") hcpc5.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc5.GetYaxis().SetTitle("count") hcpc6.GetXaxis().SetTitle("cluster per cm [1/cm]") hcpc6.GetYaxis().SetTitle("count") hdpc1.GetXaxis().SetTitle("digis per cluster") hdpc1.GetYaxis().SetTitle("count") hdpc2.GetXaxis().SetTitle("digis per cluster") hdpc2.GetYaxis().SetTitle("count") hdpc3.GetXaxis().SetTitle("digis per cluster") hdpc3.GetYaxis().SetTitle("count") hdpc4.GetXaxis().SetTitle("digis per cluster") hdpc4.GetYaxis().SetTitle("count") hdpc5.GetXaxis().SetTitle("digis per cluster") hdpc5.GetYaxis().SetTitle("count") hdpc6.GetXaxis().SetTitle("digis per cluster") hdpc6.GetYaxis().SetTitle("count") hdpcm1.GetXaxis().SetTitle("digis per cm [1/cm]") hdpcm1.GetYaxis().SetTitle("count") hdpcm2.GetXaxis().SetTitle("digis per cm [1/cm]") hdpcm2.GetYaxis().SetTitle("count") hdpcm3.GetXaxis().SetTitle("digis per cm [1/cm]") hdpcm3.GetYaxis().SetTitle("count") hdpcm4.GetXaxis().SetTitle("digis per cm [1/cm]") hdpcm4.GetYaxis().SetTitle("count") hdpcm5.GetXaxis().SetTitle("digis per cm [1/cm]") hdpcm5.GetYaxis().SetTitle("count") hdpcm6.GetXaxis().SetTitle("digis per cm [1/cm]") hdpcm6.GetYaxis().SetTitle("count") hdbc.GetXaxis().SetTitle("cm") hdbc.GetYaxis().SetTitle("count") #hdpcz.GetXaxis().SetTitle("cm") #hdpcz.GetYaxis().SetTitle("Digis/cm") #hresz.GetXaxis().SetTitle("z-pos [cm]") #hresz.GetYaxis().SetTitle("ResX [cm]") #Output-------------------------------------------------------- print "*******Events:", evnum #Fit with gauss------------------------------------------------- #gausfit1 = ROOT.TF1("gausfit1", "gaus", 0., 20.); #gausfit2 = ROOT.TF1("gausfit2", "gaus", 0., 25.); #gausfit3 = ROOT.TF1("gausfit3", "gaus", 0., 30.); #gausfit4 = ROOT.TF1("gausfit4", "gaus", 0., 35.); #gausfit5 = ROOT.TF1("gausfit5", "gaus", 0., 40.); #gausfit6 = ROOT.TF1("gausfit6", "gaus", 0., 45.); #hdpc1.Fit("gausfit1", "N+") #hdpc2.Fit("gausfit2", "N+") #hdpc3.Fit("gausfit3", "N+") #hdpc4.Fit("gausfit4", "N+") #hdpc5.Fit("gausfit5", "N+") #hdpc6.Fit("gausfit6", "N+") #"[0]*exp(-0.5*((x-[1])/[2])^2)" gaus1 = ROOT.TF1("gaus1", "gaus", 0., 30.); gaus2 = ROOT.TF1("gaus2", "gaus", 0., 30.); gaus3 = ROOT.TF1("gaus3", "gaus", 0., 30.); gaus4 = ROOT.TF1("gaus4", "gaus", 0., 30.); gaus5 = ROOT.TF1("gaus5", "gaus", 0., 30.); gaus6 = ROOT.TF1("gaus6", "gaus", 0., 30.); hdpcm1.Fit("gaus1", "N+") hdpcm2.Fit("gaus2", "N+") hdpcm3.Fit("gaus3", "N+") hdpcm4.Fit("gaus4", "N+") hdpcm5.Fit("gaus5", "N+") hdpcm6.Fit("gaus6", "N+") #Draw Histogramms----------------------------------------------- c1.cd(1) hcpt.Draw() c1.cd(2) hcpc.Draw() c1.Update() c2.cd(1) hldt.Draw() c2.cd(2) hnof.Draw() c2.Update() c3.cd(1) hcpc1.Draw() c3.cd(2) hcpc2.Draw() c3.cd(3) hcpc3.Draw() c3.cd(4) hcpc4.Draw() c3.cd(5) hcpc5.Draw() c3.cd(6) hcpc6.Draw() c3.Update() c4.cd(1) hdpc1.Draw() #gausfit1.Draw("SAME") c4.cd(2) hdpc2.Draw() #gausfit2.Draw("SAME") c4.cd(3) hdpc3.Draw() #gausfit3.Draw("SAME") c4.cd(4) hdpc4.Draw() #gausfit4.Draw("SAME") c4.cd(5) hdpc5.Draw() #gausfit5.Draw("SAME") c4.cd(6) hdpc6.Draw() #gausfit6.Draw("SAME") c4.Update() c5.cd(1) hdpcm1.Draw() gaus1.Draw("SAME") c5.cd(2) hdpcm2.Draw() gaus2.Draw("SAME") c5.cd(3) hdpcm3.Draw() gaus3.Draw("SAME") c5.cd(4) hdpcm4.Draw() gaus4.Draw("SAME") c5.cd(5) hdpcm5.Draw() gaus5.Draw("SAME") c5.cd(6) hdpcm6.Draw() gaus6.Draw("SAME") c5.Update() #c7.cd(1) #hdpcz.Draw("colz") #c7.cd(2) #hdpcz.FitSlicesY(0,0,-1,0,"QNRL") #ROOT.gDirectory.Get("hdpcz_1").Draw() #c7.Update() #c8.cd(1) #hdbc.Draw() #c8.cd(2) #hresz.Draw("colz") #c8.Update() #c8.cd(1) #hresz.Draw("colz") #c8.cd(2) #hresz.FitSlicesY(0,0,-1,0,"QNRL") #ROOT.gDirectory.Get("hresz_2").Draw() #c8.Update() #Graph mean---------------------------------------------------- #Get means clmean1 = hcpc1.GetMean() clmean2 = hcpc2.GetMean() clmean3 = hcpc3.GetMean() clmean4 = hcpc4.GetMean() clmean5 = hcpc5.GetMean() clmean6 = hcpc6.GetMean() digimean1 = hdpcm1.GetMean() digimean2 = hdpcm2.GetMean() digimean3 = hdpcm3.GetMean() digimean4 = hdpcm4.GetMean() digimean5 = hdpcm5.GetMean() digimean6 = hdpcm6.GetMean() #Create Graphs clmeangraph = ROOT.TGraph() digimeangraph = ROOT.TGraph() #Fill Graphs clmeangraph.SetPoint(0, -55.93, clmean1) clmeangraph.SetPoint(1, -43.78, clmean2) clmeangraph.SetPoint(2, -31.63, clmean3) clmeangraph.SetPoint(3, -19.48, clmean4) clmeangraph.SetPoint(4, -7.33, clmean5) clmeangraph.SetPoint(5, 4.83, clmean6) digimeangraph.SetPoint(0, -55.93, digimean1) digimeangraph.SetPoint(1, -43.78, digimean2) digimeangraph.SetPoint(2, -31.63, digimean3) digimeangraph.SetPoint(3, -19.48, digimean4) digimeangraph.SetPoint(4, -7.33, digimean5) digimeangraph.SetPoint(5, 4.83, digimean6) #Stile of the Graphs clmeangraph.SetTitle("cluster per cm - mean") digimeangraph.SetTitle("digi per cm - mean") clmeangraph.GetXaxis().SetTitle("z-pos [cm]") clmeangraph.GetYaxis().SetTitle("mean [1/cm]") clmeangraph.SetMarkerStyle(20) clmeangraph.SetMarkerSize(1) clmeangraph.SetMarkerColor(2) digimeangraph.GetXaxis().SetTitle("z-pos [cm]") digimeangraph.GetYaxis().SetTitle("mean [1/cm]") digimeangraph.SetMarkerStyle(20) digimeangraph.SetMarkerSize(1) digimeangraph.SetMarkerColor(2) Gerade = ROOT.TF1("Gerade", "[0]*x + [1]", -60., 10.); digimeangraph.Fit(Gerade,"N") steigung=Gerade.GetParameter(0) yAbsch=Gerade.GetParameter(1) text=ROOT.TLatex() text.SetTextColor(ROOT.kBlack) text.SetTextSize(0.03) text.SetNDC() y=0.8 ystep=0.05 #text.DrawLatex(0.15,y,str(steigung)) #Draw Graphs c6.cd(1) clmeangraph.Draw("AP") c6.cd(2) digimeangraph.Draw("AP") Gerade.Draw("SAME") text.DrawLatex(0.15,y,"y="+str(steigung)+"x"+"+"+str(yAbsch)) c6.Update() #Save Pictures----------------------------------------------- liste = [hcpt,hcpc,hldt,hnof,hdbc,hdpc1,hdpc2,hdpc3,hdpc4,hdpc5,hdpc6,hcpc1,hcpc2, hcpc3, hcpc4,hcpc5,hcpc6,hdpcm1,hdpcm2,hdpcm3,hdpcm4,hdpcm5,hdpcm6] if args.rootfilename!="not": clusterPlot = ROOT.TFile("pics/root/clu_plot_"+args.rootfilename+".root", "create") for i in liste: i.Write() clusterPlot.Close() else: clusterPlot = ROOT.TFile("pics/root/cluster_plot.root", "create") for i in liste: i.Write() clusterPlot.Close() if args.pdir!="not": c1.SaveAs(args.pdir+"ClusterPerTrack.eps") c2.SaveAs(args.pdir+"TrackLength.eps") c3.SaveAs(args.pdir+"ClusterPerCm.eps") c4.SaveAs(args.pdir+"DigisPerCluster.eps") c5.SaveAs(args.pdir+"DigisPerCm.eps") c6.SaveAs(args.pdir+"meangraphs.eps") c7.SaveAs(args.pdir+"Digis per cm over z.eps") c8.SaveAs(args.pdir+"Distance between cluster.eps") #------------------------------------------------------------ u = raw_input("Exit with q") if u == 'q': exit() #Auskommentiert-------------------------------------------------- # array = [] # for ihit in hitIds: # if hitindex>=resX.size(): # continue # hit=e.TpcSPHit.At(ihit) # clid=hit.getClusterID() # cl=e.TpcCluster.At(clid) # pos=cl.pos() # array.append([pos.X(), pos.Y(), pos.Z()]) # array.sort(key=lambda a: a[1]) # # for i in array: # dist = sqrt((i[0]-i+1[0])) # # if cluster == 2: # dist1 = dist # print "Dist 1:", dist1 # print "Dist:", dist