import ROOT,glob, math, sys, os, numpy, argparse sys.path.append('$PANDAPATH/macro/tpc/FOPI/mberger') #from ROOT import std from functions import sort_inner, sort_inner2, Drawh, Divideh, fillcutred, plotfullevent # {{{ argument parsing parser=argparse.ArgumentParser(description='Analyze the cosmics. Either give the path and run numbers or a specific file.') parser.add_argument("--path",help='the path where all the files to analyze are',type=str, default="./outfiles") parser.add_argument("--runs",help='which runs to analyze',type=str,default='reco') parser.add_argument("--filenames",help='name of the files to analyze',type=str,nargs='+',default="") parser.add_argument("--outfname",help='name of the output file',type=str,default="") parser.add_argument("--pdir", help="The path where to print the pictures. If defined pictures will be printed", type=str,default="./pics/tmp") parser.add_argument("--printc", help="Print the pictures either to ./pics/tmp or to the folder given with --pdir", action="store_true") parser.add_argument("--nEvents",help='number of events to analyze',type=int, default=-1) parser.add_argument("--sim", help='analyze simulated events. This changes the calculation of the SNR. Its mandatory to set this in case of simulations.',action='store_true') parser.add_argument("--raw",help='enable the MC digitizer data',type=str,default="") parser.add_argument("--chi",help='set the chi^2 value for a cut. (not used)',type=float,default=1.5) parser.add_argument("--cosmic",help='the cosmics cut. means only 1 track per event', action='store_const',const=1,default=0) parser.add_argument("--sampling",help='set the used sampling frequency (MHz)', type=float,default=15.55) parser.add_argument("--samthr",help='set the threshold for the sample amplitude', type=float,default=0.) parser.add_argument("--digithr",help='set the threshold for the digi amplitude', type=float,default=0.) parser.add_argument("--cluthr",help='set the threshold for the cluster amplitude', type=float,default=0.) parser.add_argument("--plot",help='plot everything onto the screen',action='store_true') parser.add_argument("--curv",help='enable the curvature cut (deprecated)', action='store_const',const=1,default=0) parser.add_argument("--thec",help='enable the theta cut', action='store_const',const=1,default=0) parser.add_argument("--phic",help='enable the phi cut', action='store_const',const=1,default=0) parser.add_argument("--fcut",help='DISABLE track fit flag cut' ,action='store_const',const=0,default=1) parser.add_argument("--cpcut",help='enable cluster per track cut (10)', action='store_const',const=1,default=0) parser.add_argument("--radcut",help='enable radius cut (7)', action='store_const',const=1,default=0) parser.add_argument("--fila",help='enable first and last cluster cut. first and last cluster of a track are discarded',action='store_true') parser.add_argument("--clmean",help='enable cut on cluster mean amplitude (400)', action='store_const',const=1,default=0) parser.add_argument("--chi2cdfc",help='enable chi^2 desity function cut', action='store_const',const=1,default=0) parser.add_argument("--batch",help='enable batch mode -> no plotting or what so ever',action='store_true') parser.add_argument("--sat",help='DISABLE the digi saturation cut. This cut removes all cluster with maxSampleInCluster > 1700',action='store_const',const=0,default=1) parser.add_argument("--momc",help='enable the momentum cut',action='store_const',const=1,default=0) parser.add_argument("--debug",help='turn debug output on',action='store_true') parser.add_argument("--slicef",help='use a special file for the zslices', type=str,default="") parser.add_argument("--single",help='some analysis with plots eventwise',action='store_true') parser.add_argument("--pre",help='use information from "pre" branches (before reclustering).', action='store_true') parser.add_argument("--detid",help='use the given detector id',type=int,default=-1) parser.add_argument("--plotsingle",help='plot a single histogram',type=str,default=[""],nargs='+') parser.add_argument("--hlp",help='print help',action='store_true') args=parser.parse_args() if args.hlp: parser.print_help() exit(0) thrcut=0 if (args.pdir != "./pics/tmp"): args.printc=true if (args.cluthr>0.): thrcut=1 if not args.batch: if (args.thec==0): u=raw_input("Be aware that the theta cut is disabled!") if (args.phic==0): u=raw_input("Be aware that the phi cut is disabled!") if (args.sat==0): u=raw_input("Do you really want to turn off the sample saturation cut?") if (args.fcut==0): u=raw_input("Do you really want to turn off the track fit flag cut?") if args.batch: ROOT.gROOT.SetBatch(True) Rfile=ROOT.TFile() oldpedfilenumber=0 # }}} # {{{ argument printing print "***********MISC***************" print "input file path:\t\t", args.path print "run numbers:\t\t\t",args.runs print "files to process:\t\t",args.filenames print "out file name:\t\t\t",args.outfname print "number of events:\t\t",args.nEvents print "chi2 cut for nothing:\t\t", args.chi print "sampling freq:\t\t\t",args.sampling print "print directory:\t\t",args.pdir print "detid:\t\t\t\t",args.detid print "**********THRESHOLDS**********" print "sample amp threshold:\t\t",args.samthr print "digi amp threshold:\t\t",args.digithr print "*************CUTS*************" print "cosmic flag (1 track cut):\t",args.cosmic print "cluster mean amp cut:\t\t",args.clmean print "phi cut flag:\t\t\t",args.phic print "theta cut flag:\t\t\t",args.thec print "momentum cut flag:\t\t",args.momc print "curvature cut:\t\t\t",args.curv #print "cluster saturation cut flag:\t",args.clsat print "sample saturation cut flag:\t",args.sat print "fit cut:\t\t\t",args.fcut print "cluster/track cut:\t\t",args.cpcut print "radius cut:\t\t\t",args.radcut print "first last cut:\t\t\t",args.fila print "*************FLAGS************" print "print flag:\t\t\t",args.printc print "batchmode flag:\t\t\t",args.batch print "plotting flag:\t\t\t",args.plot print "debug flag: \t\t\t",args.debug print "sim flag:\t\t\t",args.sim print "pre:\t\t\t\t",args.pre if not args.batch: nix = raw_input("Happy with this settings? (y/n)") if nix=='n': exit(0) #finished argument parsing # }}} # {{{ create runlist runList=[] colindex = args.runs.find(",") if colindex > 0: args.runs = args.runs.split(",") for i in range(len(args.runs)) : dashIndex = args.runs[i].find("-") if dashIndex > 0 : args.runs[i] = args.runs[i].split("-") for j in range(int(args.runs[i][0]), int(args.runs[i][1])+1) : runList.append(j) if dashIndex < 0 and len(args.runs) > 1 : runList.append(int(args.runs[i])) else : dashIndex = args.runs.find("-") if dashIndex > 0 : args.runs = args.runs.split("-") runList = range(int(args.runs[0]), int(args.runs[1])+1) if dashIndex < 0 and len(args.runs) > 1 : runList = [args.runs] print "runs to process:", runList # }}} # {{{ load some libraries and set common options 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+") # }}} # {{{ files to process and pedestal files files= glob.glob(args.path + "/*reco*.root") files.sort() # {{{loop over files to create filelist: totev = 0 totgoodev = 0 filecount = 0 firstloop = 1 filelist = [] filenumber = [] print "Start loop over files" if args.filenames=="": for file in files : if len(runList) >= 1 : found = 0 for i in runList : if file.find(str(i)) >= 0 : found = 1 filecount += 1 filelist.append(file) filenumber.append(i) #this filenumber has to be generated corectly for the args.filenames else: filelist=args.filenames for i in range(len(args.filenames)): num=args.filenames[i].find('runC_') if num>0: num=args.filenames[i][num+5:num+9] filenumber.append(num) print "files to process:",filelist print "run numbers:",filenumber #exit(0) # }}} 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_" } # }}} # {{{ define histos from cosmic_histos import * histsdict=cosmic_histos() histsopt=cosmic_histopt() canvasopt=cosmic_canvasopt() axlabel=cosmic_axlabel() # }}} # {{{ read in slices from slicehisto and define dicts zSlices = [] if args.slicef!="": print "slicename", args.slicef infiles=ROOT.TFile(args.slicef) zSlices = [] histsdict['hslices'] = infiles.Get("hslices") for k in range(1,7): print histsdict['hslices'].GetBinContent(k) zSlices.append(histsdict['hslices']).GetBinContent(k) else: zSlices = ( -48.0, -36.0, -24.0, -12.0, 0.0, 20.0 ) # cosmic ne #zSlices = ( 20, 30, 40, 49, 57, 70 ) # beam for i in range(len(zSlices)): histsdict['hslices'].SetBinContent(i+1,zSlices[i]) cSlices = (3,6,9,12,15,18,21,24,50) rSlices = (6.6,8.2,9.8,11.4,13,15) xySlices = (-10,-5,0,5,10,20) cpSlices = (7,14,21,28,35,42) thSlices = (90,95,100,105,110,120) # {{{ define slice dict hslices={} # {{{ unused samwidths = dict([(i, ROOT.TH1D("digiwidth"+str(zSlices[i]),"Digiwidth",100,0.,10.*(1/args.sampling)*1000)) for i in range(len(zSlices))]) samwidthxamp = dict([(i, ROOT.TH1D("digiwidthxamp"+str(zSlices[i]),"Digi width x amp",1000,0.,1000.)) for i in range(len(zSlices))]) digissliced2d = dict([(i, ROOT.TH2D("digis2d"+str(zSlices[i]),"Digi shape",20,0.,20.,2000,0.,2000.)) for i in range(len(zSlices))]) digissliced = dict([(i, ROOT.TH2D("digis"+str(zSlices[i]),"Digi shape",20,0.,20.,2000,0.,2000.)) for i in range(len(zSlices))]) ndigivsnsamp = dict([(i, ROOT.TH2D("ndigisvsnamp"+str(zSlices[i]),"ndigi vs nsamp",500,0.,500.,50,0.,50.)) for i in range(len(zSlices))]) nsamperpadpsa = dict([(i, ROOT.TH1D("nsampperpad psa"+str(zSlices[i]),"samples per pad psa",20,0.,20.)) for i in range(len(zSlices))]) nsamperdigi1 = dict([(i, ROOT.TH1D("samperdigi1"+str(zSlices[i]),"samples per digi (1 digi)",20,0.,20.)) for i in range(len(zSlices))]) nsamperdigi2 = dict([(i, ROOT.TH1D("samperdigi2"+str(zSlices[i]),"samples per digi (2 digi)",20,0.,20.)) for i in range(len(zSlices))]) # }}} resbins=200 sliceopt={} saxlabel=[] saxlabel.append({}) saxlabel.append({}) 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))]) hslices['residualsX'] = dict([(i, ROOT.TH1D("StatsResX Z:"+str(zSlices[i]),"Residuals X", resbins,-10000,10000)) for i in range(len(zSlices))]) hslices['residualsZ'] = dict([(i, ROOT.TH1D("StatsResZ Z:"+str(zSlices[i]),"Residuals Z", resbins,-10000,10000)) for i in range(len(zSlices))]) hslices['nRresidualsX'] = dict([(i, ROOT.TH1D("StatsResX R:"+str(rSlices[i]),"R Residuals X", resbins,-10000,10000)) for i in range(len(rSlices))]) hslices['nRresidualsZ'] = dict([(i, ROOT.TH1D("StatsResZ R:"+str(rSlices[i]),"R Residuals Z", resbins,-10000,10000)) for i in range(len(rSlices))]) hslices['nresperclsize'] = dict([(i, ROOT.TH1D("resperclsize"+str(cSlices[i]),"residuals for clustersize "+str(cSlices[i]),resbins,-10000,10000)) for i in range(len(cSlices))]) hslices['clustersize'] = dict([(i, ROOT.TH1D("npadspercluster3D"+str(zSlices[i]),"pads per cluster 3D",100,0.,100.)) for i in range(len(zSlices))]) hslices['clustersize2D'] = dict([(i, ROOT.TH1D("npadspercluster2D"+str(zSlices[i]),"pads per cluster 2D",100,0.,100.)) for i in range(len(zSlices))]) hslices['nsnr'] = dict([(i, ROOT.TH1D("nsnr"+str(zSlices[i]),"signal to noise",200,0.,600.)) for i in range(len(zSlices))]) hslices['csnr'] = dict([(i, ROOT.TH1D("csnr"+str(zSlices[i]),"signal to noise (cluster)",200,0.,600.)) for i in range(len(zSlices))]) hslices['nsamperpad'] = dict([(i, ROOT.TH1D("nsampperpad"+str(zSlices[i]),"samples per pad",20,0.,20.)) for i in range(len(zSlices))]) hslices['ndigiperpad'] = dict([(i, ROOT.TH1D("digiperpad"+str(zSlices[i]),"digis per pad",10,0.,10.)) for i in range(len(zSlices))]) hslices['nsamperclust'] = dict([(i, ROOT.TH1D("samperclust"+str(zSlices[i]),"samples per cluster",100,0.,100.)) for i in range(len(zSlices))]) hslices['clustermax'] = dict([(i, ROOT.TH1D("clustermax"+str(zSlices[i]),"max amp in cluster",250,0.,1000.)) for i in range(len(zSlices))]) hslices['clusteramp'] = dict([(i, ROOT.TH1D("clusteramp"+str(zSlices[i]),"cluster amp",250,0.,2000.)) for i in range(len(zSlices))]) hslices['XResidualsX'] = dict([(i, ROOT.TH1D("StatsResX X:"+str(xySlices[i]),"XResiduals X", resbins,-10000,10000)) for i in range(len(xySlices))]) hslices['YResidualsX'] = dict([(i, ROOT.TH1D("StatsResX Y:"+str(xySlices[i]),"YResiduals X", resbins,-10000,10000)) for i in range(len(xySlices))]) hslices['cpResidualsX'] = dict([(i, ROOT.TH1D("StatsResX cp:"+str(cpSlices[i]),"cpResiduals X", resbins,-10000,10000)) for i in range(len(cpSlices))]) hslices['thetaResidualsX']= dict([(i, ROOT.TH1D("StatsResX th:"+str(thSlices[i]),"thResiduals X", resbins,-10000,10000)) for i in range(len(thSlices))]) hslices['hxy res slice']=dict([(i,ROOT.TH2D("hxyres"+str(zSlices[i]),"xy res "+str(zSlices[i]),200,-16.,16.,200,-16.,16.)) for i in range(len(zSlices))]) hslices['hxy slice']=dict([(i,ROOT.TH2D("hxy"+str(zSlices[i]),"xy "+str(zSlices[i]),200,-16.,16.,200,-16.,16.)) for i in range(len(zSlices))]) sliceopt['hxy res slice']="colz" sliceopt['hxy slice']="colz" for j in range(len(zSlices)): hslices['XResidualsX Z'+str(zSlices[j])] = dict([(i, ROOT.TH1D("StatsResX X:"+str(xySlices[i])+"Z"+str(zSlices[j]),"XResiduals X", resbins,-10000,10000)) for i in range(len(xySlices))]) saxlabel[0]['hxy res slice'+str(j)]="X (cm)" saxlabel[1]['hxy res slice'+str(j)]="Y (cm)" saxlabel[0]['hxy slice'+str(j)]="X (cm)" saxlabel[1]['hxy slice'+str(j)]="Y (cm)" # }}} # }}} # {{{ define cut dict cuts={'ntrack':0,'notrack':0,'mom':0,'theta':0,'phi':0,'curv':0,'sat':0,'snr':0,'clamp':0,'fcut':0,'cp':0,'rad':0,'clmean':0,'chi2cdf':0} # }}} # {{{ main loop over the files if args.debug: print "runs",runList print "filelist",filelist print "filenumbers:",filenumber filecount = 0 for file in filelist : # {{{ get pedestal file and read in pedestals # {{{ get pedestal file number olddone=0 if not args.sim: pedfnum=pedlist.keys() pedfnum.sort() pedfnum.reverse() foundped = 0 for i in range(len(pedfnum)): if int(filenumber[filecount]) >= int(pedfnum[i]) : pedfilenumber = pedfnum[i] foundped = 1 print "found a pedestalfile" break filecount+=1 if foundped == 0 and not args.sim: print "no pedfile found!" exit(0) # }}} if pedfilenumber!=oldpedfilenumber or len(filelist)==1: # {{{ reset pedestal array print "reset pedestal array" sigmas = [None]*2 means = [None]*2 for i in range(2): sigmas[i] = [None]*6 means[i] = [None]*6 for j in range(6): sigmas[i][j] = [None]*16 means[i][j] = [None]*16 for k in range(16) : sigmas[i][j][k]=[None]*76 means[i][j][k]=[None]*76 # }}} for pedfile in glob.glob("./pedestals/*"): 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]) mean = float(words[1]) sigmas[srcid-802][adcid][chipid][channelid]=sigma means[srcid-802][adcid][chipid][channelid]=mean oldpedfilenumber=pedfilenumber if args.debug: print "used pedestal file number:",pedfilenumber # }}} # {{{ get root file if Rfile.GetName() != file : tree=None print "file:",file Rfile = ROOT.TFile.Open(file,"read") if type(Rfile)==ROOT.TFile: if Rfile.IsOpen()==ROOT.kTRUE : tree = Rfile.Get("cbmsim") if tree == None : print "No cbmsim in file", file if type(Rfile)==ROOT.TFile: Rfile.Close() continue # }}} # {{{ activate/deactivate branches and get mc or raw files tree.SetBranchStatus("*", 0) if args.pre: tree.SetBranchStatus("TpcPreSPHit.*", 1) tree.SetBranchStatus("TpcPreCluster.*", 1) tree.SetBranchStatus("TpcSample.*",1) tree.SetBranchStatus("TpcDigi.*",1) tree.SetBranchStatus("TrackPrePostFit.*",1) tree.SetBranchStatus("PreTrackFitStat_0.*",1) else: tree.SetBranchStatus("TpcSPHit.*", 1) tree.SetBranchStatus("TpcCluster.*", 1) tree.SetBranchStatus("TpcSample.*",1) tree.SetBranchStatus("TpcDigi.*",1) tree.SetBranchStatus("TrackPostFit.*",1) tree.SetBranchStatus("TrackFitStat_0.*",1) # tree.SetBranchStatus("TrackFitStat.*",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 args.raw!="": rawfile = file[0:file.find("reco.root")] rawfile += "raw.root" # rawRfile = ROOT.TFile.Open(args.raw,"read") # rawtree = rawRfile.Get("cbmsim") # rawtree.SetBranchStatus("*",0) # rawtree.SetBranchStatus("TpcSample",1) # rawtree.SetBranchStatus("TpcDigi",1) tree.AddFriend("cbmsim",args.raw) # }}} # {{{ start loop over tree filled = 0 evt = 0 goodev = 0 actev = 0 doev=tree.GetEntriesFast() if args.nEvents==-1: args.nEvents=doev print args.nEvents,"events to process,",doev ,"in this file" # {{{ stuff for single event if args.single: tcn=ROOT.TCanvas() tcn.cd() tcn.Divide(2,1) # tcn2=ROOT.TCanvas() # read in mapping file fmap=open("./tpc/FOPI/par/padplane_FOPI.dat","r") mapping=[None]*10260 for line in fmap: words=line.split() # print "words:",words if words[0]=='PadPlaneSubDivision': continue # print "pad:",int(words[0])," x:",float(words[3])," y:",float(words[4]) mapping[int(words[0])]=[] mapping[int(words[0])].append(float(words[3])) mapping[int(words[0])].append(float(words[4])) # print mapping[int(words[0])][0] # }}} for e in tree: event=[] if args.nEvents < evt: print "now should be a break" break trackbranch=ROOT.TBranch() trackfitsta=ROOT.TBranch() if args.pre: nTracks = e.TrackPrePostFit.GetEntriesFast() trackbranch=e.TrackPrePostFit trackfitstat=e.PreTrackFitStat_0 nCluster = e.TpcPreCluster.GetEntriesFast() else: nTracks = e.TrackPostFit.GetEntriesFast() trackbranch=e.TrackPostFit trackfitstat=e.TrackFitStat_0 nCluster = e.TpcCluster.GetEntriesFast() numpartracks=0 # {{{ start loop over tracks if args.debug: print "There are",nTracks," tracks in this event" print "A total of",nCluster," Cluster are found in this event" trackindex=-1 for tr in trackbranch: if type(tr)==ROOT.TpcCluster: print "now tr is a TcpCluster" continue endcap = 0 length=0 poscounter=0 meanz=0 trackindex+=1 theta=tr.getMom().Theta()*ROOT.TMath.RadToDeg() phi = tr.getMom().Phi()*ROOT.TMath.RadToDeg() cand = tr.getCand() tfs=None tfs = trackfitstat.At(trackindex) if type(tfs)!=ROOT.TrackFitStat: print "tfs is not a track fit stat object, its a:",type(tfs)," at event:",evt continue trackrep=tr.getTrackRep(0) fitflag=trackrep.getStatusFlag() resX=tfs.GetResX() resY=tfs.GetResY() resZ=tfs.GetResZ() if args.detid==-1: args.detid=cand.getDetIDs()[0] hitIds = cand.getHitIDs(args.detid) if args.debug: print "I am using detector ID:",args.detid # {{{ start loop over hitIds/cluster hits=[] hitindex=-1 clsumamp=0 clmeanamp=0 digsumamp=0 digmeanamp=0 digsum=0 sampleclsumamp=0 sampleclmeanamp=0 sampleclsum=0 clustersizes=[] clusterdists=[] if args.debug: print "There are ",len(hitIds),"cluster in this track" lastpos=ROOT.TVector3(0,0,0) for hitId in hitIds : hitindex+=1 if args.pre : hit=e.TpcPreSPHit.At(hitId) clid=hit.getClusterID() cl=e.TpcPreCluster.At(clid) else: hit=e.TpcSPHit.At(hitId) clid=hit.getClusterID() cl=e.TpcCluster.At(clid) if type(cl)==ROOT.TpcCluster: pos=cl.pos() else: print "cluster is not a cluster" continue maxSampleInCluster={'amp':0,'sourceId':0,'adcId':0,'chipId':0,'channelId':0} # {{{ start loop overt digis digis=[] samInCluster=0 digInCluster=0 meansigma=0 digperpad={} digimap=None digimap=cl.getDigiMap() sampleclsumamp=0 sampleclmeanamp=0 sampleclsum=0 for dig in digimap: digi=e.TpcDigi.At(dig.first) if digi.amp()>args.digithr: digInCluster+=1 maxSampleInDigi={'amp':0,'sourceId':0,'adcId':0,'chipId':0,'channelId':0} # {{{loop over samples samplemap=None if type(digi)!=ROOT.TpcCluster: samplemap=digi.getSampleMap() else: print "now digi is a TpcCluster" continue samInDigi=0 for sam in samplemap : sampleclsum+=1 if sam.first<65535: sample = e.TpcSample.At(sam.first) else: if args.debug: print "Bad sample number in tree",sam.first if sample.amp()>args.samthr: samInDigi+=1 if sample.amp() > maxSampleInDigi['amp']: maxSampleInDigi['amp'] = sample.amp() maxSampleInDigi['sourceId'] = sample.sourceId() maxSampleInDigi['adcId'] = sample.adcId() maxSampleInDigi['chipId'] = sample.chipId() maxSampleInDigi['channelId'] = sample.channelId() # }}} if maxSampleInDigi['amp']>maxSampleInCluster['amp']: maxSampleInCluster=maxSampleInDigi if not args.sim: meansigma+=sigmas[maxSampleInDigi['sourceId']-802][maxSampleInDigi['adcId']][maxSampleInDigi['chipId']][maxSampleInDigi['channelId']] samInCluster+=samInDigi digperpad[str(digi.padId())]=digperpad.get(str(digi.padId()),0)+1 digsumamp+=digi.amp() digsum+=1 digis.append({'amp':digi.amp(),'time':digi.t(),'pad':digi.padId(), 'maxsample':maxSampleInDigi,'samInDigi':samInDigi}) # end loop over digis # }}} if maxSampleInCluster['amp']>1700 and args.sat==1: continue if hitindex>0: dist=(lastpos-pos).Mag() if poscounter>0 and dist<5: clusterdists.append(dist) length+=dist lastpos=pos poscounter+=1 if hitIds.size() == resX.size(): res = ROOT.TVector3(resX[hitindex],resY[hitindex],resZ[hitindex]) else: res = ROOT.TVector3(9999,9999,9999) res*=10000 r=(math.sqrt(pos[0]*pos[0]+pos[1]*pos[1])) clustersizes.append(cl.size()) # {{{ calculate snr snr=0 if maxSampleInCluster['amp']>0: if not args.sim: somesigma = sigmas[maxSampleInCluster['sourceId']-802][maxSampleInCluster['adcId']][maxSampleInCluster['chipId']][maxSampleInCluster['channelId']] if somesigma > 0: snr = maxSampleInCluster['amp']/somesigma else: print "no noise info for:",maxSampleInCluster['sourceId']-802,maxSampleInCluster['adcId'],maxSampleInCluster['chipId'],maxSampleInCluster['channelId'] snr = 0 else : snr = (maxSampleInCluster['amp']*4)/20 if args.sim: meansigma=4./20. else: if len(digis)>0: meansigma/=len(digis) if cl.get2DSize()>0 and meansigma>0: meansnr=(cl.amp()/meansigma)/cl.get2DSize() else: meansnr=0 meanz+=pos.Z() # }}} # if args.sat==0 or (args.sat==1 and maxSampleInCluster['amp']<1700): hits.append({'rad':r,'pos':pos,'res':res,'size':cl.size(),'size2D':cl.get2DSize(), 'digis':digis,'maxSample':maxSampleInCluster,'snr':snr, 'digInCluster':digInCluster,'samInCluster':samInCluster, 'digPerPad':digperpad,'meansnr':meansnr,'amp':cl.amp()}) clsumamp+=cl.amp() del digperpad #end loop over hitIds/cluster # }}} if len(hits)>0: clmeanamp=clsumamp/len(hits) else: clmeanamp=-1 if digsum>0: digmeanamp=digsumamp/digsum else: digmeanamp=-1 clustersizes.sort() if len(clustersizes)>0: median=clustersizes[int(round(len(clustersizes)/2))] else: median=-1 # print clustersizes # print "median:",median if len(hits)>0: meanz/=len(hits) else: meanz=99 chi2cdf=ROOT.Math.chisquared_cdf(tr.getChiSqu(),tr.getNDF()) event.append({'hits':hits,'curv':2,'theta':theta, 'phi':phi,'mom':tr.getMom().Mag(),'length':length, 'fitflag':fitflag,'chi2':tr.getChiSqu(),'ndf':tr.getNDF(), 'clmeanamp':clmeanamp,'clsumamp':clsumamp,'digmeanamp':digmeanamp, 'digsum':digsum,'digsumamp':digsumamp,'samsum':sampleclsum,'clmedian':median, 'chi2cdf':chi2cdf,'clusterdists':clusterdists,'meanz':meanz}) # end loop over tracks # }}} # {{{ now we fill some histos histsdict['tracknum']['no cuts'].Fill(len(event)) # {{{ start second loop over event noiseev=0 trackcounter=0 if args.single: for tr in range(len(event)): if event[tr]['clmeanamp']>400 and (event[tr]['theta']<10 or event[tr]['theta']>170): # for cl in event[tr]['hits']: # if cl['res'][0]>5000: print "eventnumber (noise):",evt noiseev=1 if noiseev==0: print "eventnumber (old):",evt oldevent=event histsdict['hxy'].Reset() histsdict['hyz'].Reset() if noiseev==1: plotfullevent(oldevent) raw_input("now the actual event") plotfullevent(event) for tr in range(len(event)): if args.single: continue histsdict['cuts'].Fill(0.) fillcutred(1,event[tr],histsdict['cut reduction']) if event[tr]['fitflag']==1 and args.fcut==1: # cuts['fcut']=1+fcut histsdict['theta']['fit flag cut'].Fill(event[tr]['theta']) histsdict['phi']['fit flag cut'].Fill(event[tr]['phi']) histsdict['theta vs clmeanamp f'].Fill(event[tr]['clmeanamp'],event[tr]['theta']) histsdict['cuts'].Fill(1.) continue if event[tr]['ndf'] > 0: histsdict['trackchi2']['chindf no cuts'].Fill(event[tr]['chi2']/event[tr]['ndf']) histsdict['hclmedian'].Fill(event[tr]['clmedian']) histsdict['trackchi2']['chi no cuts'].Fill(event[tr]['chi2']) histsdict['trackchi2']['ndf no cuts'].Fill(event[tr]['chi2']) histsdict['tracklen']['no cuts'].Fill(event[tr]['length']) histsdict['theta']['no cuts'].Fill(event[tr]['theta']) histsdict['phi']['no cut'].Fill(event[tr]['phi']) histsdict['mom']['no cut'].Fill(event[tr]['mom']) histsdict['hmom'].Fill(event[tr]['mom']) histsdict['curvature']['no cut'].Fill(event[tr]['curv']) histsdict['theta vs fitf'].Fill(event[tr]['fitflag'],event[tr]['theta']) histsdict['theta vs clp'].Fill(len(event[tr]['hits']),event[tr]['theta']) histsdict['theta vs mom'].Fill(event[tr]['mom'],event[tr]['theta']) histsdict['theta vs clsumamp'].Fill(event[tr]['clsumamp'],event[tr]['theta']) histsdict['theta vs clmeanamp'].Fill(event[tr]['clmeanamp'],event[tr]['theta']) histsdict['theta vs digsum'].Fill(event[tr]['digsum'],event[tr]['theta']) histsdict['theta vs median'].Fill(event[tr]['clmedian'],event[tr]['theta']) histsdict['phi vs clsumamp'].Fill(event[tr]['clsumamp'],event[tr]['phi']) histsdict['phi vs clmeanamp'].Fill(event[tr]['clmeanamp'],event[tr]['phi']) histsdict['phi vs digsum'].Fill(event[tr]['digsum'],event[tr]['phi']) histsdict['phi vs mom'].Fill(event[tr]['mom'],event[tr]['phi']) histsdict['phi vs clp'].Fill(len(event[tr]['hits']),event[tr]['phi']) histsdict['mom vs clsumamp'].Fill(event[tr]['clsumamp'],event[tr]['mom']) histsdict['mom vs clmeanamp'].Fill(event[tr]['clmeanamp'],event[tr]['mom']) histsdict['mom vs digsum'].Fill(event[tr]['digsum'],event[tr]['mom']) histsdict['clmeanamp'].Fill(event[tr]['clmeanamp']) histsdict['clsumamp'].Fill(event[tr]['clsumamp']) histsdict['digmeanamp'].Fill(event[tr]['digmeanamp']) histsdict['digsumamp'].Fill(event[tr]['digsumamp']) histsdict['digsum'].Fill(event[tr]['digsum']) histsdict['theta vs phi'].Fill(event[tr]['phi'],event[tr]['theta']) histsdict['hclusterpertrack no cut'].Fill(len(event[tr]['hits'])) # if (len(event[tr]['hits'])>0): # histsdict['mean digis in cluster'].Fill(float(event[tr]['digsum'])/float(len(event[tr]['hits']))) # histsdict['mean samples in cluster'].Fill(float(event[tr]['samsum'])/float(len(event[tr]['hits']))) histsdict['hcl medians meana'].Fill(event[tr]['clmedian'],event[tr]['clmeanamp']) histsdict['chi2cdf'].Fill(event[tr]['chi2cdf']) # {{{ all the cuts on tracks are here histsdict['cuts'].Fill(0.) if nTracks > 1: cuts['ntrack']=1 + args.cosmic histsdict['phi']['track cut'].Fill(event[tr]['phi']) histsdict['mom']['track cut'].Fill(event[tr]['mom']) histsdict['theta']['track cut'].Fill(event[tr]['theta']) histsdict['cuts'].Fill(2.) if args.debug==1: print "number of tracks:",nTracks if event[tr]['mom'] < 0.200: histsdict['phi']['mom cut'].Fill(event[tr]['phi']) histsdict['mom']['mom cut'].Fill(event[tr]['mom']) histsdict['theta']['mom cut'].Fill(event[tr]['theta']) cuts['mom']=1 + args.momc histsdict['cuts'].Fill(4.) if (event[tr]['theta'] < 40 or event[tr]['theta'] > 140): cuts['theta']=1 + args.thec histsdict['theta']['theta cut'].Fill(event[tr]['theta']) histsdict['phi']['theta cut'].Fill(event[tr]['phi']) histsdict['mom']['theta cut'].Fill(event[tr]['mom']) histsdict['cuts'].Fill(5.) if (abs(event[tr]['phi']) < 70 or abs(event[tr]['phi']) > 110): cuts['phi']=1 + args.phic histsdict['theta']['phi cut'].Fill(event[tr]['theta']) histsdict['phi']['phi cut'].Fill(event[tr]['phi']) histsdict['mom']['phi cut'].Fill(event[tr]['mom']) histsdict['cuts'].Fill(6.) # if event[tr]['curv'] > 0.05: # histsdict['theta']['curvature cut'].Fill(event[tr]['theta']) # histsdict['phi']['curvature cut'].Fill(event[tr]['phi']) # histsdict['mom']['curvature cut'].Fill(event[tr]['mom']) # histsdict['tracklen']['curvature cut'].Fill(event[tr]['length']) # cuts['curv']=1 + args.curv # histsdict['cuts'].Fill(7.) if len(event[tr]['hits'])<10:# or len(event[tr]['hits'])>35 : histsdict['theta']['cp cut'].Fill(event[tr]['theta']) histsdict['phi']['cp cut'].Fill(event[tr]['phi']) histsdict['mom']['cp cut'].Fill(event[tr]['mom']) histsdict['cuts'].Fill(8.) cuts['cp']=1+args.cpcut if event[tr]['clmeanamp']<400: cuts['clmean']=1+args.clmean histsdict['theta']['cl mean cut'].Fill(event[tr]['theta']) histsdict['phi']['cl mean cut'].Fill(event[tr]['phi']) histsdict['mom']['cl mean cut'].Fill(event[tr]['mom']) histsdict['cuts'].Fill(9.) if event[tr]['clsumamp']<4000: histsdict['theta']['cl sum cut'].Fill(event[tr]['theta']) histsdict['phi']['cl sum cut'].Fill(event[tr]['phi']) histsdict['mom']['cl sum cut'].Fill(event[tr]['mom']) histsdict['cuts'].Fill(10.) if event[tr]['chi2cdf']<0.02: cuts['chi2cdf']=1+args.chi2cdfc histsdict['chi2cdf cutted'].Fill(event[tr]['chi2cdf']) # for clust in event[tr]['hits']: # if clust['snr']<10: # cuts['snr']=1 # }}} fillcutred(2,event[tr],histsdict['cut reduction']) if cuts['clmean']==0: fillcutred(3,event[tr],histsdict['cut reduction']) if cuts['theta']==0: fillcutred(4,event[tr],histsdict['cut reduction']) if cuts['phi']==0: fillcutred(5,event[tr],histsdict['cut reduction']) cutcounter=1 for cut in cuts: cutcounter2=1 if cuts[cut]==2: histsdict['tcuts'].Fill(cutcounter,0,1) for cutt in cuts: if cuts[cutt]==2 and cuts[cut]==2: histsdict['tcuts'].Fill(cutcounter,cutcounter2,1) cutcounter2+=1 # print "cuthist:",cut,cuts[cut],cutcounter cutcounter+=1 tmpcut=cuts['phi'] cuts['phi']=0 if cuts.values().count(2)==0: histsdict['phi']['no phi cut'].Fill(event[tr]['phi']) cuts['phi']=tmpcut tmpcut=cuts['theta'] cuts['theta']=0 if cuts.values().count(2)==0: histsdict['theta']['no theta cut'].Fill(event[tr]['theta']) cuts['theta']=tmpcut tmpcut=cuts['mom'] cuts['mom']=0 if cuts.values().count(2)==0: histsdict['mom']['no mom cut'].Fill(event[tr]['mom']) cuts['mom']=tmpcut if cuts.values().count(2)==0: numpartracks+=1 if event[tr]['ndf'] > 0: histsdict['trackchi2']['chindf all cuts'].Fill(event[tr]['chi2']/event[tr]['ndf']) histsdict['theta']['all cuts'].Fill(event[tr]['theta']) histsdict['phi']['all cuts'].Fill(event[tr]['phi']) histsdict['mom']['all cuts'].Fill(event[tr]['mom']) histsdict['tracklen']['all cuts'].Fill(event[tr]['length']) histsdict['hclusterpertrack'].Fill(len(event[tr]['hits'])) if event[tr]['length']>0: histsdict['hclusterpercm'].Fill(len(event[tr]['hits'])/event[tr]['length']) histsdict['hdigipercm cutted'].Fill(event[tr]['digsum']/event[tr]['length']) histsdict['hdigipercm cutted vs z'].Fill(event[tr]['meanz'],event[tr]['digsum']/event[tr]['length']) histsdict['curvature']['all cuts'].Fill(event[tr]['curv']) histsdict['theta vs phi cut'].Fill(event[tr]['phi'],event[tr]['theta']) histsdict['trackchi2']['chi all cuts'].Fill(event[tr]['chi2']) histsdict['trackchi2']['ndf all cuts'].Fill(event[tr]['chi2']) # histsdict['mean digis in cluster cut'].Fill(float(event[tr]['digsum'])/float(len(event[tr]['hits']))) # histsdict['mean samples in cluster cut'].Fill(float(event[tr]['samsum'])/float(len(event[tr]['hits']))) histsdict['theta vs median cut'].Fill(event[tr]['clmedian'],event[tr]['theta']) histsdict['hclmedian cut'].Fill(event[tr]['clmedian']) histsdict['chi2cdf cut'].Fill(event[tr]['chi2cdf']) for d in event[tr]['clusterdists']: histsdict['clusterdists'].Fill(d) # {{{ loop over the cluster remembersnr=0 clustercounter=-1 if args.debug==1: print "number of cluster:", len(event[tr]['hits']) for clust in event[tr]['hits']: clustercounter+=1 histsdict['cuts'].Fill(20.) amp=clust['maxSample']['amp'] # if cuts.values().count(2)==0: # histsdict['clustermax']['track cuts'].Fill(amp) histsdict['clsizeVsRad']['no cuts'].Fill(clust['rad'],clust['size']) histsdict['rad'].Fill(clust['rad']) histsdict['z vs theta'].Fill(clust['pos'][2],event[tr]['theta']) histsdict['z vs phi'].Fill(clust['pos'][2],event[tr]['phi']) histsdict['clustermax']['all'].Fill(amp) histsdict['clusteramp']['all'].Fill(clust['amp']) histsdict['hcl size vs amp'].Fill(clust['size'],clust['amp']) histsdict['hdigpercl'].Fill(len(clust['digis'])) histsdict['hcl pos z'].Fill(clust['pos'][2]) # for dig in clust['digis']: # if dig['amp'] > 1700: cuts['sat']=0 cuts['rad']=0 cuts['clamp']=0 if amp >1700: cuts['sat']=1+args.sat histsdict['cuts'].Fill(11.) histsdict['clustermax']['sat cut'].Fill(amp) histsdict['clusteramp']['sat cut'].Fill(clust['amp']) histsdict['cuts'].Fill(21.) if clust['res'][0]==99990000.0: cuts['res']=2 if amp0: histsdict['clsizeVsRad']['in curvature cut'].Fill(clust['rad']) histsdict['hxyc'].Fill(clust['pos'][0],clust['pos'][1]) histsdict['hyzc'].Fill(clust['pos'][2],clust['pos'][1]) if cuts['ntrack']>0: histsdict['clusteramp']['track cut'].Fill(clust['amp']) if cuts['mom']>0: histsdict['clusteramp']['mom cut'].Fill(clust['amp']) if cuts.values().count(2) > 0: histsdict['clustermax']['all cutted'].Fill(amp) histsdict['clusteramp']['all cutted'].Fill(clust['amp']) tmpcut=cuts['phi'] cuts['phi']=0 if cuts.values().count(2)==0: histsdict['hxy resR'].Fill(clust['pos'][0],clust['pos'][1], math.sqrt(clust['res'][0]*clust['res'][0]+ clust['res'][1]*clust['res'][1])) if tmpcut==2: histsdict['hxy resY'].Fill(clust['pos'][0],clust['pos'][1],clust['res'][1]) cuts['phi']=tmpcut #here the cuts are applied! if cuts.values().count(2)==0: histsdict['hclsizevsradius'].Fill(clust['rad'],clust['size']) histsdict['clustermax']['all cuts'].Fill(amp) histsdict['clusteramp']['all cuts'].Fill(clust['amp']) histsdict['hdigpercl cut'].Fill(len(clust['digis'])) histsdict['hcl pos z cut'].Fill(clust['pos'][2]) histsdict['hsnr'].Fill(clust['snr']) histsdict['csnr'].Fill(clust['meansnr']) histsdict['z vs theta c'].Fill(clust['pos'][2],event[tr]['theta']) histsdict['z vs phi c'].Fill(clust['pos'][2],event[tr]['phi']) if event[tr]['theta']<90: histsdict['hxy le90'].Fill(clust['pos'][0],clust['pos'][1]) histsdict['hyz le90'].Fill(clust['pos'][2],clust['pos'][1]) else: histsdict['hxy gt90'].Fill(clust['pos'][0],clust['pos'][1]) histsdict['hyz gt90'].Fill(clust['pos'][2],clust['pos'][1]) if clust['snr']<10: histsdict['hxysnr'].Fill(clust['pos'][0],clust['pos'][1]) histsdict['hyzsnr'].Fill(clust['pos'][2],clust['pos'][1]) histsdict['clustermax']['low snr'].Fill(amp) histsdict['clusteramp']['low snr'].Fill(clust['amp']) remembersnr=1 histsdict['clustersize']['3D'].Fill(clust['size']) histsdict['clustersize']['2D'].Fill(clust['size2D']) histsdict['hxy'].Fill(clust['pos'][0],clust['pos'][1]) histsdict['hyz'].Fill(clust['pos'][2],clust['pos'][1]) histsdict['clsizeVsRad']['all cuts'].Fill(clust['rad'],clust['size']) histsdict['hxy res'].Fill(clust['pos'][0],clust['pos'][1],clust['res'][0]) histsdict['hyz res'].Fill(clust['pos'][2],clust['pos'][1],clust['res'][0]) histsdict['hxy amp cut'].Fill(clust['pos'][0],clust['pos'][1],clust['amp']) if not(args.fila and (clustercounter==0 or clustercounter==len(event[tr]['hits'])-1)): histsdict['hresvsclsize'].Fill(clust['res'][0],clust['size']) histsdict['hresvsz'].Fill(clust['pos'][2],clust['res'][0]) histsdict['resxy'].Fill(clust['pos'][0],clust['pos'][1],abs(clust['res'][0])) histsdict['hresX'].Fill(clust['res'][0]) histsdict['hresY'].Fill(clust['res'][1]) histsdict['hresZ'].Fill(clust['res'][2]) #if abs(clust['res'][0])>2000: # print "bad event:",evt # print "clustercounter:",clustercounter,"length:",len(event[tr]['hits'])-1 if clustercounter==0: histsdict['hresX first'].Fill(clust['res'][0]) histsdict['hxy first'].Fill(clust['pos'][0],clust['pos'][1]) if clustercounter==len(event[tr]['hits'])-1: histsdict['hresX last'].Fill(clust['res'][0]) histsdict['hxy last'].Fill(clust['pos'][0],clust['pos'][1]) histsdict['hits per track vs res'].Fill(clust['res'][0],len(event[tr]['hits'])) #hits3d.Fill(clust['pos'][0],clust['pos'][1],clust['pos'][2]) # {{{ loop over z-slices for j in range(len(zSlices)): if clust['pos'][2]0: histsdict['theta']['low snr'].Fill(event[tr]['theta']) histsdict['phi']['low snr'].Fill(event[tr]['phi']) histsdict['tracklen']['low snr'].Fill(event[tr]['length']) if args.debug: print "number of good tracks:",numpartracks event[tr].clear() for cut in cuts: cuts[cut]=0 # }}} # }}} del event histsdict['tracknum']['all cuts'].Fill(numpartracks) if numpartracks>0: goodev+=numpartracks done=int(float(evt)/float(args.nEvents)*100) if done!=olddone: print done,"% done, Good event:",goodev,"(",totgoodev,")" olddone=done evt = evt+1 for cut in cuts: cuts[cut]=0 #end loop over event # }}} totgoodev+=goodev totev +=evt if type(Rfile)==ROOT.TFile: Rfile.Close() tree=None pos=-1 if not args.sim: histsdict['someinfo'].Fill(0.) for i in range(len(sigmas)): for j in range(len(sigmas[i])): for k in range(len(sigmas[i][j])): for l in range(len(sigmas[i][j][k])): pos+=1 if sigmas[i][j][k][l] != None: histsdict['pedsigma'].Fill(pos,sigmas[i][j][k][l]) histsdict['pedmean'].Fill(pos,means[i][j][k][l]) else: histsdict['pedsigma'].Fill(pos,-100) histsdict['pedmean'].Fill(pos,-100) # }}} # {{{ drawing saving/printing the histos if args.outfname == "" : outfolder = "anafiles/" args.outfname = outfolder + "CosmicsAna_" if filecount > 1: for i in range(len(runList)): args.outfname += str(runList[i]) if i!=len(runList)-1: args.outfname += "_" else: args.outfname += args.runs args.outfname += ".root" if args.outfname.find(".root")==-1: args.outfname+=".root" histsdict['theta']['r cl sum']=Divideh(histsdict['theta']['cl sum cut'],histsdict['theta']['no cuts'],"hrtheta1","theta r clsum") histsdict['theta']['r cl mean']=Divideh(histsdict['theta']['cl mean cut'],histsdict['theta']['no cuts'],"hrtheta2","theta r clmean") histsdict['theta']['r cl phi']=Divideh(histsdict['theta']['phi cut'],histsdict['theta']['no cuts'],"hrtheta3","theta r phi") histsdict['theta']['r cl mom']=Divideh(histsdict['theta']['mom cut'],histsdict['theta']['no cuts'],"hrtheta4","theta r mom") histsdict['phi']['r cl sum']=Divideh(histsdict['phi']['cl sum cut'],histsdict['phi']['no cut'],"hrphi1","phi r clsum") histsdict['phi']['r cl mean']=Divideh(histsdict['phi']['cl mean cut'],histsdict['phi']['no cut'],"hrphi2","phi r clmean") histsdict['phi']['r theta']=Divideh(histsdict['phi']['theta cut'],histsdict['phi']['no cut'],"hrphi3","phi r theta") histsdict['phi']['r mom']=Divideh(histsdict['phi']['mom cut'],histsdict['phi']['no cut'],"hrphi4","phi r mom") histsdict['mom']['r cl sum']=Divideh(histsdict['mom']['cl sum cut'],histsdict['mom']['no cut'],"hrmom1","mom r clsum") histsdict['mom']['r cl mean']=Divideh(histsdict['mom']['cl mean cut'],histsdict['mom']['no cut'],"hrmom2","mom r clmean") histsdict['mom']['r phi']=Divideh(histsdict['mom']['phi cut'],histsdict['mom']['no cut'],"hrmom3","mom r phi") histsdict['mom']['r mom']=Divideh(histsdict['mom']['mom cut'],histsdict['mom']['no cut'],"hrmom4","mom r mom") histsdict['mom']['r theta']=Divideh(histsdict['mom']['theta cut'],histsdict['mom']['no cut'],"hrmom5","mom r theta") outfile = ROOT.TFile(args.outfname,"recreate") outfile.cd() canvascounter=-1 histcounter = -1 canvases=[] legends=[] histsdict['hxy res norm']=Divideh(histsdict['hxy res'],histsdict['hxy'],"hxyresnorm","xy residual normalized") histsopt['hxy res norm']="colz" histsdict['hyz res norm']=Divideh(histsdict['hyz res'],histsdict['hyz'],"hyzresnorm","yz residual normalized") histsopt['hyz res norm']="colz" histsdict['hxy amp norm']=Divideh(histsdict['hxy amp cut'],histsdict['hxy'],"hxyampnorm","xy amplitude normalized") histsopt['hxy amp norm']="colz" histsdict['hxy resR norm']=Divideh(histsdict['hxy resR'],histsdict['hxy'],"hxyresRnorm","xy R residual normalized") histsopt['hxy resR norm']="colz" histsdict['hxy resY norm']=Divideh(histsdict['hxy resY'],histsdict['hxy'],"hxyresYnorm","xy Y residual normalized") histsopt['hxy resY norm']="colz" #hist=histsdict['hxy res'] #for i in range(hist.GetNbinsX()): # for j in range(histGetnBinsY()): # hist.SetBinContent(i,j,hist.GetBinContent(i,j)/ for hists in histsdict: histcounter+=1 hist=histsdict[hists] if type(hist)==dict: for h in hist: hist[h].Write() else: hist.Write() if args.plot: if histcounter%9==0: canvascounter+=1 canvases.append(ROOT.TCanvas("c"+str(canvascounter),str(canvascounter))) canvases[canvascounter].Divide(3,3) canvases[canvascounter].cd(histcounter%9+1) if type(hist)==dict: Drawh(hist,hists,legends, canvases[canvascounter].GetPad(histcounter%9+1)) else: if canvasopt.get(hists,"")=="logz": canvases[canvascounter].GetPad(histcounter%9+1).SetLogz() hist.GetXaxis().SetTitle(axlabel[0].get(hists,"")) hist.GetYaxis().SetTitle(axlabel[1].get(hists,"")) hist.Draw(histsopt.get(hists,"")) canvases[canvascounter].Update() if args.printc: if histcounter==0: tempcanv=ROOT.TCanvas("tc",hists) tempcanv.cd() if type(hist)==dict: Drawh(hist,hists,legends,tempcanv) else: hist.GetXaxis().SetTitle(axlabel[0].get(hists,"")) hist.GetYaxis().SetTitle(axlabel[1].get(hists,"")) hist.Draw(histsopt.get(hists,"")) tempcanv.SaveAs(args.pdir+"/"+hists+".eps") tempcanv.Clear() # {{{ plot all theta histos histcounter = -1 #canvascounter = -1 if args.plot: for h in histsdict['theta']: hist=histsdict['theta'][h] histcounter += 1 if histcounter%9==0: canvascounter += 1 canvases.append(ROOT.TCanvas("c"+str(canvascounter),"theta "+str(canvascounter))) canvases[canvascounter].Divide(3,3) canvases[canvascounter].cd(histcounter%9+1) hist.Draw(histsopt.get(hists,"")) canvases[canvascounter].Update() # }}} # {{{ plot all phi histos histcounter = -1 #canvascounter = -1 if args.plot: for h in histsdict['phi']: hist=histsdict['phi'][h] histcounter += 1 if histcounter%9==0: canvascounter += 1 canvases.append(ROOT.TCanvas("c"+str(canvascounter),"phi "+str(canvascounter))) canvases[canvascounter].Divide(3,3) canvases[canvascounter].cd(histcounter%9+1) hist.Draw(histsopt.get(hists,"")) canvases[canvascounter].Update() # }}} # {{{ plot all mom histos histcounter = -1 #canvascounter = -1 if args.plot: for h in histsdict['mom']: hist=histsdict['mom'][h] histcounter += 1 if histcounter%9==0: canvascounter += 1 canvases.append(ROOT.TCanvas("c"+str(canvascounter),"mom "+str(canvascounter))) canvases[canvascounter].Divide(3,3) canvases[canvascounter].cd(histcounter%9+1) hist.Draw(histsopt.get(hists,"")) canvases[canvascounter].Update() # }}} # {{{ single histo crap if hists=='tcuts': templeg=ROOT.TLegend(0.7,.10,0.9,.9,"","NDC") cutcounter=1 for cut in cuts: templeg.AddEntry(0,str(cutcounter)+" "+cut) cutcounter+=1 templeg.Draw() # }}} scanv=[] sc=-1 canvascounter+=1 hslices['hxy res norm slice']={} for i in range(len(zSlices)): hslices['hxy res norm slice'][i]=Divideh(hslices['hxy res slice'][i],hslices['hxy slice'][i],"hxyresnorm"+str(zSlices[i]),"xy residual normalized "+str(zSlices[i])) saxlabel[0]['hxy res norm slice'+str(i)]="X (cm)" saxlabel[1]['hxy res norm slice'+str(i)]="Y (cm)" hslices['hxy res norm slice'][i].GetZaxis().SetRangeUser(-1000.,1000.) hslices['hxy slice'][i].GetZaxis().SetRangeUser(-1000.,1000.) hslices['hxy res slice'][i].GetZaxis().SetRangeUser(-1000.,1000.) sliceopt['hxy res norm slice']="colz" for slicer in hslices: sc+=1 sdict=hslices[slicer] if args.plot: scanv.append(ROOT.TCanvas("c"+str(canvascounter),slicer)) scanv[sc].Divide(3,2) for i in range(len(zSlices)): sdict[i].GetXaxis().SetTitle(saxlabel[0].get(slicer+str(i),"")) sdict[i].GetYaxis().SetTitle(saxlabel[1].get(slicer+str(i),"")) if args.plot: scanv[sc].cd(i+1) sdict[i].Draw(sliceopt.get(slicer,"")) sdict[i].Write() canvascounter+=1 # }}} print 'Histos are written into file:', args.outfname if not args.batch: if args.plotsingle[0]!="": canvassss=[] for i in args.plotsingle: canvassss.append(ROOT.TCanvas()) if histsdict.get(i,"")!="": histsdict[i].Draw() nix=raw_input("Press enter to finish")