import ROOT, glob, math, sys, os, numpy from ROOT import std def sort_inner(inner): #inner is each inner list in the list of lists to be sorted #(here item at index 1 of each inner list is to be sorted) return inner[0] def sort_inner2(inner): return inner[2] dir = "/nfs/hicran/project/panda/SIM/bergerm/dev/oufiles" pdir = "pics/tmp" runs ="reco" nEvents = 10000000000 sim = 0 raw = 0 chii = 1.5 cosmic = 1 sampling = 20. histnameaddon = " data" samthresh = 0 dofit = 0 partrack = 1 nplot = 0 printc = 0 minNumHits = 10 curv = 1 outfname = "" thec = 1 phic = 0 batchmode = 0 sat = 1 momc = 1 floop = 1 slicedummy = 1 debug = 0 infilen = "anafiles/dummyslice.root" driftcut = -100 part=-1 Rfile=ROOT.TFile() # {{{ argument parsing #argument parsng: for iarg in range(len(sys.argv)): arg = sys.argv[iarg] if arg == "-path" : dir = sys.argv[iarg+1]; if arg == "-runs" : runs = sys.argv[iarg+1] if arg == "-nEvts": nEvents = int(sys.argv[iarg+1]) if arg == "-sim" : sim = 1 histnameaddon=" sim" if arg == "-raw" : raw = 1 if arg == "-chi" : chii = float(sys.argv[iarg+1]) if arg == "-b" : ROOT.gROOT.SetBatch(True) batchmode = 1 if arg == "-cosmic" : cosmic = 0 if arg == "-sampling" : sampling = float(sys.argv[iarg+1]) if arg == "-sth" : samthresh = int(sys.argv[iarg+1]) if arg == "-nfit" : dofit = 0 if arg == "-par": partrack = 0 if arg == "-nplo" : nplot = 1 if arg == "-print": printc = 1 if arg == "-pdir" : pdir = sys.argv[iarg+1] printc = 1 if arg == "-nhits" : minNumHits = int(sys.argv[iarg+1]) if arg == "-curv" : curv = 0 if arg == "-out": outfname=sys.argv[iarg+1] if arg == "-phic": phic = 1 if arg == "-momc": momc = 0 if arg == "-thec": thec = 0 if arg == "-sat": sat = 0 if arg == "-slice": floop = 0 nplot = 1 if arg == "-slicef": infilen = sys.argv[iarg+1] if arg == "-debug": debug = 1 if arg == "-zcut": driftcut = float(sys.argv[iarg+1]) thec = 0 if arg == "-part": part = (sys.argv[iarg+1]) if arg == "-h": print "-path\t define path of files to analyze" print "-runs\t define run numbers to analyze" print "-nEvts\t define number of events to analyze" print "-sim\t " print "-raw\t" print "-chi\t" print "-b \t enable batch mode, still prints pictures, but no drawing" print "-cosmic\t " print "-sampling\t set sampling frequency" print "-sth\t set sample threshold" print "-nfit\t" print "-par\t disable cut on parralel tracks" print "-nplo\t dont plot anything" exit(0) print "input file path:\t", dir print "run numbers:\t",runs print "out file name:\t",outfname print "number of events:\t",nEvents print "chi2 cut for ??:\t", chii print "batchmode flag:\t\t",batchmode print "cosmic flag (1 track cut):\t",cosmic print "sampling freq:\t\t",sampling print "fitting flag:\t\t",dofit print "plotting flag:\t\t",nplot print "print flag:\t\t",printc print "print directory:\t",pdir print "min number of hits:\t",minNumHits print "sim flag:\t\t",sim print "digi amp threshold:\t",samthresh print "z-cut:\t\t\t\t",driftcut print "parallel track cut:\t",partrack print "phi cut flag:\t\t",phic print "theta cut flag:\t\t",thec print "momentum cut flag:\t",momc print "curvature cut:\t\t",curv print "saturation cut flag:\t",sat print "debug flag: \t\t",debug nix = raw_input("Happy with this settings?") #finished argument parsing # }}} # {{{ create runlist runList=[] colindex = runs.find(",") if colindex > 0: runs = runs.split(",") for i in range(len(runs)) : dashIndex = runs[i].find("-") if dashIndex > 0 : runs[i] = runs[i].split("-") for j in range(int(runs[i][0]), int(runs[i][1])+1) : runList.append(j) # runList = range(int(runs[0]), int(runs[1])+1) if dashIndex < 0 and len(runs) > 1 : runList.append(int(runs[i])) else : dashIndex = runs.find("-") if dashIndex > 0 : runs = runs.split("-") runList = range(int(runs[0]), int(runs[1])+1) if dashIndex < 0 and len(runs) > 1 : runList = [runs] # }}} print "files to process:", runList ROOT.gROOT.ProcessLine(".x rootlogon.C") #ROOT.gROOT.ProcessLine(".x rootlogon_Bernhard.C") #ROOT.gROOT.SetBatch(True) ROOT.gROOT.ProcessLine('gSystem->Load("libPhysics")') ROOT.gROOT.ProcessLine('gStyle->SetPalette(1)') ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') #ROOT.gROOT.ProcessLine('gROOT->SetStyle("col")') ROOT.gROOT.ProcessLine('gSystem->Load("libGeom")') files= glob.glob(dir + "/*reco*.root") files.sort() pedlist = { "4219":"ped_FOPI___start~21.11.2011_16:48:27~__end~01.01.2038_00:00:00~___Src_", "4206":"ped_FOPI___start~20.11.2011_11:03:50~__end~21.11.2011_16:48:26~___Src_", "4169":"ped_FOPI___start~18.11.2011_17:04:12~__end~19.11.2011_17:59:41~___Src_", "4020":"ped_FOPI___start~18.11.2011_17:04:12~__end~19.11.2011_17:59:41~___Src_", "3979":"ped_FOPI___start~14.11.2011_08:33:53~__end~15.11.2011_15:02:28~___Src_", "3902":"ped_FOPI___start~13.11.2011_18:04:19~__end~14.11.2011_08:33:52~___Src_", "3812":"ped_FOPI___start~12.11.2011_10:12:09~__end~01.01.2038_00:00:00~___Src_", "3651":"ped_FOPI___start~09.11.2011_04:27:50~__end~12.11.2011_10:12:08~___Src_", "2613":"ped_FOPI___start~30.04.2011_16:23:37~__end~02.05.2011_00:15:00~___Src_", "2453":"ped_FOPI___start~21.04.2011_17:49:32~__end~21.04.2011_22:52:13~___Src_", "2463":"ped_FOPI___start~21.04.2011_22:52:14~__end~30.04.2011_16:23:36~___Src_", "2450":"ped_FOPI___start~21.04.2011_17:49:32~__end~21.04.2011_22:52:13~___Src_", "2341":"ped_FOPI___start~19.04.2011_20:17:10~__end~21.04.2011_12:43:36~___Src_", "2339":"ped_FOPI___start~19.04.2011_19:15:19~__end~19.04.2011_20:17:09~___Src_", "2337":"ped_FOPI___start~16.04.2011_19:33:28~__end~19.04.2011_16:50:39~___Src_", "2262":"ped_FOPI___start~13.04.2011_13:50:28~__end~14.04.2011_16:24:18~___Src_", "2230":"ped_FOPI___start~12.04.2011_13:53:04~__end~12.04.2011_21:43:13~___Src_", "2213":"ped_FOPI___start~12.04.2011_00:06:04~__end~12.04.2011_13:53:03~___Src_", "2165":"ped_FOPI___start~09.04.2011_22:54:15~__end~11.04.2011_23:59:00~___Src_" } simpath = "" for env in os.environ.keys() : if env == "PANDAPATH" and len(os.environ[env]) > 0 : simpath = os.environ[env] print "Setting pandaroot path to",simpath # {{{ define histos hSampleTot = ROOT.TH1D("hsampletot","all samples together"+histnameaddon,2024,0.,2024.) hOneDigi = ROOT.TH1D("honedigi","amp vs time of one sample"+histnameaddon,511,0.,511.) hAllDigiSamples = ROOT.TH1D("hAllDigisSamples","amp vs time"+histnameaddon,1000,0.,1000.) hOneDigiTemp = ROOT.TH1D("hOneDigiTemp","just for fitting"+histnameaddon,600,0.,600.) hDigiWidth = ROOT.TH1D("hDigiWidth","sigma of digi"+histnameaddon,100,0.,10.*(1/sampling)*1000) hDigiRMS = ROOT.TH1D("hDigiRMS","RMS of digi"+histnameaddon,100,0.,10.) hDigiMaxS = ROOT.TH1D("hDigiMaxS","Maximal sample amplitude"+histnameaddon,1024,0.,1024.) hClusterMax = ROOT.TH1D("hClusterMax","Maximal sample amplitude in Cluster"+histnameaddon,1024,0.,1024.) hintsampleamp = ROOT.TH1D("hintsampleamp","integrated sample amp"+histnameaddon,2500,0.,10000.) hchisquare = ROOT.TH1D("hchisquare","chi2 of samplefit"+histnameaddon,100,0.,10.) hsamplebigger511 = ROOT.TH1D("hsamplebigger511","digi w sample > 511"+histnameaddon,40,491.,531.) hsigmavsnofsampl = ROOT.TH2D("hsigmavsnofsampl","sigma vs nof sample"+histnameaddon,20,0.,20.,100,0.,10.*(1/sampling)*1000) hsigmavschii = ROOT.TH2D("hsigmavschii","sigma vs chii"+histnameaddon,100,0.,20.,100,0.,10.*(1/sampling)*1000) hdigi = ROOT.TH1D("hdigi","all digis same time"+histnameaddon,20,0.,20.) hndigivsnsamp = ROOT.TH2D("hndigivsnsamp","ndigi vs nsamp"+histnameaddon,500,0.,500.,50,0.,50.) hfullpad = ROOT.TH1D("hfllupad","one full pad",511,0.,511.) hnsamperpad = ROOT.TH1D("hnsamperpad","number of samples per pad",20,0.,20.) hnsamperpadpsa = ROOT.TH1D("hnsamperpadpsa","number of samples per pad psa",20,0.,20.) hnpadsndigi = ROOT.TH1D("hnpadsvsndigi","number of pads - number of digis"+histnameaddon,20,-10.,10.) hdigiperpad = ROOT.TH1D("hdigiperpad","digis per pad",10,0.,10.) hnoiseamp = ROOT.TH1D("hnoiseamp","sample amps of noise",100,0.,100.) hsamperdigi = ROOT.TH1D("hsamperdigi","samples per digi",20,0.,20.) hpadspercluster = ROOT.TH1D("hpadspercluster","pads per cluster 3D",100,0.,100.) hpadspercluster2D = ROOT.TH1D("hpadspercluster2D","pads per cluster 2D",100,0.,100.) hsnr = ROOT.TH1D("hsnr","signal to noise",200,0.,200.) hclusterpertrack = ROOT.TH1D("hclusterpertrack","cluster per track",200,0.,200.) hclusterpercm=ROOT.TH1D("hclusterpercm","Cluster per cm",1000,0.,10.) htheta = ROOT.TH1D("htheta","theta angle of tracks",180,0.,180.) hthetacut = ROOT.TH1D("hthetacut","single track theta",180,0.,180.) hthetacurvcut = ROOT.TH1D("hthetacurvcut","curv cut theta",180,0.,180.) hthetagood = ROOT.TH1D("hthetagood","good theta",180,0.,180.) hphi = ROOT.TH1D("hphi","phi angle of tracks",360,-180.,180.) hphitrackcut = ROOT.TH1D("hphitrackcut","phi angle of tracks",360,-180.,180.) hphicurvcut = ROOT.TH1D("hphicurvcut","phi angle of tracks",360,-180.,180.) hphigood = ROOT.TH1D("hphigood","good tracks",360,-180.,180.) trknumdist = ROOT.TH1D("trknumdist","Track number distribution",100,0,100) trklendist = ROOT.TH1D("trklendist","Track length distribution",150,0,30) hitspertrack=ROOT.TH1D("hitspertrack","Hts per Track",40,10.,50.) hallresi=ROOT.TH1D("hallresi","all residuals",400,-10000.,10000.) hdropout=ROOT.TH1D("hdropout","dropouts",401,0.,401.) hxy=ROOT.TH2D("hxy","xy",200,-16.,16.,200,-16.,16.) hrz=ROOT.TH2D("hrz","xz",350,-60.,20.,100,-16.,16.) hmom=ROOT.TH1D("hmom","mom",500,0,2.) hslicehist = ROOT.TH1D("hslicehist","clusterpos",80,-60.,20.) hslices = ROOT.TH1D("hslices","slicepos",6,0.,6.) hthevslen = ROOT.TH2D("hthevslen","theta vs length",360,0.,180.,800,0.,80.) htracklen = ROOT.TH1D("htracklen","track length",800,0.,80.) htracklengood = ROOT.TH1D("htracklengood","track length",800,0.,80.) hresvsclsize = ROOT.TH2D("hresvsclsize","residual vs clustersize",200,-10000.,10000.,300,0.,30.) hresX=ROOT.TH1D("hresX","residuals in X",400,-10000.,10000.) hresY=ROOT.TH1D("hresY","residuals in Y",400,-10000.,10000.) hresZ=ROOT.TH1D("hresZ","residuals in Z",400,-10000.,10000.) hclsizeradius=ROOT.TH1D("hclsizeradius","clustersizesum vs radius",150,0.,15.) hclsizeradiusccut=ROOT.TH1D("hclsizeradiusccut","clustersizesum vs radius for curvcut",150,0.,15.) hclsizeradius2D=ROOT.TH1D("hclsizeradius2D","clustersizesum vs radius",150,0.,15.) hclsizevsradius=ROOT.TH2D("hclsizevsradius","radius vs clustersize",150,0.,15.,100,0.,100.) hclsizevsradius2D=ROOT.TH2D("hclsizevsradius2D","radius vs clustersize",150,0.,15.,100,0.,100.) c99 = ROOT.TCanvas("bla","bla") c98 = ROOT.TCanvas("bla1","bla1") # }}} testfit = ROOT.TF1("trestfit","gaus") fit = ROOT.TF1("fit","gaus") # {{{loop over files to create filelist: totev = 0 filecount = 0 firstloop = 1 filelist = [] filenumber = [] print "Start loop over files" # start loop over files for file in files : # print "actual file:",file if len(runList) >= 1 : found = 0 for i in runList : if file.find(str(i)) >= 0 : found = 1 filecount += 1 # print "File number:",i print "part=",part if part>0: if file.find("reco"+part+".root")>=0 : filelist.append(file) filenumber.append(i) else: filelist.append(file) filenumber.append(i) # if found == 0 : # print "File not in runlist" # continue # }}} zSlices = [] #for floop in range(0,2) : if floop == 1 : print "slicename", infilen infile=ROOT.TFile(infilen) zSlices = [] hslices = infile.Get("hslices") for k in range(1,7): print hslices.GetBinContent(k) zSlices.append(hslices.GetBinContent(k)) #zSlices = ( -52, -41, -29, -16, -2, 20 ) # cosmic ne #zSlices = ( 20, 30, 40, 49, 57, 70 ) # beam cSlices = (3,6,9,12,15,18,21,24,50) rSlices = (6.6,8.2,9.8,11.4,13,15) # {{{ define dicts samwidths = dict([(i, ROOT.TH1D("digiwidth"+str(zSlices[i]),"Digiwidth"+histnameaddon,100,0.,10.*(1/sampling)*1000)) for i in range(len(zSlices))]) samwidthxamp = dict([(i, ROOT.TH1D("digiwidthxamp"+str(zSlices[i]),"Digi width x amp"+histnameaddon,1000,0.,1000.)) for i in range(len(zSlices))]) digissliced2d = dict([(i, ROOT.TH2D("digis2d"+str(zSlices[i]),"Digi shape"+histnameaddon,20,0.,20.,2000,0.,2000.)) for i in range(len(zSlices))]) digissliced = dict([(i, ROOT.TH2D("digis"+str(zSlices[i]),"Digi shape"+histnameaddon,20,0.,20.,2000,0.,2000.)) for i in range(len(zSlices))]) ndigivsnsamp = dict([(i, ROOT.TH2D("ndigisvsnamp"+str(zSlices[i]),"ndigi vs nsamp"+histnameaddon,500,0.,500.,50,0.,50.)) for i in range(len(zSlices))]) nsamperpad = dict([(i, ROOT.TH1D("nsampperpad"+str(zSlices[i]),"samples per pad"+histnameaddon,20,0.,20.)) for i in range(len(zSlices))]) nsamperpadpsa = dict([(i, ROOT.TH1D("nsampperpad psa"+str(zSlices[i]),"samples per pad psa"+histnameaddon,20,0.,20.)) for i in range(len(zSlices))]) ndigiperpad = dict([(i, ROOT.TH1D("digiperpad"+str(zSlices[i]),"digis per pad"+histnameaddon,10,0.,10.)) for i in range(len(zSlices))]) nsamperdigi = dict([(i, ROOT.TH1D("samperdigi"+str(zSlices[i]),"samples per digi"+histnameaddon,20,0.,20.)) for i in range(len(zSlices))]) nsamperdigi1 = dict([(i, ROOT.TH1D("samperdigi1"+str(zSlices[i]),"samples per digi (1 digi)"+histnameaddon,20,0.,20.)) for i in range(len(zSlices))]) nsamperdigi2 = dict([(i, ROOT.TH1D("samperdigi2"+str(zSlices[i]),"samples per digi (2 digi)"+histnameaddon,20,0.,20.)) for i in range(len(zSlices))]) nsnr = dict([(i, ROOT.TH1D("nsnr"+str(zSlices[i]),"signal to noise"+histnameaddon,200,0.,200.)) for i in range(len(zSlices))]) npadspercluster = dict([(i, ROOT.TH1D("npadspercluster3D"+str(zSlices[i]),"pads per cluster 3D"+histnameaddon,100,0.,100.)) for i in range(len(zSlices))]) npadspercluster2D = dict([(i, ROOT.TH1D("npadspercluster2D"+str(zSlices[i]),"pads per cluster 2D"+histnameaddon,100,0.,100.)) for i in range(len(zSlices))]) nresperclsize = dict([(i, ROOT.TH1D("resperclsize"+str(cSlices[i]),"residuals for clustersize "+str(cSlices[i]),400,-10000,10000)) for i in range(len(cSlices))]) nresvsclsize = dict([(i, ROOT.TH2D("hresvsclsize"+str(zSlices[i]),"residual vs clustersize "+str(zSlices[i]),200,-10000.,10000.,30,0.,30.)) for i in range(len(zSlices))]) if cosmic == 1 : resname="Residuals X (cm)" if sim == 1 or cosmic == 0: resname = "Residuals (cm)" print resname residuals = dict([(i, ROOT.TH1D("StatsResX Z:"+str(zSlices[i]),resname, 400,-10000,10000)) for i in range(len(zSlices))]) residualsZ = dict([(i, ROOT.TH1D("StatsResZ Z:"+str(zSlices[i]),"Residuals Z", 400,-10000,10000)) for i in range(len(zSlices))]) nRresiduals = dict([(i, ROOT.TH1D("StatsResX R:"+str(rSlices[i]),resname, 400,-10000,10000)) for i in range(len(zSlices))]) nRresidualsZ = dict([(i, ROOT.TH1D("StatsResZ R:"+str(rSlices[i]),"Residuals Z", 400,-10000,10000)) for i in range(len(zSlices))]) # }}} filecount=0 print "filelist",filelist print "filenumbers:",filenumber for file in filelist : sigmas = [None]*2 for i in range(2): sigmas[i] = [None]*6 for j in range(6): sigmas[i][j] = [None]*16 for k in range(16) : sigmas[i][j][k]=[None]*76 if sim == 0: pedfnum=pedlist.keys() print pedfnum pedfnum.sort() pedfnum.reverse() print pedfnum,len(pedfnum) foundped = 0 for i in range(len(pedfnum)): print int(pedfnum[i])+1 if int(filenumber[filecount]) >= int(pedfnum[i]) : pedfilenumber = pedfnum[i] print pedfnum[i] print pedlist[str(pedfnum[i])] foundped = 1 break filecount+=1 if foundped == 0 and sim != 1: print "no pedfile found!" sys.exit() for pedfile in glob.glob("./pedestals/*"): # if pedfile.find(pedlist[str(filenumber)]) >= 0 : if pedfile.find(pedlist[str(pedfilenumber)]) >= 0 : srcid = int(pedfile[pedfile.find("Src_")+4:pedfile.find("Src_")+7]) adcid = int(pedfile[pedfile.find("ADC")+4:pedfile.find("ADC")+5]) print "srcid=",srcid," ADC=",adcid fped = open(pedfile,"r") chipid = -1 for line in fped : if line.find("#")>=0 or line.find("Chip")>=0: chipid += 1 channelid = -1 continue channelid += 1 words = line.split() sigma = float(words[2]) sigmas[srcid-802][adcid][chipid][channelid]=sigma if Rfile.GetName() != file : print file Rfile = ROOT.TFile.Open(file,"read") # Rfile.Recover() tree = Rfile.Get("cbmsim") if tree == None : print "No cbmsim in file", file Rfile.Close() continue tree.SetBranchStatus("*", 0) tree.SetBranchStatus("TpcCluster.*", 1) # tree.SetBranchStatus("PndTrackTpc.*",1) tree.SetBranchStatus("TrackFitStat_0.*",1) # tree.SetBranchStatus("TrackFitStat.*",1) tree.SetBranchStatus("TrackPostFit.*",1) if sim == 2 : mcfile = file[0:file.find("reco.root")] mcfile += "mc.root" mcRfile = ROOT.TFile.Open(mcfile,"read") mctree = mcRfile.Get("cbmsim") mctree.SetBranchStatus("*",0) mctree.SetBranchStatus("MCTrack.*",1) if raw == 1 : rawfile = file[0:file.find("reco.root")] rawfile += "raw.root" rawRfile = ROOT.TFile.Open(rawfile,"read") rawtree = rawRfile.Get("cbmsim") rawtree.SetBranchStatus("*",0) rawtree.SetBranchStatus("TpcSample",1) filled = 0 event = 0 print nEvents,"events to process" #start loop over tree actev = 0 for e in tree : # {{{ cuts dropout = 0 if nEvents < event: print "now should be a break" break if raw == 1 : rawtree.GetEntry(event) nSamples = rawtree.TpcSample.GetEntriesFast() print nSamples," samples" donesamples = 0 totamp = 0 # nTracks = e.PndTrackTpc.GetEntriesFast() nTracks = e.TrackPostFit.GetEntriesFast() nCluster = e.TpcCluster.GetEntriesFast() hdropout.Fill(1.) # cut on single tracks if cosmic == 1 : if nTracks > 1 : dropout += 3 if nTracks == 0 : for cl in e.TpcCluster : for dig in range(cl.nDigi()) : digi = cl.getDigi(dig) for sam in range(digi.nSample()) : sample = digi.getSample(sam) hnoiseamp.Fill(sample.amp()) dropout += 4 # cut on parallelity of tracks and curvature numpartracks = 0 if 1 : for tr in e.TrackPostFit : endcap = 0 theta=tr.getMom().Theta()*ROOT.TMath.RadToDeg() phi = tr.getMom().Phi()*ROOT.TMath.RadToDeg() htheta.Fill(theta) hphi.Fill(phi) hmom.Fill(tr.getMom().Mag()) cand = tr.getCand() hitIds = cand.GetHitIDs(4) if tr.getMom().Mag() < 0.004 and momc == 1: dropout += 40 if nTracks == 1 : hthetacut.Fill(theta) hphitrackcut.Fill(phi) if ((theta < 85 and theta > 5) or (theta > 95 and theta < 175)) and thec ==1 and partrack == 1 : # if theta<90: numpartracks+=1 dropout += 10 if (abs(phi) < 80 or abs(phi) > 100) and phic == 1 and partrack ==1 : dropout += 20 cand = tr.getCand() if cand.getCurv() > 0.05 and curv == 1 and partrack ==1 : hthetacurvcut.Fill(theta) hphicurvcut.Fill(phi) dropout += 80 length=0 poscounter=0 for hitId in hitIds : cl = e.TpcCluster.At(hitId) pos=cl.pos() if pos[2]>driftcut and driftcut>-99: endcap = 1 if poscounter>0: length+=(lastpos-pos).Mag() lastpos=pos poscounter+=1 r=math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]) if cand.getCurv() > 0.05: hclsizeradiusccut.Fill(r,cl.size()) for dig in range(cl.nDigi()) : digi = cl.getDigi(dig) if digi.amp() > 1700 and sat == 1 and sim==0 : dropout += 200 if endcap == 0 and driftcut > -99: print "wrong, eveything is wrong!" dropout += 300 hthevslen.Fill(theta,length) htracklen.Fill(length) if dropout > 0 : hdropout.Fill(dropout) continue # }}} hdropout.Fill(2.) hthetagood.Fill(theta) hphigood.Fill(phi) htracklengood.Fill(length) if event % 10 == 0 : print "Event:",event if floop == 0 : for tr in e.TrackPostFit : cand = tr.getCand() hitIds = cand.GetHitIDs(4) for hitId in hitIds: # print hitId cl = e.TpcCluster.At(hitId) hslicehist.Fill(cl.pos()[2]) event+=1 continue trkindex = -1 for tr in e.TrackPostFit : trkindex += 1 resv = std.vector # {{{ get residuals numHits = tr.getNumHits() if numHits >= minNumHits : tfs = e.TrackFitStat_0.At(trkindex) resUa = tfs.GetResU() resVa = tfs.GetResV() length = 0 mom = tfs.GetMom() mom.SetMag(1.) z = tfs.GetHitPositionsZ() x = tfs.GetHitPositionsX() y = tfs.GetHitPositionsY() clsizes = tfs.GetClSizes() clsizes2D = tfs.Get2DClSizes() # build orthogonal vector to fix plane - res x mom vecX = ROOT.TVector3(1,0,0) vecZ = ROOT.TVector3(0,0,1) # define new coordinate system for this track u = mom.Cross(vecX) v = mom.Cross(u) resX = tfs.GetResX() resY = tfs.GetResY() resZ = tfs.GetResZ() for i in range(resUa.size()) : resU = resUa[i] resV = resVa[i] resX[i]*=10000 resY[i]*=10000 resZ[i]*=10000 hresX.Fill(resX[i]) hresY.Fill(resY[i]) hresZ.Fill(resZ[i]) res = ROOT.TVector3(resX[i],resY[i],resZ[i]) resUz = res.Dot(u) resVx = res.Dot(v) # change values to mum resVx*=10000 resV*=10000 resU*=10000 resUz*=10000 if cosmic == 1 : # hallresi.Fill(resVx) hallresi.Fill(resX[i]) if cosmic == 0 : hallresi.Fill(resU) hallresi.Fill(resV) hpadspercluster.Fill(clsizes[i]) hpadspercluster2D.Fill(clsizes2D[i]) for j in range(len(zSlices)) : if z[i] < zSlices[j] : if cosmic == 1 : # residuals[j].Fill(resVx) # nresvsclsize[j].Fill(resVx,clsizes[i]) residuals[j].Fill(resX[i]) residualsZ[j].Fill(resZ[i]) nresvsclsize[j].Fill(resX[i],clsizes[i]) if cosmic == 0: residuals[j].Fill(resU) residuals[j].Fill(resV) nresvsclsize[j].Fill(resU,clsizes[i]) nresvsclsize[j].Fill(resV,clsizes[i]) npadspercluster[j].Fill(clsizes[i]) npadspercluster2D[j].Fill(clsizes2D[i]) break r = math.sqrt(x[i]*x[i]+y[i]*y[i]) hclsizeradius.Fill(r,clsizes[i]) hclsizeradius2D.Fill(r,clsizes2D[i]) hclsizevsradius.Fill(r,clsizes[i]) hclsizevsradius2D.Fill(r,clsizes2D[i]) for j in range(len(rSlices)): if r < rSlices[j]: if cosmic==1: nRresiduals[j].Fill(resX[i]) nRresidualsZ[j].Fill(resZ[i]) break for j in range(len(cSlices)): if clsizes[i] <= cSlices[j] : if cosmic == 1 : # nresperclsize[j].Fill(resVx) nresperclsize[j].Fill(resX[i]) if cosmic == 0: nresperclsize[j].Fill(resU) nresperclsize[j].Fill(resV) break if i>0 : dx=x[i]-x[i-1] dy=y[i]-y[i-1] dz=z[i]-z[i-1] length += math.sqrt(dx*dx+dy*dy+dz*dz) # print "dx=",dx," dy=",dy," dz=",dz # print "length= ", length trklendist.Fill(length) # end loop over tracks # trknumdist.Fill(trkindex+1) # }}} cand = tr.getCand() if debug == 1: print "printing the track candidate" cand.Print() if sim == 0: detids = cand.GetDetIDs() hitIds = cand.GetHitIDs(4) # was former 4 elif sim == 1: hitIds = cand.GetHitIDs(10) # print "detids:",detids # for ids in range(detids.size()): # print detids[ids] # print "end detids" hclusterpertrack.Fill(hitIds.size()) if length > 0: hclusterpercm.Fill(hitIds.size()/length) #start loop over cluster if debug == 1: print "number of cluster:",hitIds.size() # print "det ids:", for hitId in hitIds : cl = e.TpcCluster.At(hitId) #loop over digis samthrcount = 0 digicount = 0 samindig = 0 samplesoverthr = [] nofofdigis = [] maxclsample = 0 hxy.Fill(cl.pos()[0],cl.pos()[1]) r=ROOT.TMath.sqrt(cl.pos()[0]*cl.pos()[0]+cl.pos()[1]*cl.pos()[1]) # if endcap == 1 : hrz.Fill(cl.pos()[2],cl.pos()[1]) # if cl.nPad()>0: # # for zcut in range(len(zSlices)) : # if cl.pos()[2] < zSlices[zcut] : # break if debug == 1: print "number of digis:",cl.nDigi() for dig in range(cl.nDigi()) : digi = cl.getDigi(dig) if digi.amp() > samthresh : nofofdigis.append([]) nofofdigis[digicount].append(digi.padId()) nofofdigis[digicount].append(digi.amp()) nofofdigis[digicount].append(digi.t()) nofofdigis[digicount].append(digi.nSample()) hsamperdigi.Fill(digi.nSample()) digicount += 1 for zcut in range(len(zSlices)) : if cl.pos()[2] < zSlices[zcut] : nsamperdigi[zcut].Fill(digi.nSample()) if cl.nDigi() == 1 : nsamperdigi1[zcut].Fill(digi.nSample()) if cl.nDigi() > 1 : nsamperdigi2[zcut].Fill(digi.nSample()) break else : continue #loop over samples maxsample = 0 intamp = 0 bigger511 = 0 isamp=0 samparr = numpy.zeros(20) maxsamplei = 0 digithrsampl = 0 snr = 0 if digi.nSample()>0: for sam in range(digi.nSample()) : sample = digi.getSample(sam) isamp += 1 if sample.amp() >= samthresh : hSampleTot.Fill(sample.amp()) samplesoverthr.append([]) samplesoverthr[samthrcount].append(sample.padId()) samplesoverthr[samthrcount].append(sample.amp()) samplesoverthr[samthrcount].append(sample.t()) samthrcount += 1 digithrsampl += 1 hAllDigiSamples.Fill(sample.t(),sample.amp()) if isamp < 20 : samparr[isamp]=sample.amp() # if sample.t()>511: # print sample.t() if sample.amp() > maxsample: maxsample = sample.amp() maxsourceid = sample.sourceId() maxadcid = sample.adcId() maxchipid = sample.chipId() maxchannelid = sample.channelId() maxsamplei = isamp intamp += sample.amp() if sample.t()>511 : bigger511 = 1 if maxsample > maxclsample : maxclsample = maxsample maxclsourceid = maxsourceid maxcladcid = maxadcid maxclchipid = maxchipid maxclchannelid = maxchannelid # print digi.nSample(), digithrsampl if digithrsampl > 0 : hnsamperpadpsa.Fill(digithrsampl) samindig += digithrsampl for zcut in range(len(zSlices)) : if cl.pos()[2] < zSlices[zcut] : nsamperpadpsa[zcut].Fill(digithrsampl) break if bigger511 == 1 : for sam in range(digi.nSample()) : sample = digi.getSample(sam) hsamplebigger511.Fill(sample.t(),sample.amp()) if intamp > 7000 : continue hintsampleamp.Fill(intamp) samplesoverthr.sort(key=sort_inner2) samplesoverthr.sort(key=sort_inner) if maxclsample > 0 : if sim == 0: somesigma = sigmas[maxclsourceid-802][maxcladcid][maxclchipid][maxclchannelid] if somesigma > 0: snr = maxclsample/somesigma else: print "no noise info for:",maxclsourceid-802,maxcladcid,maxclchipid,maxclchannelid snr = 0 else : snr = (maxclsample*4)/20 if snr > 8 : hsnr.Fill(snr) hClusterMax.Fill(maxclsample) for zcut in range(len(zSlices)) : if cl.pos()[2] < zSlices[zcut] : nsnr[zcut].Fill(snr) break meansamperpad = 0 digiperpad = 0 nofpads = 1 if len(samplesoverthr) > 0 : oldpad = samplesoverthr[0][0] samplesoverthr.append([]) samplesoverthr[len(samplesoverthr)-1].append(-99) samplesoverthr[len(samplesoverthr)-1].append(-99) samplesoverthr[len(samplesoverthr)-1].append(-99) padamp = 0 nummax = 0 sampperpad = 0 up = 1 k=0 while k < len(samplesoverthr) : actpad = samplesoverthr[k][0] if oldpad == actpad and k > 0 and samplesoverthr[k][2] == samplesoverthr[k-1][2] : samplesoverthr.pop(k) continue if actpad != oldpad : nofpads += 1 hnsamperpad.Fill(sampperpad) digiperpad = 0 for d in range(len(nofofdigis)) : if nofofdigis[d][0] == oldpad : digiperpad += 1 hdigiperpad.Fill(digiperpad) for zcut in range(len(zSlices)) : if cl.pos()[2] < zSlices[zcut] : nsamperpad[zcut].Fill(sampperpad) ndigiperpad[zcut].Fill(digiperpad) break meansamperpad =+ sampperpad sampperpad = 0 # if digiperpad == 5 : if 0 : spec = ROOT.TSpectrum() spec.SetResolution(12) nofpeaks = spec.Search(hfullpad,1) c98.cd() hfullpad.SetAxisRange(hfullpad.GetMean()-10,hfullpad.GetMean()+10) hfullpad.Draw() c98.Update() print "found ",nummax ," maxima so far" nix = raw_input("press enter") hfullpad.Reset() hfullpad.SetAxisRange(0.,511.) if samplesoverthr[k][0] >= 0 : hfullpad.Fill(samplesoverthr[k][2],samplesoverthr[k][1]) if padamp > samplesoverthr[k][1] and up == 1 : nummax += 1 up = 0 if padamp < samplesoverthr[k][1] and up == 0 : up = 1 padamp = samplesoverthr[k][1] sampperpad += 1 oldpad = actpad k += 1 meansamperpad /= nofpads meansamperdig = samindig/cl.nDigi() hnpadsndigi.Fill(nofpads-cl.nDigi()) hndigivsnsamp.Fill(len(samplesoverthr)-1,len(nofofdigis)) for zcut in range(len(zSlices)) : if cl.pos()[2] < zSlices[zcut] : ndigivsnsamp[zcut].Fill(len(samplesoverthr),len(nofofdigis)) break #end loop over cluster if 0 : # if digiperpad > 1: if endcap == 1: c99.cd() hxy.Draw("colz") c99.Update() c98.cd() hrz.Draw("colz") c98.Update() # print "snr:",snr # print "dip:",cand.getDip(),"curv:",cand.getCurv() bla=raw_input(str(theta)) hxy.Reset() hrz.Reset() event = event+1 #end loop over tree # Rfile.Close() totev +=event htheta.Draw() # nix = raw_input ("Enter to continue" ) #end loop over files if floop == 0: stotev = hslicehist.GetEntries() print "total events:", stotev stotev /= 6 print "slice events:",stotev slicesum = 0 slicen = 0 for i in range(0,79): slicesum += hslicehist.GetBinContent(i) if slicesum + hslicehist.GetBinContent(i+1) >= stotev and slicen<5: print "sum:",slicesum," slice:",i-60 if slicedummy == 0: hslices.Fill(slicen,i-60) elif slicedummy == 1 : hslices.Fill(slicen,slicen*13-60) zSlices.append(i-60) slicen+=1 slicesum = 0 zSlices.append(10) hslices.Fill(5,10) print "zslices:",zSlices print "Done Loop over files" if outfname == "" : outfolder = "anafiles/" outfname = outfolder + "sampleAnaOut" if filecount > 1: for i in range(len(runList)): outfname += str(runList[i]) outfname += "." else: outfname += runs outfname += "." if floop == 0: outfname += "slices_" outfname += str(totev) outfname += "evts" if part>=0: outfname +="_part" outfname +=part #outfname += ".thr" #outfname += str(samthresh) outfname += ".root" outfile = ROOT.TFile(outfname,"recreate") outfile.cd() toprint = [] prname = [] if nplot == 0 : if 0 : c2 = ROOT.TCanvas("c2","one digi") c4 = ROOT.TCanvas("c4","digi width") c5 = ROOT.TCanvas("c5","digi rms") c6 = ROOT.TCanvas("c6","max amplitude") c7 = ROOT.TCanvas("c7","sliced width") c8 = ROOT.TCanvas("c8","int sample amp") c9 = ROOT.TCanvas("c9","width x amp") c10 = ROOT.TCanvas("c10","chisquare") c11 = ROOT.TCanvas("c11","bigger 511") c12 = ROOT.TCanvas("c12","siz vs sig") c13 = ROOT.TCanvas("c13","chii vs sig") c14 = ROOT.TCanvas("c14","digi") c15 = ROOT.TCanvas("c15","digisliced 2d") c16 = ROOT.TCanvas("c16","digisliced") c3 = ROOT.TCanvas("c3","all samples") c1 = ROOT.TCanvas("c1","sample tot") c17 = ROOT.TCanvas("c17","ndigi vs nsamp") c18 = ROOT.TCanvas("c18","ndigi vs nsamp sliced") c19 = ROOT.TCanvas("c19","samperpad slice") c20 = ROOT.TCanvas("c20","samperpad") c21 = ROOT.TCanvas("c21","samperpad psa") c22 = ROOT.TCanvas("c22","samperpad psa slice") c23 = ROOT.TCanvas("c23","pad - dig") c24 = ROOT.TCanvas("c24","digis per pad") c25 = ROOT.TCanvas("c25","digis per pad slice") c26 = ROOT.TCanvas("c26","amplitude of noise") c27 = ROOT.TCanvas("c27","samples per digi") c28 = ROOT.TCanvas("c28","samples per digi sliced") c29 = ROOT.TCanvas("c29","snr") c30 = ROOT.TCanvas("c30","snr sliced") c31 = ROOT.TCanvas("c31","max cl sam") c32 = ROOT.TCanvas("c32","cl per track") c33 = ROOT.TCanvas("c33","multi purpose") c34 = ROOT.TCanvas("c34","residuals") c35 = ROOT.TCanvas("c35","multi purpose 2") c36 = ROOT.TCanvas("c36","pads per cluster sliced") c37 = ROOT.TCanvas("c37","res per clsize") c38 = ROOT.TCanvas("c38","res vs clsize") c39 = ROOT.TCanvas("c39","resZ vs z") c40 = ROOT.TCanvas("c40","resX vs R") c41 = ROOT.TCanvas("c41","resZ vs R") c42 = ROOT.TCanvas("c42","Cluster per cm") c43 = ROOT.TCanvas("c43","pads pes cluster 2D sliced") c44 = ROOT.TCanvas("c44","multi purpose 3") c1.cd() hSampleTot.Draw() toprint.append(c1) prname.append("sample_tot") c3.cd() hAllDigiSamples.Draw() hAllDigiSamples.Write() toprint.append(c3) prname.append("alldigisamples") if 0 : c2.cd() hOneDigi.Draw() c4.cd() hDigiWidth.Draw() c5.cd() hDigiRMS.Draw() c6.cd() hDigiMaxS.Draw() c7.Divide(3,2) c9.Divide(3,2) c15.Divide(3,2) c18.Divide(3,2) c19.Divide(3,2) c22.Divide(3,2) c25.Divide(3,2) c28.Divide(3,2) c30.Divide(3,2) c34.Divide(3,2) c36.Divide(3,2) c38.Divide(3,2) c39.Divide(3,2) for i in range(len(zSlices)): if 0 : c7.cd(i+1) samwidths[i].Draw() samwidths[i].Write() c9.cd(i+1) samwidthxamp[i].Draw() samwidthxamp[i].Write() c15.cd(i+1) digissliced2d[i].Draw("colz") digissliced2d[i].Write() c16.cd(i+1) digissliced[i].Draw() digissliced[i].Write() c18.cd(i+1) c18.GetPad(i+1).SetLogz() ndigivsnsamp[i].Write() ndigivsnsamp[i].Draw("colz") c19.cd(i+1) nsamperpad[i].Write() nsamperpad[i].Draw() c22.cd(i+1) nsamperpadpsa[i].Write() nsamperpadpsa[i].Draw() c25.cd(i+1) ndigiperpad[i].Write() ndigiperpad[i].Draw() c28.cd(i+1) nsamperdigi[i].Write() nsamperdigi[i].Draw() nsamperdigi1[i].Write() nsamperdigi1[i].SetLineColor(ROOT.kRed) nsamperdigi1[i].Draw("same") nsamperdigi2[i].Write() nsamperdigi2[i].SetLineColor(ROOT.kBlue) nsamperdigi2[i].Draw("same") c30.cd(i+1) nsnr[i].Write() nsnr[i].Draw() c34.cd(i+1) residuals[i].Write() residuals[i].Draw() c36.cd(i+1) npadspercluster[i].Draw() npadspercluster[i].Write() c38.cd(i+1) nresvsclsize[i].Draw("colz") nresvsclsize[i].Write() c39.cd(i+1) residualsZ[i].Write() residualsZ[i].Draw() c37.Divide(3,3) for i in range(len(cSlices)): c37.cd(i+1) nresperclsize[i].Draw() if floop==1: nresperclsize[i].Write() c40.Divide(3,2) toprint.append(c40) prname.append("residuals_vs_R") for i in range(len(rSlices)): c40.cd(i+1) nRresiduals[i].Draw() if floop==1: nRresiduals[i].Write() c41.Divide(3,2) toprint.append(c41) prname.append("Zresiduals_vs_R") for i in range(len(rSlices)): c41.cd(i+1) nRresidualsZ[i].Draw() if floop==1: nRresidualsZ[i].Write() c43.Divide(3,2) toprint.append(c43) prname.append("pads_per_culster_2D_sliced") for i in range(len(zSlices)): c43.cd(i+1) npadspercluster2D[i].Draw() if floop==1: npadspercluster2D[i].Write() toprint.append(c19) prname.append("samples per pad sliced") toprint.append(c22) prname.append("samples per pad psa sliced") toprint.append(c25) prname.append("digis_per pad sliced") toprint.append(c28) prname.append("samples_per_digi_sliced") toprint.append(c30) prname.append("snr_sliced") toprint.append(c36) prname.append("pads_per_cluster") toprint.append(c37) prname.append("res_per_clsize") toprint.append(c38) prname.append("res_vs_clsize_sliced") toprint.append(c39) prname.append("resZ_vs_z") if 0 : c8.cd() hintsampleamp.Draw() c10.cd() hchisquare.Draw() c11.cd() hsamplebigger511.Draw() c12.cd() hsigmavsnofsampl.Draw() c13.cd() hsigmavschii.Draw() c14.cd() hdigi.Draw() c17.cd() c17.SetLogz() hndigivsnsamp.Draw("colz") toprint.append(c17) prname.append(hndigivsnsamp.GetTitle()) c20.cd() hnsamperpad.Draw() toprint.append(c20) prname.append(hnsamperpad.GetTitle()) c21.cd() hnsamperpadpsa.Draw() toprint.append(c21) prname.append(hnsamperpadpsa.GetTitle()) c23.cd() hnpadsndigi.Draw() toprint.append(c23) prname.append(hnpadsndigi.GetTitle()) c24.cd() hdigiperpad.Draw() toprint.append(c24) prname.append(hdigiperpad.GetTitle()) c26.cd() hnoiseamp.Draw() toprint.append(c26) prname.append(hnoiseamp.GetTitle()) c27.cd() hsamperdigi.Draw() toprint.append(c27) prname.append(hsamperdigi.GetTitle()) c29.cd() hsnr.Draw() toprint.append(c29) prname.append(hsnr.GetTitle()) c31.cd() hClusterMax.Draw() toprint.append(c31) prname.append(hClusterMax.GetTitle()) c32.cd() hclusterpertrack.Draw() toprint.append(c32) prname.append(hclusterpertrack.GetTitle()) c33.Divide(3,3) c33.cd(1) thetaleg=ROOT.TLegend(0.7,.7,0.9,.9,"","NDC") htheta.Draw() thetaleg.AddEntry(htheta,"All") hthetacut.SetLineColor(2) hthetacut.Draw("same") thetaleg.AddEntry(hthetacut,"# tracks = 1") hthetacurvcut.SetLineColor(3) hthetacurvcut.Draw("same") thetaleg.AddEntry(hthetacurvcut,"curvature cut") hthetagood.SetLineColor(4) hthetagood.Draw("same") thetaleg.AddEntry(hthetagood,"survivors") thetaleg.SetFillColor(0) thetaleg.Draw() c33.cd(2) phileg=ROOT.TLegend(0.7,.7,0.9,.9,"","NDC") hphi.Draw() phileg.AddEntry(hphi,"all") hphitrackcut.SetLineColor(2) hphitrackcut.Draw("same") phileg.AddEntry(hphitrackcut,"#tracks = 1") hphicurvcut.SetLineColor(3) hphicurvcut.Draw("same") phileg.AddEntry(hphicurvcut,"curvature cut") hphigood.SetLineColor(4) hphigood.Draw("same") phileg.AddEntry(hphigood,"survivors") phileg.SetFillColor(0) phileg.Draw() c33.cd(3) hallresi.Draw() c33.cd(4) hdropout.Draw() c33.cd(5) hxy.Draw("colz") c33.cd(6) hrz.Draw("colz") c33.cd(7) hmom.Draw() c33.cd(8) hslicehist.Draw() c33.cd(9) hslices.Draw() c35.Divide(3,3) c35.cd(1) c35.GetPad(1).SetLogz(1) hthevslen.Draw("colz") c35.cd(2) htracklen.Draw() htracklengood.SetLineColor(2) htracklengood.Draw("same") c35.cd(4) hpadspercluster.Draw() c35.cd(5) hpadspercluster2D.Draw() c35.cd(7) hresX.Draw() c35.cd(8) hresY.Draw() c35.cd(9) hresZ.Draw() c44.Divide(3,3) c44.cd(1) c44.GetPad(1).SetLogz(1) hclsizevsradius.Draw("colz") c44.cd(2) c44.GetPad(2).SetLogz(1) hclsizevsradius2D.Draw("colz") c44.cd(3) hclsizeradius.Draw() c44.cd(4) hclsizeradius2D.Draw() c44.cd(5) hclsizeradiusccut.Draw() toprint.append(c33) prname.append("multipurpose") toprint.append(c35) prname.append("multipurpose2") toprint.append(c44) prname.append("multipurpose3") c42.cd() hclusterpercm.Draw() toprint.append(c42) prname.append("hcluster_per_cm") if floop == 1: hSampleTot.Write() hOneDigi.Write() hAllDigiSamples.Write() hDigiWidth.Write() hDigiRMS.Write() hDigiMaxS.Write() hintsampleamp.Write() hchisquare.Write() hsamplebigger511.Write() hsigmavsnofsampl.Write() hsigmavschii.Write() hdigi.Write() hndigivsnsamp.Write() hnsamperpad.Write() hnsamperpadpsa.Write() hnpadsndigi.Write() hdigiperpad.Write() hnoiseamp.Write() hsamperdigi.Write() hsnr.Write() hClusterMax.Write() hclusterpertrack.Write() htheta.Write() hthetacut.Write() hthetacurvcut.Write() hphi.Write() hphitrackcut.Write() hphicurvcut.Write() hallresi.Write() hdropout.Write() hxy.Write() hrz.Write() hmom.Write() hslicehist.Write() hslices.Write() hthevslen.Write() htracklen.Write() htracklengood.Write() hpadspercluster.Write() hpadspercluster2D.Write() hclusterpercm.Write() if printc == 1 : for i in range(len(toprint)): prname[i]=prname[i].replace(" ","_") toprint[i].SaveAs(pdir+prname[i]+".pdf") if floop ==0: hslicehist.Write() hslices.Write() print "saved histos in",outfname outfile.Close() if batchmode == 0 : nix = raw_input ("Enter to finish" )