diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..37cb604 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.root +**/.root +*/**/*.root \ No newline at end of file diff --git a/distribsForHNLandBG.py b/distribsForHNLandBG.py index 5be41d0..dcca120 100755 --- a/distribsForHNLandBG.py +++ b/distribsForHNLandBG.py @@ -1,3 +1,4 @@ +import sys import ROOT as r import numpy as np import funcsByBarbara as tools @@ -5,6 +6,7 @@ r.gROOT.SetBatch(r.kTRUE) r.gROOT.ProcessLine(".X lhcbStyle.C") plots = {} +pdg = r.TDatabasePDG.Instance() def addPlot(key, x_axis, tree, string, cuts=None): if not cuts: @@ -84,35 +86,35 @@ ftiny = r.TFile('../DATA/SmallHNLs/ShipAna.root', 'read') tiny = ftiny.Get('ShipAna') -tiny_geo = tools.searchForNodes3_xyz_dict('../DATA/SmallHNLs/geofile_full.10.0.Pythia8-TGeant4-0.1GeV-009-eenu.root') +tiny_geo = tools.searchForNodes3_xyz_dict('../DATA/SmallHNLs/geofile_full.10.0.Pythia8-TGeant4-0.1GeV-011.root') -fbg = r.TFile('../DATA/KLONG/ShipAna.root','read') -bg = fbg.Get('ShipAna') -bg_geo = tools.searchForNodes2_xyz_dict('../DATA/KLONG/geofile_full.10.0.Genie-TGeant4_70cc30nc.root') +#fbg = r.TFile('../DATA/KLONG/ShipAna.root','read') +#bg = fbg.Get('ShipAna') +#bg_geo = tools.searchForNodes2_xyz_dict('../DATA/KLONG/geofile_full.10.0.Genie-TGeant4_70cc30nc.root') fcosmics = r.TFile('../DATA/Cosmics/ShipAna.root','read') cosmics = fcosmics.Get('ShipAna') cosmics_geo = tools.searchForNodes2_xyz_dict('../DATA/Cosmics/geofile_full.Signal.10.0.Cosmics3.25.1-8.root') -fnubar = r.TFile('../DATA/NewKLONG/ShipAna74-75_nubar.root','read') +fnubar = r.TFile('../DATA/NewKLONG/ShipAna_nubar.root','read')#'../DATA/NewKLONG/ShipAna74-75_nubar.root','read') nubar = fnubar.Get('ShipAna') -nubar_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_neutrino74-75.root') +nubar_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nubar.root')#neutrino74-75.root') -fnu = r.TFile('../DATA/NewKLONG/ShipAna77_nu.root','read') +fnu = r.TFile('../DATA/NewKLONG/ShipAna_nu.root','read')#'../DATA/NewKLONG/ShipAna77_nu.root','read') nu = fnu.Get('ShipAna') -nu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_neutrino77.root') +nu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nu.root')#neutrino77.root') # For a correct normalisation, we subtract the events that have an HNL vertex produced outside of the 60m fiducial volume studies = { - 'pimu': {'data':pimu, 'geo':pimu_geo, 'ntot':20000-4175},#1000}, - 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'ntot':20000-4125},#7000}, - 'tiny': {'data': tiny, 'geo':tiny_geo, 'ntot':20000-4452}, - 'bg': {'data':bg, 'geo':bg_geo, 'ntot':142856},#99999}#199998 - 'cosmics': {'data':cosmics, 'geo':cosmics_geo, 'ntot':cosmics.GetEntriesFast()}, - 'nubar': {'data': nubar, 'geo':nubar_geo, 'ntot':99999*2}, - 'nu': {'data': nu, 'geo':nu_geo, 'ntot':99999} + 'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-4175},#1000}, + 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-4125},#7000}, + 'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-4498}, + #'bg': {'data':bg, 'geo':bg_geo, 'ntot':142856},#99999}#199998 + 'cosmics': {'data':cosmics, 'geo':cosmics_geo, 'seen':[], 'ntot':cosmics.GetEntriesFast()}, + 'nubar': {'data': nubar, 'geo':nubar_geo, 'seen':[], 'ntot':99999*2},#99999*2}, + 'nu': {'data': nu, 'geo':nu_geo, 'seen':[], 'ntot':99999*2}#99999} } print 'Found %s entries in the cosmics file'%cosmics.GetEntriesFast() @@ -124,6 +126,95 @@ # - DOCA < 10 cm # - has no veto? + + +class Cuts(object): + def __init__(self, cutConfigList, study): + self.cutList = [] + try: + for (index,cutConfig) in enumerate(cutConfigList): + self.cutList.append( self.Cut( self, index, cutConfig['name'], study, cutConfig['params'] ) ) + except BaseException,e: + print "cutConfigList badly set: %r" %e + def process(self, eventList): + for event in eventList: + for cut in self.cutList: + if not cut.count(event): break + return self + def select(self, event, imax): + for i in xrange(imax): + if not self.cutList[i].selection(event, self.cutList[i].study, self.cutList[i].params): + return False + return True + + class Cut(object): + def __init__(self, cutCollection, index, cutName, study, params): + self.cutCollection = cutCollection + self.index = index + self.name = cutName + self.params = params + self.study = study + self.total = 0 + self.weight = 0 + self.eventNumbers = [] + self.duplicates = [] + self.selection = self.__get(cutName) + def __get(self, cutName): + name = "cut_"+cutName + if hasattr(self,name): + try: + return getattr(self,name) + except AttributeError,e: + return None + else: + raise BaseException ( "bad cut name: %r" %cutName ) + return None + def count(self, event): + try: + if (self.selection(event,self.study,self.params)): + self.total += 1 + if (event.EventNumber not in self.eventNumbers): + self.weight += self.compute_weight(event, self.study) + self.eventNumbers.append(event.EventNumber) + else: + self.duplicates.append(event.EventNumber) + return True + else : + return False + except BaseException,e: + print "Error performing %r: %r" %(self.name,e) + def cut_fit_converged(self, event, study, params): + return fit_converged(event) + def cut_chi2(self, event, study, params): + return chi2_cut(event,params['chi2max']) + def cut_has_muons(self, event, study, params): + for station in params['stations']: + if not has_muon_station(event, study, station): + return False + return True + def cut_doca(self, event, study, params): + return doca_cut(event, params['docamax']) + def cut_z(self, event, study, params): + return z_cut(event, study) + def cut_normalization(self, event, study, params): + return goodTruthVtx(event, study) + def cut_good_tracks(self, event, study, params): + return goodTrack(event) + def cut_ip(self, event, study, params): + return ip_cut(event, params['ipmax']) + def cut_vtx_in_ellipse(self, event, study, params): + return vtxInAcceptance(event, study) + def cut_n_candidates(self, event, study, params): + #print self.cutCollection.cutList[self.index-1].name + if (event.EventNumber in self.cutCollection.cutList[self.index-1].duplicates): + #if self.cutCollection.select(event, self.index): + return nParticles_cut(event, params['nmax']) + return True + def compute_weight(self, event, study): + return weight(event, study) + + + def fit_converged(event): flags = event.DaughtersFitConverged if np.product(flags): @@ -160,8 +251,8 @@ # minz = studies[study]['geo']['lidT1I_1']['z']['pos'] + studies[study]['geo']['lidT1I_1']['z']['dim'] + 500. # maxz = studies[study]['geo']['lidT6I_1']['z']['pos'] - studies[study]['geo']['lidT6I_1']['z']['dim'] #except: - minz = studies[study]['geo']['Veto_5']['z']['pos'] + studies[study]['geo']['Veto_5']['z']['dim'] - maxz = studies[study]['geo']['Tr1_1']['z']['pos'] + studies[study]['geo']['Tr1_1']['z']['dim'] + minz = studies[study]['geo']['Veto_5']['z']['pos'] + studies[study]['geo']['Veto_5']['z']['dim'] + 0.01 + maxz = studies[study]['geo']['Tr1_1']['z']['pos'] - studies[study]['geo']['Tr1_1']['z']['dim'] - 0.01 if minz < z < maxz: return True return False @@ -172,6 +263,11 @@ return True return True +def nParticles_cut(event, nmax): + if event.NParticles > nmax: + return False + return True + # This checks that: # - there is at least one tracking station before and one after the magnet # - hits in these tracking stations are within a 5x10m ellipse @@ -213,19 +309,33 @@ return tools.inEllipse(x,y,500./2,1000./2) def weight(event, study): - w = event.HNLw / studies[study]['ntot'] - if ('bg' in study) or (study == 'nu' or 'study' == 'nubar'): + #if event.EventNumber not in studies[study]['seen']: + #studies[study]['seen'].append(event.EventNumber) + if ('bg' in study) or (study == 'nu' or study == 'nubar'): w = event.NuWeight - if 'cosmics' in study: + elif 'cosmics' in study: w = event.HNLw * 1.e6 / 200. + elif study in ['pimu', 'mumunu', 'tiny']: + w = event.HNLw / studies[study]['ntot'] + else: + print 'WEIGHTING ERROR: ', study + sys.exit() + #else: + # w = 0. + #if study in ['pimu', 'mumunu', 'tiny']: + # print study, event.EventNumber, event.NParticles return w -def computeResidualWeight(filename): +def computeResidualWeight(filename=None): if not filename: filename = "forElena_nu_notVetoed.txt" fScan = file(filename) lines = fScan.readlines() - + sample = nu + s = 'nu' + if ('bar' in filename) or ('anti' in filename): + sample = nubar + s = 'nubar' lines = lines[3:-1] events = [] for line in lines: @@ -234,19 +344,14 @@ count = 0 weighted_count = 0. print "unvetoed: ", len(events) - print "total Reconstructed: ", nu.GetEntriesFast() - #for entry in events: - # bg.GetEntry(entry) - # particles = bg.Particles - # for p in particles: - # count+=1 - # weighted_count += bg.NuWeight - for event in nu: + print "total Reconstructed: ", sample.GetEntriesFast() + for event in sample: if event.EventNumber in events: - print event.EventNumber - if True:#fit_converged(event): - count_reco +=1 - if True:#fit_converged(event): #and doca_cut(event, 50.) and ip_cut(event, 250.) and z_cut(event, 'nu') and vtxInAcceptance(event, 'nu'): + count_reco +=1 + if (fit_converged(event) and chi2_cut(event, 5.) and z_cut(event, s) and vtxInAcceptance(event, s) + and goodTrack(event) and doca_cut(event, 30) and has_muon_station(event, s, 1) + and has_muon_station(event, s, 2) and ip_cut(event, 250.)): + #if True:#fit_converged(event): #and doca_cut(event, 50.) and ip_cut(event, 250.) and z_cut(event, 'nu') and vtxInAcceptance(event, 'nu'): #print event.EventNumber, 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1] print event.EventNumber print 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1] @@ -268,6 +373,80 @@ else: return float_str +cutConfigList = [ {'name': 'normalization', 'params':{}}, + {'name': 'fit_converged', 'params':{}}, + {'name': 'chi2', 'params':{'chi2max':5.}}, + {'name': 'z', 'params':{}}, + {'name': 'vtx_in_ellipse','params':{}}, + {'name': 'good_tracks', 'params':{}}, + {'name': 'doca', 'params':{'docamax':50.}}, + {'name': 'ip', 'params':{'ipmax':500.}}, + {'name': 'has_muons', 'params':{'stations':[1,2]}}, + {'name': 'doca', 'params':{'docamax':30.}}, + {'name': 'ip', 'params':{'ipmax':250.}}, + {'name': 'n_candidates', 'params':{'nmax':1}}, +] + +def numbers_v2(): + mycuts, results = {}, {} + for study in studies: + print 'Reading %s...'%study + mycuts[study] = Cuts(cutConfigList, study).process(studies[study]['data']) + for cut in mycuts[study].cutList: + print '%s\t%s\t%s\t%s'%(cut.name, cut.total, len(cut.eventNumbers), cut.weight) + results[study] = {} + results[study]['numbers']=[ cut.total for cut in mycuts[study].cutList[1:] ] + results[study]['numbersWithWeights']=[ cut.weight for cut in mycuts[study].cutList[1:] ] + print + print 'Events with more than one particle:' + for entry in mycuts[study].cutList[-2].duplicates: #[ x for x in mycuts[study].cutList[-1].eventNumbers if mycuts[study].cutList[-1].eventNumbers.count(x) > 1 ]: + for tr in studies[study]['data']: + if tr.EventNumber == entry: + if mycuts[study].select(tr, len(mycuts[study].cutList)-1): + print tr.EventNumber, tr.vtxz, tr.IP0, [pdg.GetParticle(code).GetName() for code in tr.DaughtersTruthPDG] + print + print + return results + +def makeTex(results): + print + print 'Background:' + print '\t\\begin{table}' + print '\t\\centering Backgrounds: \\tiny' + print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}' + print '\t\t\t\\hline' + print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\' + print '\t\t\t\\hline' + cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm', + 'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$'] + for i in xrange(len(results['nu']['numbers'])): + print '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['nu']['numbersWithWeights'][i]), ' & ', '%s'%results['nu']['numbers'][i], ' & ',\ + '%s'%latex_float(results['nubar']['numbersWithWeights'][i]), ' & ', '%s'%results['nubar']['numbers'][i], ' & ',\ + '%s'%latex_float(results['cosmics']['numbersWithWeights'][i]), ' & ', '%s'%results['cosmics']['numbers'][i], '\\\\ ' + print '\t\t\t\\hline' + print '\t\t\\end{tabular}' + print '\t\\end{table}' + print + print 'Signal:' + print '\t\\begin{table}' + print '\t\\centering Signal: \\tiny' + print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}' + print '\t\t\t\\hline' + print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\' + print '\t\t\t\\hline' + cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm', + 'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$'] + for i in xrange(len(results['nu']['numbers'])): + print '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['pimu']['numbersWithWeights'][i]), ' & ', '%s'%results['pimu']['numbers'][i], ' & ',\ + '%s'%latex_float(results['mumunu']['numbersWithWeights'][i]), ' & ', '%s'%results['mumunu']['numbers'][i], ' & ',\ + '%s'%latex_float(results['tiny']['numbersWithWeights'][i]), ' & ', '%s'%results['tiny']['numbers'][i], '\\\\ ' + print '\t\t\t\\hline' + print '\t\t\\end{tabular}' + print '\t\\end{table}' + print + + + def numbers_for_TP(): results = {} for s in ['pimu', 'mumunu', 'tiny', 'nu', 'nubar', 'cosmics']:#studies: @@ -276,11 +455,14 @@ ntot = t.GetEntriesFast() nAfterCuts = [0]*10 nAfterCuts_weighted = [0.]*10 + initialEventNums = [] + eventNums = [] for event in t: - w = weight(event,s) + w = weight(event,s) # e' importante fare questa cosa a monte, x come e' definita la funzione! if fit_converged(event) and goodTruthVtx(event,s): nAfterCuts[0] +=1 nAfterCuts_weighted[0] += w + initialEventNums.append(int(event.EventNumber)) if chi2_cut(event, 5.): nAfterCuts[1]+=1 nAfterCuts_weighted[1]+=w @@ -313,6 +495,8 @@ if ip_cut(event, 250.): nAfterCuts[8+1] +=1 nAfterCuts_weighted[8+1] +=w + eventNums.append(int(event.EventNumber)) + studies[s]['seen'].append(event.EventNumber) print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot print '%s \t weighted survived particles: '%s, nAfterCuts_weighted, ' of %s total'%ntot @@ -320,13 +504,21 @@ results[s]['ntot']=studies[s]['ntot'] results[s]['numbers']=nAfterCuts results[s]['numbersWithWeights']=nAfterCuts_weighted + results[s]['listOfSurvivors']=list(set(eventNums)) + studies[s]['listOfSurvivors']=list(set(eventNums)) + results[s]['listOfRecoed']=list(set(initialEventNums)) + studies[s]['listOfRecoed']=list(set(initialEventNums)) print for s in ['pimu', 'mumunu', 'tiny']:#studies.keys():#['pimu', 'mumunu', 'bg']: - for list in results[s]: - print s,list,results[s][list] + for lst in results[s]: + if not ('Survivors' in lst) and not ('Recoed' in lst): + print s,list,results[s][lst] + print + print 'Events with surviving nu: %s \t with nubar: %s'%(len(results['nu']['listOfSurvivors']), len(results['nubar']['listOfSurvivors'])) + print print 'Background:' print '\t\\begin{table}' - print '\t\\centering\\footnotesize' + print '\t\\centering Backgrounds: \\tiny' print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}' print '\t\t\t\\hline' print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\' @@ -342,7 +534,7 @@ print '\t\\end{table}' print 'Signal:' print '\t\\begin{table}' - print '\t\\centering\\footnotesize' + print '\t\\centering Signal: \\tiny' print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}' print '\t\t\t\\hline' print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\' @@ -356,8 +548,74 @@ print '\t\t\t\\hline' print '\t\t\\end{tabular}' print '\t\\end{table}' + print + out = {'pimu':[], 'mumunu':[], 'tiny':[]} + for s in ['pimu', 'mumunu', 'tiny']: + for event in studies[s]['data']: + if event.EventNumber in studies[s]['listOfSurvivors']: + if event.NParticles > 1: + out[s].append(event.EventNumber) + #print s, event.EventNumber, event.NParticles, r.gGeoManager.FindNode(event.vtxx, event.vtxy, event.vtxz).GetName() + print s, list(set(out[s])) +def badEvents(): + numbers_for_TP() + hnpnu_recoed = r.TH1F('np nu recoed', 'np nu recoed', 60, 0, 59) + hnpnu_survivors = r.TH1F('np nu survived', 'np nu survived', 60, 0, 59) + seen = [] + for part in studies['nu']['data']: + if not (part.EventNumber in seen): + if part.EventNumber in studies['nu']['listOfSurvivors']: + hnpnu_survivors.Fill(part.NParticles) + if part.EventNumber in studies['nu']['listOfRecoed']: + hnpnu_recoed.Fill(part.NParticles) + seen.append(part.EventNumber) + #c1 = r.TCanvas(); c1.cd() + #hnpnu.Draw() + #c1.Print('badEvents_nu_nparticles.png') + hnpnubar_recoed = r.TH1F('np nubar recoed', 'np nubar recoed', 60, 0, 59) + hnpnubar_survivors = r.TH1F('np nubar survived', 'np nubar survived', 60, 0, 59) + seen = [] + for part in studies['nubar']['data']: + if not (part.EventNumber in seen): + if part.EventNumber in studies['nubar']['listOfSurvivors']: + hnpnubar_survivors.Fill(part.NParticles) + if part.EventNumber in studies['nubar']['listOfRecoed']: + hnpnubar_recoed.Fill(part.NParticles) + seen.append(part.EventNumber) + #c1 = r.TCanvas(); c1.cd() + #hnpnubar.Draw() + #c1.Print('badEvents_nubar_nparticles.png') + hnpnu_recoed.SetLineColor(r.kBlack); hnpnu_recoed.SetMarkerColor(r.kBlack) + hnpnu_survivors.SetLineColor(r.kMagenta); hnpnu_survivors.SetMarkerColor(r.kMagenta) + c1 = r.TCanvas(); c1.cd() + hnpnu_recoed.DrawNormalized('p'); hnpnu_survivors.DrawNormalized('p same') + c1.Print('nuAndNubar/nu_recoed_and_survived.png') + c1.Print('nuAndNubar/nu_recoed_and_survived.pdf') + c1.Close() + hnpnubar_recoed.SetLineColor(r.kBlack); hnpnubar_recoed.SetMarkerColor(r.kBlack) + hnpnubar_survivors.SetLineColor(r.kMagenta); hnpnubar_survivors.SetMarkerColor(r.kMagenta) + c1 = r.TCanvas(); c1.cd() + hnpnubar_recoed.DrawNormalized('p'); hnpnubar_survivors.DrawNormalized('p same') + c1.Print('nuAndNubar/nubar_recoed_and_survived.png') + c1.Print('nuAndNubar/nubar_recoed_and_survived.pdf') + c1.Close() + hnpnu_recoed.SetLineColor(r.kBlue); hnpnu_recoed.SetMarkerColor(r.kBlue) + hnpnubar_recoed.SetLineColor(r.kGreen); hnpnubar_recoed.SetMarkerColor(r.kGreen) + c1 = r.TCanvas(); c1.cd() + hnpnu_recoed.DrawNormalized('p'); hnpnubar_recoed.DrawNormalized('p same') + c1.Print('nuAndNubar/nu_and_nubar_recoed.png') + c1.Print('nuAndNubar/nu_and_nubar_recoed.pdf') + c1.Close() + hnpnu_survivors.SetLineColor(r.kBlue); hnpnu_survivors.SetMarkerColor(r.kBlue) + hnpnubar_survivors.SetLineColor(r.kGreen); hnpnubar_survivors.SetMarkerColor(r.kGreen) + c1 = r.TCanvas(); c1.cd() + hnpnu_survivors.DrawNormalized('p'); hnpnubar_survivors.DrawNormalized('p same') + c1.Print('nuAndNubar/nu_and_nubar_survived.png') + c1.Print('nuAndNubar/nu_and_nubar_survived.pdf') + c1.Close() + def do_my_studies(): diff --git a/lookAtEvents.py b/lookAtEvents.py new file mode 100644 index 0000000..ebf4775 --- /dev/null +++ b/lookAtEvents.py @@ -0,0 +1,112 @@ +import ROOT as r +import funcsByBarbara as tools +cSaver = [] + +def getName(code): + try: + mother = pdg.GetParticle( code ).GetName() + except: + mother = 'beam' + if code == 9900015: + mother = 9900015 + return mother + +def numOfFakeCandidates(anaFile = '../DATA/NewPIMU/ShipAna.root'): + f = r.TFile(anaFile, 'read') + t = f.Get('ShipAna') + numTot = 0 + numFake = 0 + for HNL in t: + mums = HNL.DaughtersTruthMotherPDG + for mum in mums: + if mum != 9900015: + numFake += 1 + break + numTot +=1 + cSaver.append(r.TCanvas()) + t.Draw("vtxy:vtxx","DaughtersTruthMotherPDG != 9900015 && DOCA < 30. && IP0 < 250. && vtxz > -1958. && vtxz < 2578.") + cSaver.append(r.TCanvas()) + t.Draw("vtxz","DaughtersTruthMotherPDG != 9900015 && DOCA < 30. && IP0 < 250. && vtxz > -1958. && vtxz < 2578.") + cSaver.append(r.TCanvas()) + t.Draw("veto5_y:veto5_x","DaughtersTruthMotherPDG != 9900015 && DOCA < 30. && IP0 < 250. && vtxz > -1958. && vtxz < 2578.") + #t.GetEntry(811) + #print int(t.EventNumber) + #print len(t.straw_x) + #print len(t.veto5_x) + #print + #t.GetEntry(1298) + #print int(t.EventNumber) + #print len(t.straw_x) + #print len(t.veto5_x) + #print + #t.GetEntry(2) + #print int(t.EventNumber) + #print len(t.straw_x) + #print len(t.veto5_x) + #print + cutVeto = "!(Alt$(veto5_z,0))" + cSaver.append(r.TCanvas()) + t.Draw("vtxz","DaughtersTruthMotherPDG != 9900015 && DOCA < 30. && IP0 < 250. && vtxz > -1958. && vtxz < 2578. && %s"%cutVeto) + cutLS = "!(Alt$(liquidscint_z,0))" + cSaver.append(r.TCanvas()) + t.Draw("vtxz","DaughtersTruthMotherPDG != 9900015 && DOCA < 30. && IP0 < 250. && vtxz > -1958. && vtxz < 2578. && %s && %s"%(cutVeto, cutLS)) + cutEllipse = "( (vtxx*vtxx)/(250.*250.) + (vtxy*vtxy)/(500.*500.) ) < 1." + cutVtx = "!(BadTruthVtx)" + badMum = "DaughtersTruthMotherPDG[][0] != 9900015 && DaughtersTruthMotherPDG[][1] != 9900015" + goodMum = "DaughtersTruthMotherPDG[][0] == 9900015 && DaughtersTruthMotherPDG[][1] == 9900015" + docaVtxIp = "DOCA < 30. && IP0 < 250. && vtxz > -1958. && vtxz < 2578." + trackQ = "(DaughtersChi2/DaughtersNPoints)<5. && DaughtersAlwaysIn" + #t.Scan("EventNumber","%s && %s && %s && %s && %s && %s && %s"%(badMum, docaVtxIp, cutVeto, cutLS, cutEllipse, cutVtx, trackQ)) + #t.Scan("EventNumber","%s && %s && %s && %s && %s && %s && %s"%(goodMum, docaVtxIp, cutVeto, cutLS, cutEllipse, cutVtx, trackQ)) + t.Draw(">>elistAll","%s && %s && %s && %s && %s && %s"%(docaVtxIp, cutVeto, cutLS, cutEllipse, cutVtx, trackQ)) + elistAll = r.gDirectory.Get("elistAll") + print "Total candidates after the cuts: ", elistAll.GetN() + t.Draw(">>elistBad","%s && %s && %s && %s && %s && %s && %s"%(badMum, docaVtxIp, cutVeto, cutLS, cutEllipse, cutVtx, trackQ)) + elistBad = r.gDirectory.Get("elistBad") + print "Bad candidates after the cuts: ", elistBad.GetN() + t.Draw(">>elistGood","%s && %s && %s && %s && %s && %s && %s"%(goodMum, docaVtxIp, cutVeto, cutLS, cutEllipse, cutVtx, trackQ)) + elistGood = r.gDirectory.Get("elistGood") + print "Good candidates after the cuts: ", elistGood.GetN() + return numFake, numTot + + +def lookAtSingleEvent(evNumber=5771, inputfile='../DATA/NewPIMU/ship.10.0.Pythia8-TGeant4-1.0GeV-008_rec.root'): + #inputfile = '../DATA/NewPIMU/ship.10.0.Pythia8-TGeant4-1.0GeV-008_rec.root' + f = r.TFile(inputfile,'read') + t = f.Get('cbmsim') + tGeom = r.TGeoManager("Geometry", "Geane geometry") + gMan = tGeom.Import(inputfile.replace('ship','geofile_full').replace('_rec','')) + fGeo = r.gGeoManager + pdg = r.TDatabasePDG().Instance() + t.GetEntry(evNumber) + nCandidates = 0 + for p in t.Particles: + nCandidates += 1 + print + print '\t %s candidates in this event.'%nCandidates + print + print 'Tracks in this event: name, fStartZ, node, mother' + print + for tr in t.MCTrack: + name = getName(tr.GetPdgCode()) + try: + mother = getName(t.MCTrack[tr.GetMotherId()].GetPdgCode()) + except: + mother = 'beam' + print '%s \t %.2f \t %s \t %s'%(name, tr.GetStartZ(), fGeo.FindNode(tr.GetStartX(),tr.GetStartY(),tr.GetStartZ()).GetName(), mother) + print + print 'Candidates in this event: daugther name, daughter node, daugther mother (x2)' + print + for p in t.Particles: + t1,t2 = p.GetDaughter(0),p.GetDaughter(1) + pdg1, mumPdg1, truthX1, truthY1, truthZ1 = tools.retrieveMCParticleInfo(t, t1) + pdg2, mumPdg2, truthX2, truthY2, truthZ2 = tools.retrieveMCParticleInfo(t, t1) + name1 = getName(pdg1) + name2 = getName(pdg2) + mother1 = getName(mumPdg1) + mother2 = getName(mumPdg2) + node1 = fGeo.FindNode(truthX1, truthY1, truthZ1).GetName() + node2 = fGeo.FindNode(truthX2, truthY2, truthZ2).GetName() + print '%s\t%s\t%s\t%s\t\t%s\t%s\t%s\t%s'%( truthZ1, name1, node1, mother1, truthZ2, name2, node2, mother2 ) + + diff --git a/lookAtGeo.py b/lookAtGeo.py index 6f51cc6..b9e3de7 100755 --- a/lookAtGeo.py +++ b/lookAtGeo.py @@ -1,137 +1,154 @@ -import ROOT,os,sys,getopt,time -import shipunit as u -import shipRoot_conf -from ShipGeoConfig import ConfigRegistry -from ROOT import * - -# function to prepare two dictionaries: one that has the node names as key and the related intex -# in the nodes array as value and the other one with key and value swapped. -def prepareNodeDictionaries(nodes): - # key: index, value: name of the node - nodeDict_index = {} - # key: name of the node, value: index - nodeDict_name = {} - for (i,n) in enumerate(nodes): - #print i,n.GetName() - nodeDict_index[i]=n.GetName() - nodeDict_name[n.GetName()]=i - return {'nodeDict_index':nodeDict_index, 'nodeDict_name':nodeDict_name} - -# function to print the node names and their index. Usefull to search the name of the element you would like -# to study if you don't know it by heart. Then the name can be used with the other functions for the specific studies. -def searchForNodes(inputFile): - r = loadGeometry(inputFile) - fGeo = r['fGeo'] - ## Get the top volume - fGeo = ROOT.gGeoManager - tv = fGeo.GetTopVolume() - nodes = tv.GetNodes() - for (i,n) in enumerate(nodes): - print i,n.GetName() - -def searchForNodes2(inputFile): - r = loadGeometry(inputFile) - fGeo = r['fGeo'] - ## Get the top volume - #fGeo = ROOT.gGeoManager - tv = fGeo.GetTopVolume() - nodes = tv.GetNodes() - for (i,n) in enumerate(nodes): - print i,n.GetName(),findPositionElement(n)['z'],findDimentionBoxElement(n)['z'] - -def searchForNodes2_xyz(inputFile): - r = loadGeometry(inputFile) - fGeo = r['fGeo'] - ## Get the top volume - #fGeo = ROOT.gGeoManager - tv = fGeo.GetTopVolume() - nodes = tv.GetNodes() - for (i,n) in enumerate(nodes): - print i,n.GetName() - print " x: ", findPositionElement(n)['x'],findDimentionBoxElement(n)['x'] - print " y: ", findPositionElement(n)['y'],findDimentionBoxElement(n)['y'] - print " z: ", findPositionElement(n)['z'],findDimentionBoxElement(n)['z'] - -# basic function to be called to load the geometry from a file -def loadGeometry(inputFile): - try: - dy = float( inputFile.split("/")[-1].split(".")[1]) - except: - print '\tloadGeometry WARNING: could not read dy from the file name, assuming 10.0!' - dy = 10.0 - - ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) - # init geometry and mag. field - ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) - # -----Create geometry---------------------------------------------- - import shipDet_conf - - tgeom = ROOT.TGeoManager("Geometry", "Geane geometry") - geofile = inputFile.replace('ship.','geofile_full.').replace('_rec.','.') - print geofile - gMan = tgeom.Import(geofile) - fGeo = ROOT.gGeoManager - return {'fGeo':fGeo,'gMan':gMan, 'ShipGeo':ShipGeo} - -def getNode(nodeName,fGeo=None): - if fGeo is None: - r = loadGeometry(inputFile) - fGeo = r['fGeo'] - tv = fGeo.GetTopVolume() - nodes = tv.GetNodes() - dicts = prepareNodeDictionaries(nodes) - return nodes[dicts['nodeDict_name'][nodeName]] - -def findDimentionBoxElement(node): - ## GetDZ() etc gives you half of the dimention - sh = node.GetVolume().GetShape() - return {'x':sh.GetDX(), - 'y':sh.GetDY(), - 'z':sh.GetDZ()} - -def findPositionElement(node): - pos = node.GetMatrix().GetTranslation() - - return {'x':pos[0], - 'y':pos[1], - 'z':pos[2]} - -## inputFile: file used from which I would like to retrive the geometry used -## myNodes_name: geometrical elements I would like to analyse (find the position) -def findPositionGeoElement(inputFile, myNodes_name): - - r = loadGeometry(inputFile) - fGeo = r['fGeo'] - #volumes = gMan.GetListOfVolumes() - #for v in volumes: - # print v.GetName(), v.GetNumber() - - ## Get the top volume - - tv = fGeo.GetTopVolume() - nodes = tv.GetNodes() - - tmp = prepareNodeDictionaries(nodes) - nodeDict_name = tmp['nodeDict_name'] - res = {} - for nd_name in myNodes_name: - nd_index = nodeDict_name[nd_name] - nd = nodes[nd_index] - pos = findPositionElement(nd) - dims = findDimentionBoxElement(nd) - res[nd_name]={'z':pos['z'], - 'dimZ':dims['z'], - 'node':nd} - return res - -#inputFile = "../data/neutrino661/ship.10.0.Genie-TGeant4_D.root" -##searchForNodes(inputFile) - -### nodes name to be looked at -#myNodes_name = ["volLayer2_%s"%i for i in xrange(0,12)] -#myGeoEl = findPositionGeoElement(inputFile, myNodes_name) -#print myGeoEl - -##print myGeoEl["volLayer2_0"]['z'],myGeoEl["volLayer2_11"]['z'] -##print myGeoEl["volLayer2_11"]['z']-myGeoEl["volLayer2_0"]['z'] - +import ROOT,os,sys,getopt,time +import shipunit as u +import shipRoot_conf +from ShipGeoConfig import ConfigRegistry +from ROOT import * + +# function to prepare two dictionaries: one that has the node names as key and the related intex +# in the nodes array as value and the other one with key and value swapped. +def prepareNodeDictionaries(nodes): + # key: index, value: name of the node + nodeDict_index = {} + # key: name of the node, value: index + nodeDict_name = {} + for (i,n) in enumerate(nodes): + #print i,n.GetName() + nodeDict_index[i]=n.GetName() + nodeDict_name[n.GetName()]=i + return {'nodeDict_index':nodeDict_index, 'nodeDict_name':nodeDict_name} + +# function to print the node names and their index. Usefull to search the name of the element you would like +# to study if you don't know it by heart. Then the name can be used with the other functions for the specific studies. +def searchForNodes(inputFile, volName = None): + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + ## Get the top volume + fGeo = ROOT.gGeoManager + + if volName is None: + tv = fGeo.GetTopVolume() + else: + tv = fGeo.GetVolume(volName) + + nodes = tv.GetNodes() + for (i,n) in enumerate(nodes): + print i,n.GetName() + +def searchForNodes2(inputFile,volName=None): + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + ## Get the top volume + #fGeo = ROOT.gGeoManager + if volName is None: + tv = fGeo.GetTopVolume() + else: + tv = fGeo.GetVolume(volName) + + nodes = tv.GetNodes() + for (i,n) in enumerate(nodes): + print i,n.GetName(),findPositionElement(n)['z'],findDimentionBoxElement(n)['z'] + +def searchForNodes2_xyz(inputFile, volName=None): + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + ## Get the top volume + #fGeo = ROOT.gGeoManager + if volName is None: + tv = fGeo.GetTopVolume() + else: + tv = fGeo.GetVolume(volName) + nodes = tv.GetNodes() + for (i,n) in enumerate(nodes): + print i,n.GetName() + print " x: ", findPositionElement(n)['x'],findDimentionBoxElement(n)['x'] + print " y: ", findPositionElement(n)['y'],findDimentionBoxElement(n)['y'] + print " z: ", findPositionElement(n)['z'],findDimentionBoxElement(n)['z'] + +# basic function to be called to load the geometry from a file +def loadGeometry(inputFile): + dy = float( inputFile.split("/")[-1].split(".")[1]) + + ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) + # init geometry and mag. field + ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) + # -----Create geometry---------------------------------------------- + import shipDet_conf + + tgeom = ROOT.TGeoManager("Geometry", "Geane geometry") + geofile = inputFile.replace('ship.','geofile_full.').replace('_rec.','.') + print geofile + gMan = tgeom.Import(geofile) + fGeo = ROOT.gGeoManager + return {'fGeo':fGeo,'gMan':gMan, 'ShipGeo':ShipGeo} + +def getNode(nodeName,fGeo=None, volName=None): + if fGeo is None: + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + if volName is None: + tv = fGeo.GetTopVolume() + else: + tv = fGeo.GetVolume(volName) + + nodes = tv.GetNodes() + dicts = prepareNodeDictionaries(nodes) + return nodes[dicts['nodeDict_name'][nodeName]] + +def findDimentionBoxElement(node): + ## GetDZ() etc gives you half of the dimention + sh = node.GetVolume().GetShape() + return {'x':sh.GetDX(), + 'y':sh.GetDY(), + 'z':sh.GetDZ()} + +def findPositionElement(node): + pos = node.GetMatrix().GetTranslation() + + return {'x':pos[0], + 'y':pos[1], + 'z':pos[2]} + +## inputFile: file used from which I would like to retrive the geometry used +## myNodes_name: geometrical elements I would like to analyse (find the position) +def findPositionGeoElement(inputFile, myNodes_name, volName=None): + + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + #volumes = gMan.GetListOfVolumes() + #for v in volumes: + # print v.GetName(), v.GetNumber() + + ## Get the top volume + if volName is None: + tv = fGeo.GetTopVolume() + else: + tv = fGeo.GetVolume(volName) + + nodes = tv.GetNodes() + + tmp = prepareNodeDictionaries(nodes) + nodeDict_name = tmp['nodeDict_name'] + res = {} + for nd_name in myNodes_name: + nd_index = nodeDict_name[nd_name] + nd = nodes[nd_index] + pos = findPositionElement(nd) + dims = findDimentionBoxElement(nd) + res[nd_name]={'z':pos['z'], + 'dimZ':dims['z'], + 'node':nd, + 'dimX':dims['x'], + 'dimY':dims['y']} + return res + +#inputFile = "../data/neutrino661/ship.10.0.Genie-TGeant4_D.root" +##searchForNodes(inputFile) + +### nodes name to be looked at +#myNodes_name = ["volLayer2_%s"%i for i in xrange(0,12)] +#myGeoEl = findPositionGeoElement(inputFile, myNodes_name) +#print myGeoEl + +##print myGeoEl["volLayer2_0"]['z'],myGeoEl["volLayer2_11"]['z'] +##print myGeoEl["volLayer2_11"]['z']-myGeoEl["volLayer2_0"]['z'] + diff --git a/lookAtGeo_OLD.py b/lookAtGeo_OLD.py new file mode 100755 index 0000000..6f51cc6 --- /dev/null +++ b/lookAtGeo_OLD.py @@ -0,0 +1,137 @@ +import ROOT,os,sys,getopt,time +import shipunit as u +import shipRoot_conf +from ShipGeoConfig import ConfigRegistry +from ROOT import * + +# function to prepare two dictionaries: one that has the node names as key and the related intex +# in the nodes array as value and the other one with key and value swapped. +def prepareNodeDictionaries(nodes): + # key: index, value: name of the node + nodeDict_index = {} + # key: name of the node, value: index + nodeDict_name = {} + for (i,n) in enumerate(nodes): + #print i,n.GetName() + nodeDict_index[i]=n.GetName() + nodeDict_name[n.GetName()]=i + return {'nodeDict_index':nodeDict_index, 'nodeDict_name':nodeDict_name} + +# function to print the node names and their index. Usefull to search the name of the element you would like +# to study if you don't know it by heart. Then the name can be used with the other functions for the specific studies. +def searchForNodes(inputFile): + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + ## Get the top volume + fGeo = ROOT.gGeoManager + tv = fGeo.GetTopVolume() + nodes = tv.GetNodes() + for (i,n) in enumerate(nodes): + print i,n.GetName() + +def searchForNodes2(inputFile): + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + ## Get the top volume + #fGeo = ROOT.gGeoManager + tv = fGeo.GetTopVolume() + nodes = tv.GetNodes() + for (i,n) in enumerate(nodes): + print i,n.GetName(),findPositionElement(n)['z'],findDimentionBoxElement(n)['z'] + +def searchForNodes2_xyz(inputFile): + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + ## Get the top volume + #fGeo = ROOT.gGeoManager + tv = fGeo.GetTopVolume() + nodes = tv.GetNodes() + for (i,n) in enumerate(nodes): + print i,n.GetName() + print " x: ", findPositionElement(n)['x'],findDimentionBoxElement(n)['x'] + print " y: ", findPositionElement(n)['y'],findDimentionBoxElement(n)['y'] + print " z: ", findPositionElement(n)['z'],findDimentionBoxElement(n)['z'] + +# basic function to be called to load the geometry from a file +def loadGeometry(inputFile): + try: + dy = float( inputFile.split("/")[-1].split(".")[1]) + except: + print '\tloadGeometry WARNING: could not read dy from the file name, assuming 10.0!' + dy = 10.0 + + ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) + # init geometry and mag. field + ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) + # -----Create geometry---------------------------------------------- + import shipDet_conf + + tgeom = ROOT.TGeoManager("Geometry", "Geane geometry") + geofile = inputFile.replace('ship.','geofile_full.').replace('_rec.','.') + print geofile + gMan = tgeom.Import(geofile) + fGeo = ROOT.gGeoManager + return {'fGeo':fGeo,'gMan':gMan, 'ShipGeo':ShipGeo} + +def getNode(nodeName,fGeo=None): + if fGeo is None: + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + tv = fGeo.GetTopVolume() + nodes = tv.GetNodes() + dicts = prepareNodeDictionaries(nodes) + return nodes[dicts['nodeDict_name'][nodeName]] + +def findDimentionBoxElement(node): + ## GetDZ() etc gives you half of the dimention + sh = node.GetVolume().GetShape() + return {'x':sh.GetDX(), + 'y':sh.GetDY(), + 'z':sh.GetDZ()} + +def findPositionElement(node): + pos = node.GetMatrix().GetTranslation() + + return {'x':pos[0], + 'y':pos[1], + 'z':pos[2]} + +## inputFile: file used from which I would like to retrive the geometry used +## myNodes_name: geometrical elements I would like to analyse (find the position) +def findPositionGeoElement(inputFile, myNodes_name): + + r = loadGeometry(inputFile) + fGeo = r['fGeo'] + #volumes = gMan.GetListOfVolumes() + #for v in volumes: + # print v.GetName(), v.GetNumber() + + ## Get the top volume + + tv = fGeo.GetTopVolume() + nodes = tv.GetNodes() + + tmp = prepareNodeDictionaries(nodes) + nodeDict_name = tmp['nodeDict_name'] + res = {} + for nd_name in myNodes_name: + nd_index = nodeDict_name[nd_name] + nd = nodes[nd_index] + pos = findPositionElement(nd) + dims = findDimentionBoxElement(nd) + res[nd_name]={'z':pos['z'], + 'dimZ':dims['z'], + 'node':nd} + return res + +#inputFile = "../data/neutrino661/ship.10.0.Genie-TGeant4_D.root" +##searchForNodes(inputFile) + +### nodes name to be looked at +#myNodes_name = ["volLayer2_%s"%i for i in xrange(0,12)] +#myGeoEl = findPositionGeoElement(inputFile, myNodes_name) +#print myGeoEl + +##print myGeoEl["volLayer2_0"]['z'],myGeoEl["volLayer2_11"]['z'] +##print myGeoEl["volLayer2_11"]['z']-myGeoEl["volLayer2_0"]['z'] + diff --git a/ntuple_nuVeto_small.py b/ntuple_nuVeto_small.py new file mode 100644 index 0000000..2035142 --- /dev/null +++ b/ntuple_nuVeto_small.py @@ -0,0 +1,613 @@ +import sys, re, getopt +from array import array +from ROOT import * + +fileName, fileNameGeo, outputFileName, sampleType = None, None, None, None +maxevents = 0 +verbose = True + +# Parse commandline inputs +try: + opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:v:", []) +except getopt.GetoptError: + # print help information and exit: + print '\tUSAGE: -f inputfile.root -o outputfile.root -t sampleType (sig, nuBg or cosmics) -n maxevents -v verbose (0 or 1)' + sys.exit() +for o, a in opts: + if o in ("-f"): + fileName = a + if o in ("-o"): + outputFileName = a + if o in ("-n"): + maxevents = int(a) + if o in ("-t"): + sampleType = str(a) + if o in ("-v"): + verbose = bool(int(a)) + +fileNameGeo = fileName.replace('ship', 'geofile_full') + +def namestr(obj, namespace=globals()): + return [name for name in namespace if namespace[name] is obj] + +for item in [fileName, fileNameGeo, outputFileName, sampleType]: + if not item: + print '\tFATAL! %s not defined!'%namestr(item)[0] + sys.exit() + +for item in [fileName, fileNameGeo, outputFileName, sampleType, maxevents, verbose]: + print '\tINFO: %s set to %s'%(namestr(item)[0], item) + +# Load SHiP (LHCb) style +gROOT.ProcessLine(".x mystyle.C") +# Obsolete debug flags +oldGeo = False +ON = True +# Load PDG database +pdg = TDatabasePDG.Instance() +# Add elements to PDG database +import pythia8_conf +pythia8_conf.addHNLtoROOT() +# TODO: MAKE THIS A DICTIONARY AND DE-VERBOSIZE IT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 +if not(pdg.GetParticle(1000010020)): + pdg.AddParticle("Deuteron","Deuteron", 1.875613e+00, kTRUE, 0,3,"Ion",1000010020) + pdg.AddAntiParticle("AntiDeuteron", - 1000010020) +if not(pdg.GetParticle(1000010030)): + pdg.AddParticle("Triton","Triton", 2.808921e+00, kFALSE, 3.885235e+17,3,"Ion",1000010030); + pdg.AddAntiParticle("AntiTriton", - 1000010030); +if not(pdg.GetParticle(1000020040) ): + pdg.AddParticle("Alpha","Alpha",3.727379e+00,kFALSE,0,6,"Ion",1000020040); + pdg.AddAntiParticle("AntiAlpha", - 1000020040); +if not(pdg.GetParticle(1000020030) ): + pdg.AddParticle("HE3","HE3",2.808391e+00,kFALSE,0,6,"Ion",1000020030); + pdg.AddAntiParticle("AntiHE3", -1000020030); +if not(pdg.GetParticle(1000030070) ): + print "TO BE CHECKED the data insert for Li7Nucleus" + pdg.AddParticle("Li7Nucleus","Li7Nucleus",3.727379e+00/4.*7.,kFALSE,0,0,"Boson",1000030070); +if not(pdg.GetParticle(1000060120) ): + print "ERROR: random values insert for C12Nucleus" + pdg.AddParticle("C12Nucleus","C12Nucleus",0.1,kFALSE,0,0,"Isotope",1000060120); +if not(pdg.GetParticle(1000030060) ): + print "ERROR: random values insert for Li6Nucleus" + pdg.AddParticle("Li6Nucleus","Li6Nucleus",6.015121,kFALSE,0,0,"Isotope",1000030060); +if not(pdg.GetParticle(1000070140) ): + print "ERROR: random values insert for for N14" + pdg.AddParticle("N14","N14",0.1,kTRUE,0,0,"Isotope",1000070140); +if not(pdg.GetParticle(1000050100) ): + print "TO BE CHECKED the data insert for B10" + pdg.AddParticle("B10","B10",10.0129370,kTRUE,0,0,"Isotope",1000050100); +if not(pdg.GetParticle(1000020060) ): + print "TO BE CHECKED the data insert for He6" + pdg.AddParticle("He6","He6",6.0188891,kFALSE,806.7e-3,0,"Isotope",1000020060); +if not(pdg.GetParticle(1000040080) ): + print "TO BE CHECKED the data insert for Be8" + pdg.AddParticle("Be8","Be8",8.00530510,kFALSE,6.7e-17,0,"Isotope",1000040080); +if not(pdg.GetParticle(1000030080) ): + print "TO BE CHECKED the data insert for Li8" + pdg.AddParticle("Li8","Li8",8.002248736,kTRUE,178.3e-3,0,"Isotope",1000030080); +if not(pdg.GetParticle(1000040170) ): + print "ERROR: didn't find what is it particle with code 1000040170, random number inserted!" + pdg.AddParticle("None","None",0.1,kFALSE,0,0,"Isotope",1000040170); +if not(pdg.GetParticle(1000040100) ): + print "TO BE CHECKED the data insert for Be10" + pdg.AddParticle("Be10","Be10",10.0135338,kTRUE,5.004e+9,0,"Isotope",1000040100); +if not(pdg.GetParticle(1000040070) ): + print "TO BE CHECKED the data insert for Be7" + pdg.AddParticle("Be7","Be7",11.021658,kTRUE,13.81,0,"Isotope",1000040070); +if not(pdg.GetParticle(1000230470) ): + print "ERROR: didn't find what is it particle with code 1000230470, random number inserted!" + pdg.AddParticle("None2","None2",0.1,kFALSE,0,0,"Isotope",1000230470); +if not(pdg.GetParticle(1000080170) ): + print "ERROR: didn't find what is it particle with code 1000080170, random number inserted!" + pdg.AddParticle("None3","None3",0.1,kFALSE,0,0,"Isotope",1000080170); +if not(pdg.GetParticle(1000240500) ): + print "ERROR: didn't find what is it particle with code 1000240500, random number inserted!" + pdg.AddParticle("None4","None4",0.1,kFALSE,0,0,"Isotope",1000240500); +if not(pdg.GetParticle(1000210450) ): + print "ERROR: didn't find what is it particle with code 1000210450, random number inserted (Sc45Nucleus)!" + pdg.AddParticle("Sc45Nucleus","Sc45Nucleus",0.1,kFALSE,0,0,"Isotope",1000210450); +if not(pdg.GetParticle(1000040090) ): + print "TO BE CHECKED the data insert for Be9" + pdg.AddParticle("Be9","Be9",9.0121822,kFALSE,0,0,"Isotope",1000040090); +if not(pdg.GetParticle(1000080160) ): + print "TO BE CHECKED the data insert for O16" + pdg.AddParticle("O16","O16",15.99491461956,kFALSE,0,0,"Isotope",1000080160); +if not(pdg.GetParticle(1000220460) ): + print "TO BE CHECKED the data insert for Ar40" + pdg.AddParticle("Ar40","Ar40",39.9623831225,kFALSE,0,0,"Isotope",1000220460); +if not(pdg.GetParticle(1000050110) ): + print "TO BE CHECKED the data insert for B11" + pdg.AddParticle("B11","B11",11.0093054,kFALSE,0,0,"Isotope",1000050110); + +def PrintEventPart(particles,pdg): + print "Particles in the event:" + for (pid,part) in enumerate(particles): + print pid, pdg.GetParticle(part.GetPdgCode()).GetName(), part.GetMotherId() + print + +def PrintRPChits(rpc,pdg): + print "Hits in:" + for (ix,RPCpt) in enumerate(rpc): + tmpName = pdg.GetParticle(RPCpt.PdgCode()) + print ix,RPCpt.GetZ(), tmpName, RPCpt.GetTrackID() , fGeo.FindNode(RPCpt.GetX(),RPCpt.GetY(),RPCpt.GetZ()).GetVolume().GetName() + print + +# weight computation moved to offlineForBarbara.py + +def addVect(t,name,vectType): + vect = vector(vectType)() + t.Branch( name, vect ) + return t, vect + +def addVar(t,name,varType): + var = array(varType,[-999]) + t.Branch(name,var,name+"/"+varType.upper()) + return t, var + +def putToZero(var): + var[0] = 0 + +def getPartName(partId): + return pdg.GetParticle(partId).GetName() + +def lookingForDecay(particles): + decayMotherList = [] + for p in particles: + mID = p.GetMotherId() + if mID>=0 and not (mID in decayMotherList): + decayMotherList.append(mID) + decayPartIndex.push_back(mID) + decayPartID.push_back(particles[mID].GetPdgCode()) + decayPartStrID.push_back(pdg.GetParticle(particles[mID].GetPdgCode()).GetName()) + decayPos_x.push_back(p.GetStartX()) + decayPos_y.push_back(p.GetStartY()) + decayPos_z.push_back(p.GetStartZ()) + decayMother.push_back(particles[mID].GetMotherId()) + decayP.push_back(p.GetP()) + +def wasFired(indexKids, detPoints, detPos, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0): + # + def lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId): + charges = [] + trackids = [] + if hitCharges is None: + assert(hitTrackId is None) + for pos in detPos: + foundStat = False + foundStat_eff = False + for hit in detPoints: + if (indexKids is None) or (hit.GetTrackID() in indexKids): + # check if it is in one of the considered active layer + if pos[0]<=hit.GetZ()<=pos[1]: + Eloss += hit.GetEnergyLoss() + hitID = hit.GetTrackID() + if not hitID in trackids and not hitCharges is None: + if hit.PdgCode()>100000: + charges.append(9) + else: + charges.append(int(pdg.GetParticle(hit.PdgCode()).Charge())) + trackids.append(hitID) + hitCharges.push_back(charges[-1]) + hitTrackId.push_back(trackids[-1]) + if pointVects is not None: + pointVects[0].push_back(hit.GetX()) + pointVects[1].push_back(hit.GetY()) + pointVects[2].push_back(hit.GetZ()) + flag = True + foundStat = True + eff_val = gRandom.Uniform(0.,1.) + if eff_val<0.9: + flag_eff = flag_eff or True + foundStat_eff = True + else: + flag_eff = flag_eff or False + if foundStat: + nstat+=1 + if foundStat_eff: + nstat_eff+=1 + particles = 0 + flag_Ethr = Eloss>=Ethr + return flag, flag_eff, nstat, nstat_eff, Eloss, flag_Ethr,particles + # + # Now in partKidTrackID I should have the trackID connected to my charged particles + flag_eff = False + flag = False + nstat = 0 + nstat_eff = 0 + Eloss = 0 + flag,flag_eff,nstat,nstat_eff,Eloss,flag_Ethr,particles = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId) + if flag==False and checkOn: + print "To be study event %s"%entry + return flag, flag_eff, nstat, nstat_eff, flag_Ethr, particles + +def convertion(tmpName): + if "Rib" in tmpName: tmpName = "Rib" + elif "LiSc" in tmpName: tmpName= "LiSc" + elif "Tr" in tmpName: tmpName = "Tr" + elif "Startplate" in tmpName: tmpName = "Startplate" + elif re.search("T\d+O", tmpName): tmpName = "TO" + elif re.search("T\d+I", tmpName): tmpName = "TI" + elif "Endplate" in tmpName: tmpName = "Endplate" + elif "volCoil" in tmpName: tmpName = "volCoil" + elif "volFeYoke" in tmpName: tmpName = "volFeYoke" + elif "gas" in tmpName: tmpName = "gas" + elif "wire" in tmpName: tmpName == "wire" + elif "VetoTimeDet" in tmpName: tmpName = "upstreamVeto" + elif "Veto" in tmpName: tmpName = "strawVeto" + return tmpName + +def wasFired_node(indexKids, detPoints, detVolNames, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0): + # + def lookingForHits(detPoints,flag,hitCharges, hitTrackId): + if hitCharges is None: + assert(hitTrackId is None) + for hit in detPoints: + if (indexKids is None) or (hit.GetTrackID() in indexKids): + tmpName = fGeo.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName() + tmpName = convertion(tmpName) + if hit.GetZ()>8000: + continue + if tmpName in detVolNames: + ## trick to remove the case of the tracking system the hits in the gas or straw from the straw-veto: + if 'trackingSystem' in detVolNames and hit.GetZ()<0: + continue + if "strawVetoSystem" in detVolNames and hit.GetZ()>0: + continue + if "upstreamVeto" in detVolNames and not(hit.GetZ()>= upstreamVetoPos[0][0] and hit.GetZ()<=upstreamVetoPos[0][1]): + continue + # check if it is in one of the considered active layer + if True: #pos[0]<=hit.GetZ()<=pos[1]: + flag = True + return flag + # + # Now in partKidTrackID I should have the trackID connected to my charged particles + flag = False + flag = lookingForHits(detPoints,flag,hitCharges, hitTrackId) + if flag==False and checkOn: + print "To be study event %s"%entry + return flag + + +""" Main program """ + +# Open file +f = TFile(fileName) +t = f.Get("cbmsim") +totentries = t.GetEntries() +if verbose: + print + print + print "Program started, reading %s from %s"%(t.GetName(), fileName) +if (maxevents>0) and (maxevents < totentries): + entries = maxevents +else: + entries = totentries +if verbose: + print "Requested to process %s events out of %s in tree"%(entries, totentries) + +# Load geometry +from lookAtGeo import * +r = loadGeometry(fileNameGeo) +fGeo = r['fGeo'] + +# Functions by Elena for offline selection +# (Also includes another dictionary of all geometry nodes) +import offlineForBarbara as offline +offline.loadNodes(fGeo) +offline.ShipGeo = r['ShipGeo'] + +# Handle geometry +# Names of useful geometry nodes +myNodes_name = ["VetoTimeDet_1"] +myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)] +myNodes_name += ["Veto_5"] +# Positions of selected nodes +myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name,None) +# Places where a neutrino can interact +lastPassive_nodeName = ["volIron_23"] +OPERA_nodeName = ["volIron", "volRPC", "volHPT","volArm2MS","volFeYokes","volCoil","volFeYoke"]#["volIron_%s"%i for i in xrange(12,24)]+["volRpc_%s"%i for i in xrange(11,22)]+["volArm2MS_1"] +entrance_nodeName = ["T1Lid"]#['T1lid_1'] +volumeIn_nodeName = ["TI"]#['T%sI'%i for i in [1,2,3,5]] +volumeOut_nodeName = ["TO"]#['T%sO'%i for i in [1,2,3,5]] +OPERA_volNames = ["volIron","volRpc","volHPT"] +strawVeto_volNames = ["Veto"] +# Z positions of the VETO and tracking systems +Tracking = [myGeoEl["Tr1_1"]['z']-myGeoEl["Tr1_1"]['dimZ'],myGeoEl["Tr4_4"]['z']+myGeoEl["Tr4_4"]['dimZ']] +upstreamVeto = [myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']] +trackStationsPos = [[myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ'],myGeoEl["Tr%s_%s"%(i,i)]['z']+myGeoEl["Tr%s_%s"%(i,i)]['dimZ']] for i in xrange(1,5)] +strawVetoPos = [[myGeoEl["Veto_5"]['z']-myGeoEl["Veto_5"]['dimZ'],myGeoEl["Veto_5"]['z']+myGeoEl["Veto_5"]['dimZ']]] +vetoWall = [[myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+0.001,myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ']]]#myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+6000.]] +upstreamVetoPos = [[myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]] +RPCstationsPos = [[-3500,myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ']-0.5]] +# Places where a neutrino can interact (some duplicates?) +trackStationsPos_node = ["Tr", "gas", "straw", 'trackingSystem'] +## should be that this does not work since it could be that hits are also assigned to gas or straw +strawVetoPos_node = ["strawVeto", "gas", "straw", "strawVetoSystem"] +vetoWall_node = ["LiSc","cave","Rib","TI","TO","Tr","strawVeto"] +upstreamVetoPos_node = ["upstreamVeto","cave"] +RPCstationsPos_node = ["volRpc","volMagneticSpectrometer","volIron","volHPT"] +# Printout positions +print "OPERApos: ",RPCstationsPos +print "trackingPos: ",trackStationsPos +print "strawVetoPos: ",strawVetoPos +print "scintPos: ",vetoWall +print "upstreamVetoPos: ",upstreamVetoPos + +# Prepare output ntuple +nf = TFile(outputFileName,"RECREATE") +nt = TTree("t","t") +# Event number -------------------------------------------- +nt, event = addVar(nt, 'event','i') +# Neutrino information ------------------------------------ +nt, nNu = addVar(nt, 'nNu', 'i') +nt, isPrimary = addVar(nt, 'isPrimary','i') +nt, startZ_nu = addVar(nt, 'startZ_nu', 'f') +nt, startY_nu = addVar(nt, 'startY_nu', 'f') +nt, startX_nu = addVar(nt, 'startX_nu', 'f') +nt, nuE = addVar(nt, 'nuE', 'f') +nt, weight = addVar(nt, 'weight', 'f') +# Interaction element code -------------------------------- +## 0: last passive material opera-mu-system +## 1: two windows around liquid scintillator +## 2: along vaccum tank +## 3: along all OPERA +## 4: between the two entrance windows +## 5: along vacuum tank outer window +## 6: along vacuum tank inner window +## 999: wrong OPERA place ## needed until we have the new production +## -1: anything else, but what???? #it should never been present +nt, interactionElement = addVar(nt, 'interactionElement','i') +# Neutrino daughters -------------------------------------- +# Do the particles come from a neutrino? +nt, nPart_fromNu = addVar(nt, 'nPart_fromNu', 'i') +nt, idPart_fromNu = addVect(nt, 'idPart_fromNu', 'int') +#nt, idStrPart_fromNu = addVect(nt, 'idStrPart_fromNu', 'string') +# Do charged particles come from a neutrino? +nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i') +nt, idChargedPart_fromNu = addVect(nt, 'idChargedPart_fromNu', 'int') +#nt, idStrChargedPart_fromNu = addVect(nt, 'idStrChargedPart_fromNu', 'string') +# Do neutral particles come from a neutrino? +nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i') +nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int') +#nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string') +# RPC ----------------------------------------------------- +# In case something (no matter if it comes from the nu-interaction kids) fired the RPC +nt, RPCany = addVar(nt,'RPCany', 'i' ) +# accounting for an eff. of each station of 90% +nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' ) +# upstreamVeto -------------------------------------------- +nt, upstreamVetoany = addVar(nt,'upstreamVetoAny', 'i' ) +# accounting for an eff. of each station of 90% +nt, upstreamVetoany_eff = addVar(nt,'upstreamVetoAny_eff', 'i' ) +# scintVeto ----------------------------------------------- +## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls +nt, scintVetoAny = addVar(nt,'scintVetoAny', 'i' ) +nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' ) +# strawVeto ----------------------------------------------- +# In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto +nt, strawVetoAny = addVar(nt,'strawVetoAny', 'i' ) +nt, strawVetoAny_eff = addVar(nt,'strawVetoAny_eff', 'i' ) +# TrackSyst ----------------------------------------------- +# only the case of anything fired it is considered, obviously +nt, TrackSyst = addVar(nt, "TrackSyst", 'i') +nt, TrackSyst_eff = addVar(nt, "TrackSyst_eff", 'i') +# VETO and tracking systems with thresholds --------------- +nt, strawVetoAny_Ethr = addVar(nt,"strawVetoAny_Ethr","i") +nt, scintVetoAny_Ethr = addVar(nt,"scintVetoAny_Ethr","i") +nt, upstreamVetoAny_Ethr = addVar(nt,"upstreamVetoAny_Ethr","i") +nt, RPCany_Ethr = addVar(nt,"RPCany_Ethr","i") +nt, TrackSyst_Ethr = addVar(nt, "TrackSyst_Ethr", "i") +# Decay information --------------------------------------- +nt, decayPartIndex = addVect(nt, "decayPartIndex", "float") +nt, decayPartID = addVect(nt, "decayPartID", "float") +nt, decayPartStrID = addVect(nt, "decayPartStrID", "string") +nt, decayPos_x = addVect(nt, "decayPos_x", "float") +nt, decayPos_y = addVect(nt, "decayPos_y", "float") +nt, decayPos_z = addVect(nt, "decayPos_z", "float") +nt, decayMother = addVect(nt, "decayMother", "float") +nt, decayP = addVect(nt,"decayP","float") +# Is this a NC interaction? ------------------------------- +nt, NC = addVar(nt, "NC", "i") +# Store where the neutrino interacted --------------------- +#nt, nuInteractionNode = addVect(nt, "nuInteractionNode", "string") +nt, nuIntNumSimpl = addVar(nt,"nuIntNumSimpl","i") +dictNodeNames = {'volIron':0, 'cave':1, 'LiSc':2, 'Startplate':3, 'TI':4, 'rockD':5, 'Endplate':6, 'Rib':7, 'volFeYoke':8, 'volHPT':9, 'TO':10, 'volRpc':11, 'volCoil':12, 'T1Lid':13, 'straw':14, 'strawVeto':15, 'volBase':16, 'Tr':17, 'gas':18, 'wire':19, 'others':20} +# Was the event reconstructed? ---------------------------- +nt, recoed = addVar(nt, "recoed","i") +# How many candidates in the reco event? ------------------ +nt, nRecoed = addVar(nt, "nRecoed","i") +# Add stuff for online selection -------------------------- +offline.addOfflineToTree(nt) + +# Now loop on the tree +# t = original tree +# nt = new tree +for entry in xrange(entries): + if not (entry%1000): + print "\tProcessing entry %s of %s..."%(entry,entries) + t.GetEntry(entry) + event[0] = entry + particles = t.MCTrack + rpc = t.ShipRpcPoint + scint = t.vetoPoint + strawPoints = t.strawtubesPoint + vetoScint = t.vetoPoint + recoParts = t.Particles + # initialize containers + putToZero(nNu) + putToZero(nPart_fromNu) + putToZero(nChargedPart_fromNu) + putToZero(nNeutrPart_fromNu) + idPart_fromNu.clear() + idChargedPart_fromNu.clear() + idNeutrPart_fromNu.clear() + #idStrPart_fromNu.clear() + #idStrChargedPart_fromNu.clear() + #idStrNeutrPart_fromNu.clear() + decayPartIndex.clear() + decayPartID.clear() + decayPartStrID.clear() + decayPos_x.clear() + decayPos_y.clear() + decayPos_z.clear() + decayMother.clear() + decayP.clear() + #nuInteractionNode.clear() + primaryDone=False + isPrimary[0] = int(False) + lookingForDecay(particles) + nRecoed[0] = 0 + recoed[0] = 0 + # Which one is the "primary" particle? + if sampleType=="nuBg" or sampleType == "cosmics": + startPartRange = [0] + primaryMum = -1 + elif sampleType=="sig": + startPartRange = [1] + primaryMum = 0 + # Look at the primary particle + for ip in startPartRange: + # Were there any hit in the tracking stations? + TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, None, None, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] ) + # If no, skip this event + if TrackSyst[0]==0: + skipEvent=True + continue + skipEvent=False + # Select the primary MC particle + part = particles[ip] + pdgPart = pdg.GetParticle(part.GetPdgCode()) + if sampleType=="nuBg": + assert("nu" in pdgPart.GetName()) + # Add information for offline selection + # (mainly signal normalisation and zeroing of particle-wise arrays) + offline.pushOfflineByEvent(t, sampleType, verbose) + ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother + #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena + if True: + ## Starting the counter of how many particles were produced by the interaction of this specific nu + for recoP in recoParts: + nRecoed[0] += 1 + offline.pushOfflineByParticle(t, recoP) + if nRecoed[0]>0: + recoed[0]=1 + #else: + #skipEvent=True + #continue + nNu[0]+=1 + if part.GetMotherId()==primaryMum: + assert(primaryDone==False) + primaryDone=True + isPrimary[0] = int(True) + assert(len(idPart_fromNu)==0) + assert(len(idNeutrPart_fromNu)==0) + else: + continue + # Where did this guy interact? + tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName() + #nuInteractionNode.push_back(tmpName) + tmpName = convertion(tmpName) + if not tmpName in dictNodeNames: + tmpName = "others" + nuIntNumSimpl[0] = dictNodeNames[tmpName] + # Store interaction point coordinates + startZ_nu[0] = part.GetStartZ() + startY_nu[0] = part.GetStartY() + startX_nu[0] = part.GetStartX() + # Primary particle information + nuE[0] = part.GetEnergy() + # Looking for particles produced by the neutrino interaction + interacted = False + partKidsId = [] + NC[0] = 0 + # Loop on MCTracks + for ip2 in xrange(0,len(particles)): + part2 = particles[ip2] + # exit if we have reached the empty part of the array + if not (type(part2)==type(ShipMCTrack())): + break + if part2.GetMotherId()==ip: + interacted = True + part2Id = part2.GetPdgCode() + partKidsId.append(part2Id) + # Was this a neutral current interaction? + if not (pdg.GetParticle(part2.GetPdgCode()) == None) and ("nu" in pdg.GetParticle(part2.GetPdgCode()).GetName()) and (part2.GetMotherId()==0): + NC[0] = 1 + nPart_fromNu[0]+=1 + idPart_fromNu.push_back(part2Id) + if pdg.GetParticle(part2.GetPdgCode()) == None: + part2Name = str(part2Id) + else: + part2Name = getPartName(part2Id) + #idStrPart_fromNu.push_back(part2Name) + # Counting how many particles produced by the nu-interaction are charged or neutral + if (pdg.GetParticle(part2.GetPdgCode())) == None or (int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0)): + nNeutrPart_fromNu[0]+=1 + idNeutrPart_fromNu.push_back(part2Id) + #idStrNeutrPart_fromNu.push_back(part2Name) + else: + nChargedPart_fromNu[0]+=1 + idChargedPart_fromNu.push_back(part2Id) + #idStrChargedPart_fromNu.push_back(part2Name) + # How many are produced by the interaction in the last passive element of the "opera-mu system"? + # Finding out the "interaction element code" + part2Z = part2.GetStartZ() + part2X = part2.GetStartX() + part2Y = part2.GetStartY() + somewhere = False + nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName() + if nodeName in lastPassive_nodeName: + somewhere = True + interactionElement[0] = 0 + intElName = 'OPERA' + # To know if it was in the full OPERA-system (excluded last passive) + if tmpName in OPERA_nodeName and somewhere==False: + somewhere = True + interactionElement[0] = 3 + intElName = 'OPERA' + # To know if it was between the two windows + if tmpName in entrance_nodeName: + assert(somewhere==False) + somewhere = True + interactionElement[0] = 1 + # Vacuum tank outer window + if tmpName in volumeIn_nodeName and somewhere==False: + interactionElement[0] = 5 + somewhere = True + # Vacuum tank inner window + elif nodeName in volumeOut_nodeName and somewhere==False: + interactionElement[0] = 6 + somewhere = True + # Vacuum tank ribs + if "Rib" in tmpName: + interactionElement[0] = 7 + # Somewhere else + if somewhere==False: + interactionElement[0] = -1 + # End loop on MCTracks + # Should be changed to avoid entering in the loop if there isn't anything in the tracking + if skipEvent: + continue + strawVetoAny[0],strawVetoAny_eff[0],dummy,dummy,strawVetoAny_Ethr[0],dummy = wasFired(None, strawPoints, strawVetoPos,None,None, checkOn=False, pointVects=None)#[strawVetoPoint_x, strawVetoPoint_y, strawVetoPoint_z]) + upstreamVetoany[0], upstreamVetoany_eff[0], dummy,dummy, upstreamVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, upstreamVetoPos, None, None, checkOn=False, pointVects=None)#[upstreamVetoPoint_x, upstreamVetoPoint_y, upstreamVetoPoint_z]) + RPCany[0], RPCany_eff[0],dummy,dummy, RPCany_Ethr[0],dummy = wasFired(None, rpc, RPCstationsPos, None, None, checkOn=False, pointVects=None)#[rpcPoint_x, rpcPoint_y, rpcPoint_z]) + scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=0.015) + #assert(len(idStrPart_fromNu)==len(idPart_fromNu)) + #assert(len(idStrChargedPart_fromNu)==len(idChargedPart_fromNu)) + #assert(len(idStrNeutrPart_fromNu)==len(idNeutrPart_fromNu)) + #assert(len(idStrPart_fromNu)==nChargedPart_fromNu[0]+nNeutrPart_fromNu[0]) + #assert(len(idStrChargedPart_fromNu)==nChargedPart_fromNu[0]) + #assert(len(idStrNeutrPart_fromNu)==nNeutrPart_fromNu[0]) + if not primaryDone: + print "It looks like there is no primary nu interacting" + PrintEventPart(particles,pdg) + sys.exit() + weight[0] = offline.findWeight(sampleType, NC[0], nuE[0], part, entries, pdgPart.GetName(), ON)#calcWeight(NC[0], nuE[0], part.GetWeight(), entries, pdgPart.GetName(), ON) + nt.Fill() + # End loop on tree +nt.Write() +nf.Save() +nf.Close() +f.Close() +print "Output wrote to %s"%outputFileName +print "Number of events with vertex outside Veto_5+500 - Tr4_4: %s"%offline.num_bad_z diff --git a/ntuple_nuVeto_small_BACKUP.py b/ntuple_nuVeto_small_BACKUP.py new file mode 100644 index 0000000..8aecd56 --- /dev/null +++ b/ntuple_nuVeto_small_BACKUP.py @@ -0,0 +1,911 @@ +from ROOT import * +from array import array +import re +import getopt + +try: + opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:", []) +except getopt.GetoptError: + # print help information and exit: + print '\tUSAGE: -f inputfile.root -o outputfile.root -t sampleType (sig, nuBg or cosmics) -n maxevents' + sys.exit() +for o, a in opts: + if o in ("-f"): + fileName = a + if o in ("-o"): + outputFileName = a + if o in ("-n"): + maxEvents = int(a) + if o in ("-t"): + sampleType = str(a) + +fileNameGeo = fileName.replace('ship', 'geofile_full') +gROOT.ProcessLine(".x mystyle.C") +pdg = TDatabasePDG.Instance() + +sampleType = 'nuBg' + +oldGeo = False +## are the weights to be switched on? +ON = True +nu = True + +#where = "../../data/74x_75x/" +#where = "/afs/cern.ch/work/b/bstora/data/" +# +#fileNameGeo = "%sgeofile_full.10.0.Genie-TGeant4.root"%where +# +#if nu: +# fileName = "%sship.10.0.Genie-TGeant4_rec_neutrino_newReco_nu_77_79.root"%where +# outputFileName = "%sfromAnalysis/ntupleStudy_highFlux_nu_77_77_higherThr.root"%where +#else: +# fileName = "%sship.10.0.Genie-TGeant4_rec_neutrino_newReco_antinu_76_78.root"%where +# outputFileName = "%sfromAnalysis/ntupleStudy_highFlux_antinu_76_78_higherThr.root"%where + + +if not(pdg.GetParticle(1000010020)): + pdg.AddParticle("Deuteron","Deuteron", 1.875613e+00, kTRUE, 0,3,"Ion",1000010020) + pdg.AddAntiParticle("AntiDeuteron", - 1000010020) + +if not(pdg.GetParticle(1000010030)): + pdg.AddParticle("Triton","Triton", 2.808921e+00, kFALSE, 3.885235e+17,3,"Ion",1000010030); + pdg.AddAntiParticle("AntiTriton", - 1000010030); + +if not(pdg.GetParticle(1000020040) ): + #pdg.AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,khShGev/(12.33*kYear2Sec),6,"Ion",ionCode); + pdg.AddParticle("Alpha","Alpha",3.727379e+00,kFALSE,0,6,"Ion",1000020040); + pdg.AddAntiParticle("AntiAlpha", - 1000020040); + +if not(pdg.GetParticle(1000020030) ): + pdg.AddParticle("HE3","HE3",2.808391e+00,kFALSE,0,6,"Ion",1000020030); + pdg.AddAntiParticle("AntiHE3", -1000020030); + +if not(pdg.GetParticle(1000030070) ): + print "TO BE CHECKED the data insert for Li7Nucleus" + pdg.AddParticle("Li7Nucleus","Li7Nucleus",3.727379e+00/4.*7.,kFALSE,0,0,"Boson",1000030070); + +if not(pdg.GetParticle(1000060120) ): + print "ERROR: random values insert for C12Nucleus" + pdg.AddParticle("C12Nucleus","C12Nucleus",0.1,kFALSE,0,0,"Isotope",1000060120); + +if not(pdg.GetParticle(1000030060) ): + print "ERROR: random values insert for Li6Nucleus" + pdg.AddParticle("Li6Nucleus","Li6Nucleus",6.015121,kFALSE,0,0,"Isotope",1000030060); + +if not(pdg.GetParticle(1000070140) ): + print "ERROR: random values insert for for N14" + pdg.AddParticle("N14","N14",0.1,kTRUE,0,0,"Isotope",1000070140); + +if not(pdg.GetParticle(1000050100) ): + print "TO BE CHECKED the data insert for B10" + pdg.AddParticle("B10","B10",10.0129370,kTRUE,0,0,"Isotope",1000050100); + +if not(pdg.GetParticle(1000020060) ): + print "TO BE CHECKED the data insert for He6" + pdg.AddParticle("He6","He6",6.0188891,kFALSE,806.7e-3,0,"Isotope",1000020060); + +if not(pdg.GetParticle(1000040080) ): + print "TO BE CHECKED the data insert for Be8" + pdg.AddParticle("Be8","Be8",8.00530510,kFALSE,6.7e-17,0,"Isotope",1000040080); + +if not(pdg.GetParticle(1000030080) ): + print "TO BE CHECKED the data insert for Li8" + pdg.AddParticle("Li8","Li8",8.002248736,kTRUE,178.3e-3,0,"Isotope",1000030080); + +if not(pdg.GetParticle(1000040170) ): + print "ERROR: didn't find what is it particle with code 1000040170, random number inserted!" + pdg.AddParticle("None","None",0.1,kFALSE,0,0,"Isotope",1000040170); + +if not(pdg.GetParticle(1000040100) ): + print "TO BE CHECKED the data insert for Be10" + pdg.AddParticle("Be10","Be10",10.0135338,kTRUE,5.004e+9,0,"Isotope",1000040100); + +if not(pdg.GetParticle(1000040070) ): + print "TO BE CHECKED the data insert for Be7" + pdg.AddParticle("Be7","Be7",11.021658,kTRUE,13.81,0,"Isotope",1000040070); + +if not(pdg.GetParticle(1000230470) ): + print "ERROR: didn't find what is it particle with code 1000230470, random number inserted!" + pdg.AddParticle("None2","None2",0.1,kFALSE,0,0,"Isotope",1000230470); + +if not(pdg.GetParticle(1000080170) ): + print "ERROR: didn't find what is it particle with code 1000080170, random number inserted!" + pdg.AddParticle("None3","None3",0.1,kFALSE,0,0,"Isotope",1000080170); + +if not(pdg.GetParticle(1000240500) ): + print "ERROR: didn't find what is it particle with code 1000240500, random number inserted!" + pdg.AddParticle("None4","None4",0.1,kFALSE,0,0,"Isotope",1000240500); + +if not(pdg.GetParticle(1000210450) ): + print "ERROR: didn't find what is it particle with code 1000210450, random number inserted (Sc45Nucleus)!" + pdg.AddParticle("Sc45Nucleus","Sc45Nucleus",0.1,kFALSE,0,0,"Isotope",1000210450); + +if not(pdg.GetParticle(1000040090) ): + print "TO BE CHECKED the data insert for Be9" + pdg.AddParticle("Be9","Be9",9.0121822,kFALSE,0,0,"Isotope",1000040090); + +if not(pdg.GetParticle(1000080160) ): + print "TO BE CHECKED the data insert for O16" + pdg.AddParticle("O16","O16",15.99491461956,kFALSE,0,0,"Isotope",1000080160); + +if not(pdg.GetParticle(1000220460) ): + print "TO BE CHECKED the data insert for Ar40" + pdg.AddParticle("Ar40","Ar40",39.9623831225,kFALSE,0,0,"Isotope",1000220460); + +if not(pdg.GetParticle(1000050110) ): + print "TO BE CHECKED the data insert for B11" + pdg.AddParticle("B11","B11",11.0093054,kFALSE,0,0,"Isotope",1000050110); + + + +def PrintEventPart(particles,pdg): + print "Particles in the event:" + for (pid,part) in enumerate(particles): + print pid, pdg.GetParticle(part.GetPdgCode()).GetName(), part.GetMotherId() + print + +def PrintRPChits(rpc,pdg): + print "Hits in:" + for (ix,RPCpt) in enumerate(rpc): + #if pdg.GetParticle(RPCpt.PdgCode())<10000: + # tmpName = pdg.GetParticle(RPCpt.PdgCode()).GetName() + #else: + tmpName = pdg.GetParticle(RPCpt.PdgCode()) + #if RPCpt.GetZ()>=-2484.0 and RPCpt.GetZ()<=-2482.0: + print ix,RPCpt.GetZ(), tmpName, RPCpt.GetTrackID() , fGeo.FindNode(RPCpt.GetX(),RPCpt.GetY(),RPCpt.GetZ()).GetVolume().GetName() + print + +KsID = 310 +KLID = 130 +NID = 2112 + +# weight computation moved to offlineForBarbara.py +""" +ff = TFile("histoForWeights.root") +h_GioHans = ff.Get("h_Gio") +# t.Draw("interactionElement", "(isPrimary==1)*weight") + +def calcWeight(NC, E, w, entries, nuName, ON=True): + if not ON: + return 1 + + if "bar" in nuName: + reduction = 0.5 + #flux = 1.88e+11 * 2.e+20 / 5.e+13 + flux = 6.98e+11 * 2.e+20 / 5.e+13 + else: + reduction = 1. + #flux = 2.88e+11 * 2.e+20/ 5.e+13 + flux = 1.09e+12 * 2.e+20/ 5.e+13 + + #reduction = 1. + #flux = 2.e+18 + #entries = 99999+99999. + crossSec = 6.7e-39*E * reduction + + + + #flux = 2.e+18 + NA = 6.022e+23 + + binN = h_GioHans.GetXaxis().FindBin(E) + return crossSec * flux * h_GioHans.GetBinContent(binN) * w * NA / entries +""" + +def addVect(t,name,vectType): + vect = vector(vectType)() + t.Branch( name, vect ) + return t, vect + +def addVar(t,name,varType): + var = array(varType,[-999]) + t.Branch(name,var,name+"/"+varType.upper()) + return t, var + +def putToZero(var): + var[0] = 0 + +def getPartName(partId): + return pdg.GetParticle(partId).GetName() + +def lookingForDecay(particles): + decayMotherList = [] + + for p in particles: + mID = p.GetMotherId() + if mID>=0 and not (mID in decayMotherList): + decayMotherList.append(mID) + decayPartIndex.push_back(mID) + decayPartID.push_back(particles[mID].GetPdgCode()) + decayPartStrID.push_back(pdg.GetParticle(particles[mID].GetPdgCode()).GetName()) + decayPos_x.push_back(p.GetStartX()) + decayPos_y.push_back(p.GetStartY()) + decayPos_z.push_back(p.GetStartZ()) + decayMother.push_back(particles[mID].GetMotherId()) + decayP.push_back(p.GetP()) + +def wasFired(indexKids, detPoints, detPos, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0): + + def lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId): + charges = [] + trackids = [] + if hitCharges is None: + assert(hitTrackId is None) + + for pos in detPos: + foundStat = False + foundStat_eff = False + + for hit in detPoints: + if (indexKids is None) or (hit.GetTrackID() in indexKids): + # check if it is in one of the considered active layer + if pos[0]<=hit.GetZ()<=pos[1]: + Eloss += hit.GetEnergyLoss() + hitID = hit.GetTrackID() + if not hitID in trackids and not hitCharges is None: + if hit.PdgCode()>100000: + charges.append(9) + else: + charges.append(int(pdg.GetParticle(hit.PdgCode()).Charge())) + trackids.append(hitID) + hitCharges.push_back(charges[-1]) + hitTrackId.push_back(trackids[-1]) + if pointVects is not None: + pointVects[0].push_back(hit.GetX()) + pointVects[1].push_back(hit.GetY()) + pointVects[2].push_back(hit.GetZ()) + + flag = True + #nstat+=1 + foundStat = True + eff_val = gRandom.Uniform(0.,1.) + if eff_val<0.9: + flag_eff = flag_eff or True + foundStat_eff = True + #nstat_eff+=1 + else: + flag_eff = flag_eff or False + if foundStat: + nstat+=1 + if foundStat_eff: + nstat_eff+=1 + + #charges.sort() + particles = 0 + #for i in xrange(int(len(charges)/2.)): + # if charges[0+i]*charges[-1-i]<0: + # particles+=1 + # else: + # break + flag_Ethr = Eloss>=Ethr + return flag, flag_eff, nstat, nstat_eff, Eloss, flag_Ethr,particles + + + # Now in partKidTrackID I should have the trackID connected to my charged particles + + flag_eff = False + flag = False + nstat = 0 + nstat_eff = 0 + Eloss = 0 + + flag,flag_eff,nstat,nstat_eff,Eloss,flag_Ethr,particles = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId) + + if flag==False and checkOn: + print "To be study event %s"%entry + + return flag, flag_eff, nstat, nstat_eff, flag_Ethr, particles + +def convertion(tmpName): + if "Rib" in tmpName: tmpName = "Rib" + elif "LiSc" in tmpName: tmpName= "LiSc" + elif "Tr" in tmpName: tmpName = "Tr" + elif "Startplate" in tmpName: tmpName = "Startplate" + elif re.search("T\d+O", tmpName): tmpName = "TO" + elif re.search("T\d+I", tmpName): tmpName = "TI" + elif "Endplate" in tmpName: tmpName = "Endplate" + #if "VetoTimeDet" in tmpName: tmpName = "upstreamVeto" + elif "volCoil" in tmpName: tmpName = "volCoil" + elif "volFeYoke" in tmpName: tmpName = "volFeYoke" + elif "gas" in tmpName: tmpName = "gas" + elif "wire" in tmpName: tmpName == "wire" + elif "VetoTimeDet" in tmpName: tmpName = "upstreamVeto" + elif "Veto" in tmpName: tmpName = "strawVeto" + return tmpName + +def wasFired_node(indexKids, detPoints, detVolNames, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0): + + def lookingForHits(detPoints,flag,hitCharges, hitTrackId): + if hitCharges is None: + assert(hitTrackId is None) + + for hit in detPoints: + if (indexKids is None) or (hit.GetTrackID() in indexKids): + #print hit.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName(), detVolNames + + tmpName = fGeo.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName() + tmpName = convertion(tmpName) + if hit.GetZ()>8000: + continue + #print tmpName, detVolNames + if tmpName in detVolNames: + ## trick to remove the case of the tracking system the hits in the gas or straw from the straw-veto: + if 'trackingSystem' in detVolNames and hit.GetZ()<0: + continue + if "strawVetoSystem" in detVolNames and hit.GetZ()>0: + continue + #if "upstreamVeto" in detVolNames and not(hit.GetZ()>=-2484.0 and hit.GetZ()<=-2482.0): + #upstreamVetoPos + if "upstreamVeto" in detVolNames and not(hit.GetZ()>= upstreamVetoPos[0][0] and hit.GetZ()<=upstreamVetoPos[0][1]): + continue + #print RPChit.GetTrackID(), t_ID + # check if it is in one of the considered active layer + if True: #pos[0]<=hit.GetZ()<=pos[1]: + flag = True + + return flag + + + # Now in partKidTrackID I should have the trackID connected to my charged particles + + flag = False + + flag = lookingForHits(detPoints,flag,hitCharges, hitTrackId) + + if flag==False and checkOn: + print "To be study event %s"%entry + + return flag + + + +f = TFile(fileName) +t = f.Get("cbmsim") + +entries = t.GetEntries() +from lookAtGeo import * + + +r = loadGeometry(fileNameGeo) +fGeo = r['fGeo'] + +# stuff by elena +from offlineForBarbara import * + +myNodes_name = ["VetoTimeDet_1"] +myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)] +myNodes_name += ["Veto_5"] + +myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name,None) + +lastPassive_nodeName = ["volIron_23"] +OPERA_nodeName = ["volIron", "volRPC", "volHPT","volArm2MS","volFeYokes","volCoil","volFeYoke"]#["volIron_%s"%i for i in xrange(12,24)]+["volRpc_%s"%i for i in xrange(11,22)]+["volArm2MS_1"] +entrance_nodeName = ["T1Lid"]#['T1lid_1'] +volumeIn_nodeName = ["TI"]#['T%sI'%i for i in [1,2,3,5]] +volumeOut_nodeName = ["TO"]#['T%sO'%i for i in [1,2,3,5]] + +OPERA_volNames = ["volIron","volRpc","volHPT"] +strawVeto_volNames = ["Veto"] + +top=fGeo.GetTopVolume() + +Tracking = [myGeoEl["Tr1_1"]['z']-myGeoEl["Tr1_1"]['dimZ'],myGeoEl["Tr4_4"]['z']+myGeoEl["Tr4_4"]['dimZ']] +upstreamVeto = [myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']] + +trackStationsPos = [[myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ'],myGeoEl["Tr%s_%s"%(i,i)]['z']+myGeoEl["Tr%s_%s"%(i,i)]['dimZ']] for i in xrange(1,5)] +strawVetoPos = [[myGeoEl["Veto_5"]['z']-myGeoEl["Veto_5"]['dimZ'],myGeoEl["Veto_5"]['z']+myGeoEl["Veto_5"]['dimZ']]] +vetoWall = [[myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+0.001,myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ']]]#myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+6000.]] +upstreamVetoPos = [[myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]] +RPCstationsPos = [[-3500,myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ']-0.5]] + + +trackStationsPos_node = ["Tr", "gas", "straw", 'trackingSystem'] +## should be that this does not work since it could be that hits are also assigned to gas or straw +strawVetoPos_node = ["strawVeto", "gas", "straw", "strawVetoSystem"] +vetoWall_node = ["LiSc","cave","Rib","TI","TO","Tr","strawVeto"] +upstreamVetoPos_node = ["upstreamVeto","cave"] +RPCstationsPos_node = ["volRpc","volMagneticSpectrometer","volIron","volHPT"] + +print "OPERApos: ",RPCstationsPos +print "trackingPos: ",trackStationsPos +print "strawVetoPos: ",strawVetoPos +print "scintPos: ",vetoWall +print "upstreamVetoPos: ",upstreamVetoPos + + +### Ntuple elements: +nf = TFile(outputFileName,"RECREATE") + +nt = TTree("t","t") +nt, event = addVar(nt, 'event','i') +nt, nNu = addVar(nt, 'nNu', 'i') +nt, isPrimary = addVar(nt, 'isPrimary','i') +nt, startZ_nu = addVar(nt, 'startZ_nu', 'f') +nt, startY_nu = addVar(nt, 'startY_nu', 'f') +nt, startX_nu = addVar(nt, 'startX_nu', 'f') +nt, nuE = addVar(nt, 'nuE', 'f') +nt, weight = addVar(nt, 'weight', 'f') +#nt, hansW = addVar(nt, 'hansW', 'f') + +## 0: last passive material opera-mu-system +## 1: two windows around liquid scintillator +## : along vaccum tank +## 3: along all OPERA +## 4: between the two entrance windows +## 5: along vacuum tank outer window +## 6: along vacuum tank inner window +## 999: wrong OPERA place ## needed until we have the new production +## -1: anything else, but what???? #it should never been present + +nt, interactionElement = addVar(nt, 'interactionElement','i') +### From Neutrino +nt, nPart_fromNu = addVar(nt, 'nPart_fromNu', 'i') +nt, idPart_fromNu = addVect(nt, 'idPart_fromNu', 'int') +nt, idStrPart_fromNu = addVect(nt, 'idStrPart_fromNu', 'string') + +nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i') +nt, idChargedPart_fromNu = addVect(nt, 'idChargedPart_fromNu', 'int') +nt, idStrChargedPart_fromNu = addVect(nt, 'idStrChargedPart_fromNu', 'string') + +nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i') +nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int') +nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string') + +''' +nt, rpcPoint_x = addVect(nt, 'rpcPoint_x', 'float') +nt, rpcPoint_y = addVect(nt, 'rpcPoint_y', 'float') +nt, rpcPoint_z = addVect(nt, 'rpcPoint_z', 'float') + +nt, strawPoint_x = addVect(nt, 'strawPoint_x', 'float') +nt, strawPoint_y = addVect(nt, 'strawPoint_y', 'float') +nt, strawPoint_z = addVect(nt, 'strawPoint_z', 'float') + +nt, strawVetoPoint_x = addVect(nt, 'strawVetoPoint_x', 'float') +nt, strawVetoPoint_y = addVect(nt, 'strawVetoPoint_y', 'float') +nt, strawVetoPoint_z = addVect(nt, 'strawVetoPoint_z', 'float') + + +nt, scintWPoint_x = addVect(nt, 'upstreamVetoPoint_x', 'float') +nt, scintWPoint_y = addVect(nt, 'upstreamVetoPoint_y', 'float') +nt, scintWPoint_z = addVect(nt, 'upstreamVetoPoint_z', 'float') +''' +### RPC ### +# was at least an RPC fired by something coming from the nu interaction? +nt, RPC = addVar(nt,'RPC', 'i' ) +# accounting for an eff. of each station of 90% +nt, RPC_eff = addVar(nt,'RPC_eff', 'i' ) + +#nt, nRPCstations = addVar(nt, 'nRPCstations', 'i') +#nt, nRPCstations_eff = addVar(nt, 'nRPCstations_eff', 'i') + +## In case something (no matter if it comes from the nu-interaction kids) fired the RPC +#nt, nRPCstationsAny = addVar(nt, 'nRPCstationsAny', 'i') +nt, RPCany = addVar(nt,'RPCany', 'i' ) +#nt, nRPCstationsAny_eff = addVar(nt, 'nRPCstationsAny_eff', 'i') +nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' ) +############ + +### scintW ### +# was at least an scintW fired by something coming from the nu interaction? +nt, scintW = addVar(nt,'upstreamVeto', 'i' ) +nt, scintWany = addVar(nt,'upstreamVetoAny', 'i' ) + +# accounting for an eff. of each station of 90% +nt, scintW_eff = addVar(nt,'upstreamVeto_eff', 'i' ) +nt, scintWany_eff = addVar(nt,'upstreamVetoAny_eff', 'i' ) + +# was at least an scintW fired by something coming from the nu interaction? +nt, scintVeto = addVar(nt,'scintVeto', 'i' ) +# accounting for an eff. of each station of 90% +nt, scintVeto_eff = addVar(nt,'scintVeto_eff', 'i' ) + + +## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls +#nt, nScintVetoStationsAny = addVar(nt, 'nScintVetoStationsAny', 'i') +nt, scintVetoAny = addVar(nt,'scintVetoAny', 'i' ) +nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' ) + +## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls +nt, scintVetoAny = addVar(nt,'upstreamVetoAny', 'i' ) +nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' ) + +nt, hitCharges = addVect(nt, "hitCharges", 'int') +nt, hitTrackId = addVect(nt, "hitTrackId", 'int') + +############### + +### strawVeto ### +# was at least an strawVeto fired by something coming from the nu interaction? +nt, strawVeto = addVar(nt,'strawVeto', 'i' ) +# accounting for an eff. of each station of 90% +nt, strawVeto_eff = addVar(nt,'strawVeto_eff', 'i' ) + + +## In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto +#nt, nStrawVetoStationsAny = addVar(nt, 'nStrawVetoStationsAny', 'i') +nt, strawVetoAny = addVar(nt,'strawVetoAny', 'i' ) +nt, strawVetoAny_eff = addVar(nt,'strawVetoAny_eff', 'i' ) +############### + + +### TrackSyst ### +## only the case of anything fired it is considered, obviously +nt, TrackSyst = addVar(nt, "TrackSyst", 'i') +nt, TrackSyst_eff = addVar(nt, "TrackSyst_eff", 'i') + +################## + + +nt, decayPartIndex = addVect(nt, "decayPartIndex", "float") +nt, decayPartID = addVect(nt, "decayPartID", "float") +nt, decayPartStrID = addVect(nt, "decayPartStrID", "string") +nt, decayPos_x = addVect(nt, "decayPos_x", "float") +nt, decayPos_y = addVect(nt, "decayPos_y", "float") +nt, decayPos_z = addVect(nt, "decayPos_z", "float") +nt, decayMother = addVect(nt, "decayMother", "float") +nt, decayP = addVect(nt,"decayP","float") + +nt, strawVeto_Ethr = addVar(nt,"strawVeto_Ethr","i") +nt, scintVeto_Ethr = addVar(nt,"scintVeto_Ethr","i") +nt, scintW_Ethr = addVar(nt,"upstreamVeto_Ethr","i") +nt, RPC_Ethr = addVar(nt,"RPC_Ethr","i") +nt, strawVetoAny_Ethr = addVar(nt,"strawVetoAny_Ethr","i") +nt, scintVetoAny_Ethr = addVar(nt,"scintVetoAny_Ethr","i") +nt, scintWAny_Ethr = addVar(nt,"upstreamVetoAny_Ethr","i") +nt, RPCany_Ethr = addVar(nt,"RPCany_Ethr","i") +nt, TrackSyst_Ethr = addVar(nt, "TrackSyst_Ethr", "i") + +nt, NC = addVar(nt, "NC", "i") + +#nt, nParticles = addVar(nt, "nParticles", "i") + +nt, nuInteractionNode = addVect(nt, "nuInteractionNode", "string") +#nt, nuInteractionNode_sel = addVar(nt, "nuInteractionNode_sel","i") +nt, nuInteractionNodeSimpl = addVect(nt, "nuInteractionNodeSimpl", "string") + +nt, nuIntNumSimpl = addVar(nt,"nuIntNumSimpl","i") + +nt, TrackSyst_node = addVar(nt, "TrackSyst_node", "i") +nt, strawVetoAny_node = addVar(nt, "strawVetoAny_node", "i") +nt, scintWany_node = addVar(nt, "scintWany_node", "i") +nt, RPCany_node = addVar(nt, "RPCany_node", "i") +nt, scintVetoAny_node = addVar(nt, "scintVetoAny_node", "i") + +nt, recoed = addVar(nt, "recoed","i") +nt, nRecoed = addVar(nt, "nRecoed","i") + +dictNodeNames = {'volIron':0, 'cave':1, 'LiSc':2, 'Startplate':3, 'TI':4, 'rockD':5, 'Endplate':6, 'Rib':7, 'volFeYoke':8, 'volHPT':9, 'TO':10, 'volRpc':11, 'volCoil':12, 'T1Lid':13, 'straw':14, 'strawVeto':15, 'volBase':16, 'Tr':17, 'gas':18, 'wire':19, 'others':20} + +entries = 10000 +for entry in xrange(entries): + if not (entry%1000): + print "%s / %s"%(entry,entries) + + t.GetEntry(entry) + + event[0] = entry + particles = t.MCTrack + rpc = t.ShipRpcPoint + scint = t.vetoPoint + strawPoints = t.strawtubesPoint + vetoScint = t.vetoPoint + recoParts = t.Particles + + ## initialization + putToZero(nNu) + putToZero(nPart_fromNu) + putToZero(nChargedPart_fromNu) + putToZero(nNeutrPart_fromNu) + + '''putToZero(nNeutrPart_withKs) + putToZero(nChargedPart_withKs) + putToZero(nPart_withKs) + + putToZero(nNeutrPart_withKL) + putToZero(nChargedPart_withKL) + putToZero(nPart_withKL) + + putToZero(nNeutrPart_withN) + putToZero(nChargedPart_withN) + putToZero(nPart_withN) + putToZero(hasKs) + putToZero(hasKL) + putToZero(hasN) + ''' + hitCharges.clear() + hitTrackId.clear() + + idPart_fromNu.clear() + idChargedPart_fromNu.clear() + idNeutrPart_fromNu.clear() + idStrPart_fromNu.clear() + idStrChargedPart_fromNu.clear() + idStrNeutrPart_fromNu.clear() + + ''' + strawVetoPoint_x.clear() + strawVetoPoint_y.clear() + strawVetoPoint_z.clear() + + strawPoint_x.clear() + strawPoint_y.clear() + strawPoint_z.clear() + + rpcPoint_x.clear() + rpcPoint_y.clear() + rpcPoint_z.clear() + + scintWPoint_x.clear() + scintWPoint_y.clear() + scintWPoint_z.clear() + ''' + decayPartIndex.clear() + decayPartID.clear() + decayPartStrID.clear() + decayPos_x.clear() + decayPos_y.clear() + decayPos_z.clear() + decayMother.clear() + decayP.clear() + + nuInteractionNode.clear() + + primaryDone=False + isPrimary[0] = int(False) + + lookingForDecay(particles) + nRecoed[0] = 0 + recoed[0] = 0 + #for (ip,part) in enumerate(particles): + + if sampleType=="nuBg" or sampleType == "cosmics": + starPartRange = [0] + elif sampleType=="sig": + starPartRange = [1] + + for ip in startPartRange: + + TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, hitCharges, hitTrackId, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] ) + if TrackSyst[0]==0: + skipEvent=True + continue + skipEvent=False + + part = particles[ip] + + pdgPart = pdg.GetParticle(part.GetPdgCode()) + if sampleType=="nuBg": + assert("nu" in pdgPart.GetName()) + + pushOfflineByEvent(t, sampleType) + ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother + #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena + if True: + ## Starting the counter of how many particles were produced by the interaction of this specific nu + #fromThisNu.push_back(0) + for recoP in recoParts: + nRecoed[0] += 1 + pushOfflineByParticle(t, recoP) + + if nRecoed[0]>0: + recoed[0]=1 + #else: + #skipEvent=True + #continue + + nNu[0]+=1 + if part.GetMotherId()==-1: + assert(primaryDone==False) + primaryDone=True + isPrimary[0] = int(True) + assert(len(idStrPart_fromNu)==0) + assert(len(idStrNeutrPart_fromNu)==0) + else: + continue + + + tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName() + nuInteractionNode.push_back(tmpName) + + tmpName = convertion(tmpName) + + nuInteractionNodeSimpl.push_back(tmpName) + if not tmpName in dictNodeNames: tmpName = "others" + + nuIntNumSimpl[0] = dictNodeNames[tmpName] + + startZ_nu[0] = part.GetStartZ() + startY_nu[0] = part.GetStartY() + startX_nu[0] = part.GetStartX() + + nuE[0] = part.GetEnergy() + # Re-initialization for each interacting neutrino + ## ---> <--- ## + + ## Looking for particles produced by the neutrino interaction + interacted = False + partKidsId = [] + #indexPartKids = [] + + #maxPcharged[0] = 0 + #maxPneutr[0] = 0 + NC[0] = 0 + for ip2 in xrange(0,len(particles)): + part2 = particles[ip2] + ## exit if we have reached the empty part of the array + if not (type(part2)==type(ShipMCTrack())): + break + if part2.GetMotherId()==ip: + interacted = True + part2Id = part2.GetPdgCode() + partKidsId.append(part2Id) + #indexPartKids.append(ip2) + + if not pdg.GetParticle(part2.GetPdgCode()) == None and "nu" in pdg.GetParticle(part2.GetPdgCode()).GetName() and part2.GetMotherId()==0: + NC[0] = 1 + # kid_decay_x.push_back(part2.GetX()) + #kid_decay_y.push_back(part2.GetY()) + #kid_decay_z.push_back(part2.GetZ()) + + + tmpP = part2.GetP() + #print "pID",part2Id, pdg.GetParticle(part2Id).GetName() + + nPart_fromNu[0]+=1 + idPart_fromNu.push_back(part2Id) + if pdg.GetParticle(part2.GetPdgCode()) == None: + part2Name = str(part2Id) + else: + part2Name = getPartName(part2Id) + idStrPart_fromNu.push_back(part2Name) + + ## Counting how many particles produced by the nu-interaction are charged or neutral + if pdg.GetParticle(part2.GetPdgCode()) == None or int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0): + nNeutrPart_fromNu[0]+=1 + idNeutrPart_fromNu.push_back(part2Id) + idStrNeutrPart_fromNu.push_back(part2Name) + #if fabs(maxPcharged[0])lastPassiveEl[0] and part2Z ",nodeName, lastPassive_nodeName + + if nodeName in lastPassive_nodeName: + somewhere = True + interactionElement[0] = 0 + intElName = 'OPERA' + #RpcPassiveFlag = True + #lastPassive[0] = int(True) + + ## To know if it was in the full OPERA-system (excluded last passive) + #if part2Z>OPERA[0] and part2ZOPERA_wrong[0] and part2ZentranceWindows[0][1] and part2Z=-2501.8 && startZ_nu<=-2478.8") +''' + + +'''### Looking for the first RPC station that can be fired (i.e. after the interaction point) + leftRPC = None#tmp_firstRPCStation = None + rightRPC = None + for (indexRPC,RPCs) in enumerate(RPCstationsPos): + if part2ZRPCstationsPos[potRPCindex][0]: + RPCflag = True + + ## For each station after the interaction point I simulate a 90% efficiency of the detector + ## only if one station is still fired I put a positive flag + eff_flag = gRandom.Uniform(0.,1.) + if eff_flag<0.9: + RPCflag_eff = True + print RPCflag, RPCflag_eff + + for (pid,part) in enumerate(particles): + print pid, pdg.GetParticle(part.GetName(), part.GetMotherId() + assert(False) +''' diff --git a/offlineForBarbara.py b/offlineForBarbara.py new file mode 100644 index 0000000..4415784 --- /dev/null +++ b/offlineForBarbara.py @@ -0,0 +1,432 @@ +from lookAtGeo import * +import tools +import shipunit as u +from ShipGeoConfig import ConfigRegistry +from operator import mul, add + +dy = 10. +# init geometry and mag. field +#ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) + + +def searchForNodes3_xyz_dict(fGeo, verbose=False): + from tools import findPositionElement, findDimentionBoxElement, findPositionElement2 + d = {} + #r = loadGeometry(inputFile) + #fGeo = r['fGeo'] + ## Get the top volume + #fGeo = ROOT.gGeoManager + tv = fGeo.GetTopVolume() + topnodes = tv.GetNodes() + for (j,topn) in enumerate(topnodes): + # top volumes + if verbose: + print j, topn.GetName() + print " x: ", findPositionElement(topn)['x'],findDimentionBoxElement(topn)['x'] + print " y: ", findPositionElement(topn)['y'],findDimentionBoxElement(topn)['y'] + print " z: ", findPositionElement(topn)['z'],findDimentionBoxElement(topn)['z'] + d[topn.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}} + d[topn.GetName()]['x']['pos'] = findPositionElement(topn)['x'] + d[topn.GetName()]['x']['dim'] = findDimentionBoxElement(topn)['x'] + d[topn.GetName()]['y']['pos'] = findPositionElement(topn)['y'] + d[topn.GetName()]['y']['dim'] = findDimentionBoxElement(topn)['y'] + d[topn.GetName()]['z']['pos'] = findPositionElement(topn)['z'] + d[topn.GetName()]['z']['dim'] = findDimentionBoxElement(topn)['z'] + if topn.GetVolume().GetShape().IsCylType(): + d[topn.GetName()]['r']['pos'] = findPositionElement(topn)['r'] + d[topn.GetName()]['r']['dim'] = findDimentionBoxElement(topn)['r'] + else: + d[topn.GetName()]['r']['pos'] = 0. + d[topn.GetName()]['r']['dim'] = 0. + # First children + nodes = topn.GetNodes() + if nodes: + topPos = topn.GetMatrix().GetTranslation() + for (i,n) in enumerate(nodes): + if verbose: + print j, topn.GetName(), i, n.GetName() + print " x: ", findPositionElement2(n,topPos)['x'],findDimentionBoxElement(n)['x'] + print " y: ", findPositionElement2(n,topPos)['y'],findDimentionBoxElement(n)['y'] + print " z: ", findPositionElement2(n,topPos)['z'],findDimentionBoxElement(n)['z'] + d[n.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}} + d[n.GetName()]['x']['pos'] = findPositionElement2(n,topPos)['x'] + d[n.GetName()]['x']['dim'] = findDimentionBoxElement(n)['x'] + d[n.GetName()]['y']['pos'] = findPositionElement2(n,topPos)['y'] + d[n.GetName()]['y']['dim'] = findDimentionBoxElement(n)['y'] + d[n.GetName()]['z']['pos'] = findPositionElement2(n,topPos)['z'] + d[n.GetName()]['z']['dim'] = findDimentionBoxElement(n)['z'] + if n.GetVolume().GetShape().IsCylType(): + d[n.GetName()]['r']['pos'] = findPositionElement2(n,topPos)['r'] + d[n.GetName()]['r']['dim'] = findDimentionBoxElement(n)['r'] + else: + d[n.GetName()]['r']['pos'] = 0. + d[n.GetName()]['r']['dim'] = 0. + # Second children + cnodes = n.GetNodes() + if cnodes: + localpos = n.GetMatrix().GetTranslation() + localToGlobal = [] + for i in xrange(3): + localToGlobal.append(localpos[i]+topPos[i]) + for (k,cn) in enumerate(cnodes): + if verbose: + print j, topn.GetName(), i, n.GetName(), k, cn.GetName() + print " x: ", findPositionElement2(cn,localToGlobal)['x'],findDimentionBoxElement(cn)['x'] + print " y: ", findPositionElement2(cn,localToGlobal)['y'],findDimentionBoxElement(cn)['y'] + print " z: ", findPositionElement2(cn,localToGlobal)['z'],findDimentionBoxElement(cn)['z'] + d[cn.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}} + d[cn.GetName()]['x']['pos'] = findPositionElement2(cn,localToGlobal)['x'] + d[cn.GetName()]['x']['dim'] = findDimentionBoxElement(cn)['x'] + d[cn.GetName()]['y']['pos'] = findPositionElement2(cn,localToGlobal)['y'] + d[cn.GetName()]['y']['dim'] = findDimentionBoxElement(cn)['y'] + d[cn.GetName()]['z']['pos'] = findPositionElement2(cn,localToGlobal)['z'] + d[cn.GetName()]['z']['dim'] = findDimentionBoxElement(cn)['z'] + if cn.GetVolume().GetShape().IsCylType(): + d[cn.GetName()]['r']['pos'] = findPositionElement2(cn,localToGlobal)['r'] + d[cn.GetName()]['r']['dim'] = findDimentionBoxElement(cn)['r'] + else: + d[cn.GetName()]['r']['pos'] = 0. + d[cn.GetName()]['r']['dim'] = 0. + return d + + +ff = ROOT.TFile("histoForWeights.root") +h_GioHans = ff.Get("h_Gio") +def calcWeightNu(NC, E, w, entries, nuName, ON=True): + # Only for neutrinos and antineutrinos + if not ON: + return 1 + if "bar" in nuName: + reduction = 0.5 + flux = 6.98e+11 * 2.e+20 / 5.e+13 + else: + reduction = 1. + flux = 1.09e+12 * 2.e+20/ 5.e+13 + crossSec = 6.7e-39*E * reduction + NA = 6.022e+23 + binN = h_GioHans.GetXaxis().FindBin(E) + return crossSec * flux * h_GioHans.GetBinContent(binN) * w * NA / entries + + +def findWeight(sampleType, NC, E, MCTrack, entries, nuName, ON): + if sampleType == 'nuBg': + return calcWeightNu(NC, E, MCTrack.GetWeight(), entries, nuName, ON) + elif sampleType == 'sig': + return MCTrack.GetWeight() # for the acceptance, multiply by normalization + elif sampleType == 'cosmics': + return MCTrack.GetWeight() # multiply by 1.e6 / 200. + + + +def hasStrawStations(event, trackId, listOfWantedStraws): + ok = [False]*len(listOfWantedStraws) + positions = [ (nodes[det]['z']['pos'] - nodes[det]['z']['dim'], nodes[det]['z']['pos'] + nodes[det]['z']['dim'] ) for det in listOfWantedStraws ] + for hit in event.strawtubesPoint: + if hit.GetTrackID() == trackId: + for (i,det) in enumerate(listOfWantedStraws): + if (positions[i][0] < hit.GetZ() < positions[i][1]) and tools.inEllipse(hit.GetX(), hit.GetY(), 250., 500.): + ok[i] = True + return bool(reduce(mul, ok, 1)) + +def hasGoodStrawStations(event, trackId): + #ok = [False]*2 + okupstream = [False]*2 + okdownstream = [False]*2 + upstream = [ (nodes[det]['z']['pos'] - nodes[det]['z']['dim'], nodes[det]['z']['pos'] + nodes[det]['z']['dim'] ) for det in ['Tr1_1', 'Tr2_2'] ] + downstream = [ (nodes[det]['z']['pos'] - nodes[det]['z']['dim'], nodes[det]['z']['pos'] + nodes[det]['z']['dim'] ) for det in ['Tr3_3', 'Tr4_4'] ] + for hit in event.strawtubesPoint: + if hit.GetTrackID() == trackId: + for i in xrange(2): + if (upstream[i][0] < hit.GetZ() < upstream[i][1]) and tools.inEllipse(hit.GetX(), hit.GetY(), 250., 500.): + okupstream[i] = True + if (downstream[i][0] < hit.GetZ() < downstream[i][1]) and tools.inEllipse(hit.GetX(), hit.GetY(), 250., 500.): + okdownstream[i] = True + ok = [ bool(reduce(mul, l, 1)) for l in [okupstream, okdownstream] ] + return bool(reduce(add, ok, 0)) + +def findHNLvertex(event): + for t in event.MCTrack: + if t.GetMotherId() == 1: + return t.GetStartZ() + return False + +def has_muon_station(event, trackId, station): + zIn = nodes['muondet%s_1'%(station-1)]['z']['pos'] - nodes['muondet%s_1'%(station-1)]['z']['dim'] + zOut = nodes['muondet%s_1'%(station-1)]['z']['pos'] + nodes['muondet%s_1'%(station-1)]['z']['dim'] + for hit in event.muonPoint: + if hit.GetTrackID() == trackId: + if zIn <= hit.GetZ() <= zOut: + return True + return False + +def hasEcalDeposit(event, trackId, ELossThreshold): + ELoss = 0. + for hit in event.EcalPoint: + if hit.GetTrackID() == trackId: + ELoss += hit.GetEnergyLoss() + if ELoss >= ELossThreshold: + return True + return False + +def hasMuons(event, trackId): + m1 = 0 + m2 = 0 + m3 = 0 + m4 = 0 + for ahit in event.muonPoint: + if ahit.GetTrackID() == trackId: + detID = ahit.GetDetectorID() + if(detID == 476) : + m1 += 1 + if(detID == 477) : + m2 += 1 + if(detID == 478) : + m3 += 1 + if(detID == 479) : + m4 += 1 + return [bool(m1), bool(m2), bool(m3), bool(m4)] + +def myVertex(t1,t2,PosDir): + # closest distance between two tracks + # d = |pq . u x v|/|u x v| + a = ROOT.TVector3(PosDir[t1][0](0) ,PosDir[t1][0](1), PosDir[t1][0](2)) + u = ROOT.TVector3(PosDir[t1][1](0),PosDir[t1][1](1),PosDir[t1][1](2)) + c = ROOT.TVector3(PosDir[t2][0](0) ,PosDir[t2][0](1), PosDir[t2][0](2)) + v = ROOT.TVector3(PosDir[t2][1](0),PosDir[t2][1](1),PosDir[t2][1](2)) + pq = a-c + uCrossv = u.Cross(v) + dist = pq.Dot(uCrossv)/(uCrossv.Mag()+1E-8) + # u.a - u.c + s*|u|**2 - u.v*t = 0 + # v.a - v.c + s*v.u - t*|v|**2 = 0 + E = u.Dot(a) - u.Dot(c) + F = v.Dot(a) - v.Dot(c) + A,B = u.Mag2(), -u.Dot(v) + C,D = u.Dot(v), -v.Mag2() + t = -(C*E-A*F)/(B*C-A*D) + X = c.x()+v.x()*t + Y = c.y()+v.y()*t + Z = c.z()+v.z()*t + # sT = ROOT.gROOT.FindAnything('cbmsim') + #print 'test2 ',X,Y,Z,dist + #print 'truth',sTree.MCTrack[2].GetStartX(),sTree.MCTrack[2].GetStartY(),sTree.MCTrack[2].GetStartZ() + return X,Y,Z,abs(dist) + +def addFullInfoToTree(elenaTree): + elenaTree, DaughtersPt = tools.AddVect(elenaTree, 'DaughtersPt', 'float') + elenaTree, DaughtersChi2 = tools.AddVect(elenaTree, 'DaughtersChi2', 'float') + elenaTree, DaughtersNPoints = tools.AddVect(elenaTree, 'DaughtersNPoints', 'int') + elenaTree, DaughtersTruthProdX = tools.AddVect(elenaTree, 'DaughtersTruthProdX', 'float') + elenaTree, DaughtersTruthProdY = tools.AddVect(elenaTree, 'DaughtersTruthProdY', 'float') + elenaTree, DaughtersTruthProdZ = tools.AddVect(elenaTree, 'DaughtersTruthProdZ', 'float') + elenaTree, DaughtersTruthPDG = tools.AddVect(elenaTree, 'DaughtersTruthPDG', 'int') + elenaTree, DaughtersTruthMotherPDG = tools.AddVect(elenaTree, 'DaughtersTruthMotherPDG', 'int') + elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') + elenaTree, straw_x = tools.AddVect(elenaTree, 'straw_x', 'float') + elenaTree, straw_y = tools.AddVect(elenaTree, 'straw_y', 'float') + elenaTree, straw_z = tools.AddVect(elenaTree, 'straw_z', 'float') + elenaTree, muon_x = tools.AddVect(elenaTree, 'muon_x', 'float') + elenaTree, muon_y = tools.AddVect(elenaTree, 'muon_y', 'float') + elenaTree, muon_z = tools.AddVect(elenaTree, 'muon_z', 'float') + elenaTree, ecal_x = tools.AddVect(elenaTree, 'ecal_x', 'float') + elenaTree, ecal_y = tools.AddVect(elenaTree, 'ecal_y', 'float') + elenaTree, ecal_z = tools.AddVect(elenaTree, 'ecal_z', 'float') + elenaTree, hcal_x = tools.AddVect(elenaTree, 'hcal_x', 'float') + elenaTree, hcal_y = tools.AddVect(elenaTree, 'hcal_y', 'float') + elenaTree, hcal_z = tools.AddVect(elenaTree, 'hcal_z', 'float') + elenaTree, veto5_x = tools.AddVect(elenaTree, 'veto5_x', 'float') + elenaTree, veto5_y = tools.AddVect(elenaTree, 'veto5_y', 'float') + elenaTree, veto5_z = tools.AddVect(elenaTree, 'veto5_z', 'float') + elenaTree, liquidscint_x = tools.AddVect(elenaTree, 'liquidscint_x', 'float') + elenaTree, liquidscint_y = tools.AddVect(elenaTree, 'liquidscint_y', 'float') + elenaTree, liquidscint_z = tools.AddVect(elenaTree, 'liquidscint_z', 'float') + elenaTree, DOCA = tools.AddVar(elenaTree, 'DOCA', 'float') + elenaTree, vtxx = tools.AddVar(elenaTree, 'vtxx', 'float') + elenaTree, vtxy = tools.AddVar(elenaTree, 'vtxy', 'float') + elenaTree, vtxz = tools.AddVar(elenaTree, 'vtxz', 'float') + elenaTree, IP0 = tools.AddVar(elenaTree, 'IP0', 'float') + elenaTree, Mass = tools.AddVar(elenaTree, 'Mass', 'float') + elenaTree, Pt = tools.AddVar(elenaTree, 'Pt', 'float') + elenaTree, P = tools.AddVar(elenaTree, 'P', 'float') + elenaTree, NParticles = tools.AddVar(elenaTree, 'NParticles', 'int') + elenaTree, HNLw = tools.AddVar(elenaTree, 'HNLw', 'float') + elenaTree, NuWeight = tools.AddVar(elenaTree, 'NuWeight', 'float') + elenaTree, EventNumber = tools.AddVar(elenaTree, 'EventNumber', 'int') + elenaTree, DaughterMinPt = tools.AddVar(elenaTree, 'DaughterMinPt', 'float') + elenaTree, DaughterMinP = tools.AddVar(elenaTree, 'DaughterMinP', 'float') + elenaTree, DaughtersAlwaysIn = tools.AddVar(elenaTree, 'DaughtersAlwaysIn', 'int') + elenaTree, BadTruthVtx = tools.AddVar(elenaTree, 'BadTruthVtx', 'int') + +DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal = None, None, None, None, None, None, None +DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2 = None, None, None, None, None, None +MaxDaughtersRedChi2, MinDaughtersNdf = None, None + +def addOfflineToTree(elenaTree): + global DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal + global DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2 + global MaxDaughtersRedChi2, MinDaughtersNdf, HNLw, NuWeight + elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') # + elenaTree, DOCA = tools.AddVect(elenaTree, 'DOCA', 'float') # + elenaTree, vtxx = tools.AddVect(elenaTree, 'vtxx', 'float') # + elenaTree, vtxy = tools.AddVect(elenaTree, 'vtxy', 'float') # + elenaTree, vtxz = tools.AddVect(elenaTree, 'vtxz', 'float') # + elenaTree, IP0 = tools.AddVect(elenaTree, 'IP0', 'float') # + #elenaTree, NParticles = tools.AddVar(elenaTree, 'NParticles', 'int') # + elenaTree, DaughtersAlwaysIn = tools.AddVect(elenaTree, 'DaughtersAlwaysIn', 'int') # + elenaTree, BadTruthVtx = tools.AddVect(elenaTree, 'BadTruthVtx', 'int') # + elenaTree, Has1Muon1 = tools.AddVect(elenaTree, 'Has1Muon1', 'int') # + elenaTree, Has1Muon2 = tools.AddVect(elenaTree, 'Has1Muon2', 'int') # + elenaTree, Has2Muon1 = tools.AddVect(elenaTree, 'Has2Muon1', 'int') # + elenaTree, Has2Muon2 = tools.AddVect(elenaTree, 'Has2Muon2', 'int') # + elenaTree, HasEcal = tools.AddVect(elenaTree, 'HasEcal', 'int') # + elenaTree, MaxDaughtersRedChi2 = tools.AddVect(elenaTree, 'MaxDaughtersRedChi2', 'float') # + elenaTree, MinDaughtersNdf = tools.AddVect(elenaTree, 'MinDaughtersNdf', 'int') # + +nodes = None +def loadNodes(fGeo): + global nodes + nodes = searchForNodes3_xyz_dict(fGeo) + +num_bad_z = 0 + +def signalNormalisationZ(tree, datatype, verbose): + # By event + # Uses MC truth!! + global BadTruthVtx, num_bad_z + tools.PutToZero(BadTruthVtx) + z_hnl_vtx = findHNLvertex(tree) + bad_z = False + if not z_hnl_vtx: + if "sig" in datatype: + print 'ERROR: hnl vertex not found!' + ii = 0 + for g in tree.MCTrack: + ii +=1 + if ("sig" in datatype) and ii < 3: + pass + elif ("sig" in datatype) and ii >= 3: + sys.exit() + if not (nodes['Veto_5']['z']['pos']-nodes['Veto_5']['z']['dim']-500. < z_hnl_vtx < nodes['Tr4_4']['z']['pos']+nodes['Tr4_4']['z']['dim']): + bad_z = True + num_bad_z += 1 + if "sig" in datatype: + print z_hnl_vtx + tools.Push(BadTruthVtx, int(bad_z)) + +def nParticles(tree): + # By event + global NParticles + np = 0 + for HNL in tree.Particles: + np += 1 + tools.PutToZero(NParticles); tools.Push(NParticles, np) + +def hasEcalAndMuons(tree, particle): + # By particle + global Has1Muon1, Has1Muon2, Has2Muon1 + global Has2Muon2, HasEcal + flag2Muon1 = False + flag2Muon2 = False + flag1Muon1 = False + flag1Muon2 = False + flagEcal = False + t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) + # AND or OR? + if ( has_muon_station(tree, t1, 1) and has_muon_station(tree, t2, 1) ): + flag2Muon1 = True + if ( has_muon_station(tree, t1, 2) and has_muon_station(tree, t2, 2) ): + flag2Muon2 = True + if ( has_muon_station(tree, t1, 1) or has_muon_station(tree, t2, 1) ): + flag1Muon1 = True + if ( has_muon_station(tree, t1, 2) or has_muon_station(tree, t2, 2) ): + flag1Muon2 = True + if ( hasEcalDeposit(tree, t1, 150.*u.MeV) or hasEcalDeposit(tree, t2, 150.*u.MeV) ): + flagEcal = True + tools.Push(Has2Muon1, int(flag2Muon1)) + tools.Push(Has2Muon2, int(flag2Muon2)) + tools.Push(Has1Muon1, int(flag1Muon1)) + tools.Push(Has1Muon2, int(flag1Muon2)) + tools.Push(HasEcal, int(flagEcal)) + +def chi2Ndf(tree, particle): + # By particle + global MaxDaughtersRedChi2, MinDaughtersNdf, DaughtersFitConverged + t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) + reducedChi2 = [] + ndfs = [] + converged = [] + for tr in [t1,t2]: + x = tree.FitTracks[tr] + ndfs.append( int(round(x.getFitStatus().getNdf())) ) + reducedChi2.append( x.getFitStatus().getChi2()/x.getFitStatus().getNdf() ) + converged.append( x.getFitStatus().isFitConverged() ) + tools.Push(MaxDaughtersRedChi2, max(reducedChi2)) + tools.Push(MinDaughtersNdf, min(ndfs)) + tools.Push( DaughtersFitConverged, int(converged[0]*converged[1]) ) + + +def goodBehavedTracks(tree, particle): + # By particle + # Uses MC truth!! + global DaughtersAlwaysIn + t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) + accFlag = True + for tr in [t1,t2]: + mctrid = tree.fitTrack2MC[tr] + if not hasGoodStrawStations(tree, mctrid): + accFlag = False + tools.Push(DaughtersAlwaysIn, int(accFlag)) + +def vertexInfo(tree, particle): + # By particle + global vtxx, vtxy, vtxz + global IP0, DOCA + t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) + PosDir = {} + for tr in [t1,t2]: + xx = tree.FitTracks[tr].getFittedState() + PosDir[tr] = [xx.getPos(),xx.getDir()] + xv,yv,zv,doca = myVertex(t1,t2,PosDir) + tools.Push(DOCA, doca) + tools.Push(vtxx, xv); tools.Push(vtxy, yv); tools.Push(vtxz, zv) + # impact parameter to target + HNLPos = ROOT.TLorentzVector() + particle.ProductionVertex(HNLPos) + HNLMom = ROOT.TLorentzVector() + particle.Momentum(HNLMom) + tr = ROOT.TVector3(0,0,ShipGeo.target.z0) + t = 0 + for i in range(3): t += HNLMom(i)/HNLMom.P()*(tr(i)-HNLPos(i)) + ip = 0 + for i in range(3): ip += (tr(i)-HNLPos(i)-t*HNLMom(i)/HNLMom.P())**2 + ip = ROOT.TMath.Sqrt(ip) + tools.Push(IP0, ip) + + +def prepareFillingsByParticle(): + # By event + global DaughtersAlwaysIn, DaughtersFitConverged, MinDaughtersNdf, MaxDaughtersRedChi2 + global Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2, HasEcal + global vtxx, vtxy, vtxz, IP0, DOCA + tools.PutToZero(DaughtersAlwaysIn) + tools.PutToZero(Has2Muon1); tools.PutToZero(Has2Muon2); tools.PutToZero(HasEcal) + tools.PutToZero(Has1Muon1); tools.PutToZero(Has1Muon2) + tools.PutToZero(DOCA) + tools.PutToZero(vtxx); tools.PutToZero(vtxy); tools.PutToZero(vtxz) + tools.PutToZero(IP0) + tools.PutToZero(MinDaughtersNdf); tools.PutToZero(MaxDaughtersRedChi2) + tools.PutToZero(DaughtersFitConverged) + + +def pushOfflineByEvent(tree, datatype, verbose): + # True HNL decay vertex (only for signal normalisation) + signalNormalisationZ(tree, datatype, verbose) + ## Number of particles + #nParticles(tree) + # Empties arrays filled by particle + prepareFillingsByParticle() + +def pushOfflineByParticle(tree, particle): + hasEcalAndMuons(tree, particle) + goodBehavedTracks(tree, particle) + chi2Ndf(tree, particle) + vertexInfo(tree, particle) \ No newline at end of file diff --git a/tools.py b/tools.py new file mode 100755 index 0000000..367e9f2 --- /dev/null +++ b/tools.py @@ -0,0 +1,281 @@ +#from elena import * +import math +from array import array +from collections import OrderedDict +import ROOT + +def findDimentionBoxElement(node): + sh = node.GetVolume().GetShape() + if sh.IsCylType(): + #print node.GetName(), " is a cylinder", sh.GetRmin(), sh.GetRmax() + if sh.HasRmin(): + r = (sh.GetRmax() - sh.GetRmin()) + else: + r = 0. + return {'x':sh.GetDX(), + 'y':sh.GetDY(), + 'z':sh.GetDZ(), + 'r':r} + return {'x':sh.GetDX(), + 'y':sh.GetDY(), + 'z':sh.GetDZ()} + +def findPositionElement(node): + pos = node.GetMatrix().GetTranslation() + sh = node.GetVolume().GetShape() + if sh.IsCylType(): + if sh.HasRmin(): + r = (sh.GetRmax() - sh.GetRmin())/2. + else: + r = sh.GetRmax() + return {'x':pos[0], + 'y':pos[1], + 'z':pos[2], + 'r':r} + return {'x':pos[0], + 'y':pos[1], + 'z':pos[2]} + +def findPositionElement2(node,topPos): + local_pos = node.GetMatrix().GetTranslation() + pos = [] + for i in xrange(3): + pos.append(local_pos[i]+topPos[i]) + sh = node.GetVolume().GetShape() + if sh.IsCylType(): + if sh.HasRmin(): + r = (sh.GetRmax() - sh.GetRmin())/2. + else: + r = sh.GetRmax() + return {'x':pos[0], + 'y':pos[1], + 'z':pos[2], + 'r':r} + return {'x':pos[0], + 'y':pos[1], + 'z':pos[2]} + +def findDimentionBoxElement2(node,top): + sh = node.GetVolume().GetShape() + if sh.IsCylType(): + #print node.GetName(), " is a cylinder", sh.GetRmin(), sh.GetRmax() + if sh.HasRmin(): + r = (sh.GetRmax() - sh.GetRmin()) + else: + r = 0. + return {'x':sh.GetDX(), + 'y':sh.GetDY(), + 'z':sh.GetDZ(), + 'r':r} + return {'x':sh.GetDX(), + 'y':sh.GetDY(), + 'z':sh.GetDZ()} + +# Custom sorting +def sortNodesBy(nodes, coord='z', what='pos', printdict=True): + d = OrderedDict(sorted(nodes.items(), key = lambda item: item[1][coord][what])) + if printdict: + print + print + print "Nodes by %s %s:"%(coord, what) + for node in d.keys(): + print node, '\t\t', (d[node][coord]['pos']-d[node][coord]['dim'])/100., ' m', '\t', 2*d[node][coord]['dim']/100., ' m' + return d + +def retrieveMCParticleInfo(sTree, key): + mcPartKey = sTree.fitTrack2MC[key] + mcPart = sTree.MCTrack[mcPartKey] + if not mcPart: + print "Error in retrieveMCParticleInfo" + sys.exit() + mother = sTree.MCTrack[mcPart.GetMotherId()] + #print mcPartKey, mcPart, mcPart.GetMotherId(), mother + try: + mumId = mother.GetPdgCode() + except: + mumId = 0 + return int(mcPart.GetPdgCode()), int(mumId), mcPart.GetStartX(), mcPart.GetStartY(), mcPart.GetStartZ() + + +def AddVect(t,name,vectType): + vect = ROOT.vector(vectType)() + t.Branch( name, vect ) + return t, vect + +def AddVar(t,name,varType): + var = array(varType[0],[-999]) + t.Branch(name,var,name+"/"+varType.upper()) + return t, var + +def PutToZero(var): + if not isinstance(var, array): + var.clear() + else: + var[0] = 0 + +def Push(leaf, value): + if not isinstance(leaf, array): + leaf.push_back(value) + else: + leaf[0] = value + + +def wasFired(indexKids, detPoints, detPos, pointsVects=None, Ethr=0): + def lookingForHits(detPoints,flag,nstat,nstat_eff,indexKids,Eloss,Ethr): + for hit in detPoints: + if (indexKids is None) or (hit.GetTrackID() in indexKids): + #print RPChit.GetTrackID(), t_ID + # check if it is in one of the considered active layer + for pos in detPos: + if pos[0]<=hit.GetZ()<=pos[1]: + Eloss += hit.GetEnergyLoss() + #print pos, hit.GetZ() + if pointsVects is not None: + pointsVects[0].push_back(hit.GetX()) + pointsVects[1].push_back(hit.GetY()) + pointsVects[2].push_back(hit.GetZ()) + + flag = True + nstat+=1 + #eff_val = gRandom.Uniform(0.,1.) + #if eff_val<0.9: + # flag_eff = flag_eff or True + # nstat_eff+=1 + #else: + # flag_eff = flag_eff or False + flag_Ethr = Eloss>=Ethr + return flag, nstat, nstat_eff, Eloss, flag_Ethr + # Now in partKidTrackID I should have the trackID connected to my charged particles + #flag_eff = False + flag = False + nstat = 0 + nstat_eff = 0 + Eloss = 0 + flag,nstat,nstat_eff,Eloss,flag_Ethr = lookingForHits(detPoints,flag,nstat,nstat_eff,indexKids,Eloss,Ethr) + #if flag==False and checkOn: + #print "To be study event %s"%entry + #PrintEventPart(particles,pdg) + #PrintRPChits(detPoints,pdg) + #print + #print + #assert(False) + #return flag, flag_eff, nstat, nstat_eff, flag_Ethr + return flag, flag_Ethr + + +#ff = ROOT.TFile("histoForWeights.root","read") +#h_GioHans = ff.Get("h_Gio") + +def calcWeight(E, w, entries, nuName, weightHist): + if "bar" in nuName: + reduction = 0.5 + flux = 6.98e+11 * 2.e+20 / 5.e+13 + else: + reduction = 1. + flux = 1.09e+12 * 2.e+20/ 5.e+13 + crossSec = 6.7e-39*E * reduction + NA = 6.022e+23 + binN = weightHist.GetXaxis().FindBin(E) + return crossSec * flux * weightHist.GetBinContent(binN) * w * NA / entries + + +def calcWeightOldNtuple(x,y,z, E, nodes, entries, weightHist, datatype): + if 'sig' in datatype: + return -999 + params = {'OPERA': {'material':'Fe','lenght':60.,}, + 'window1': {'material':'Fe','lenght':nodes["lidT1O_1"]['z']['dim']}, + 'window2': {'material':'Al','lenght':nodes["lidT1I_1"]['z']['dim']}, + 'volumeOut': {'material':'Fe','lenght':30.,}, + 'volumeIn': {'material':'Al','lenght':8.}} + materialPars = {'Fe':{'A':55.847, #g/mol # fuori + 'rho': 7874 * 1.e+3/1.e+6,#g/cm3 + }, + 'Al':{'A':26.98, #g/mol # dentro + 'rho': 2700 * 1.e+3/1.e+6 #g/cm3 + }} + perc = {'OPERA':0.33, + 'window1':0.015, + 'window2':0.015, + 'volumeOut':0.32, + 'volumeIn':0.32} + intElName = interactionElement(x,y,z, nodes) + #print intElName + NA = 6.022e+23 + material = params[intElName]['material'] + crossSec = 8.8e-39*E #1.5e-38*fabs(E) ##3./2*(1.2+0.35)*1.e-38*fabs(E) + flux = 2.e+18#5.e+16 + L = params[intElName]['lenght'] + rho = materialPars[material]['rho'] + n_element = perc[intElName]*entries + binN = weightHist.GetXaxis().FindBin(E) + weight = crossSec * NA * rho * L * flux / n_element * weightHist.GetBinContent(binN) + return weight + +def interactionElement(x,y,z, nodes): + # window1 + if nodes["lidT1O_1"]['z']['pos']-nodes["lidT1O_1"]['z']['dim'] < z < nodes["lidT1O_1"]['z']['pos']+nodes["lidT1O_1"]['z']['dim']: + return "window1" + # window2 + if nodes["lidT1I_1"]['z']['pos']-nodes["lidT1I_1"]['z']['dim'] < z < nodes["lidT1I_1"]['z']['pos']+nodes["lidT1I_1"]['z']['dim']: + return "window2" + # vacuum tank borders + if nodes["lidT1O_1"]['z']['pos']+nodes["lidT1O_1"]['z']['dim'] < z < nodes["lidT6O_1"]['z']['pos']+nodes["lidT6O_1"]['z']['dim']: # include le lid x la parte fuori ma non x la parte dentro + if inInnerLid(x,y,z, nodes): + return "volumeIn" + else: + return "volumeOut" + # opera + try: + minz, maxz = nodes["volIron_12"]['z']['pos']-nodes["volIron_12"]['z']['dim'], nodes["volIron_24"]['z']['pos']+nodes["volIron_24"]['z']['dim'] + except KeyError: + try: + minz, maxz = nodes["volLayer_11"]['z']['pos']-nodes["volLayer_11"]['z']['dim'], nodes["volLayer_0"]['z']['pos']+nodes["volLayer_0"]['z']['dim'] + except KeyError: + # Pazienza esaurita, ora lo hanno chiamato "volMagneticSpectrometer_1" ed e' lungo piu' di 9m!!! + try: + minz, maxz = nodes["volMagneticSpectrometer_1"]['z']['pos']-nodes["volMagneticSpectrometer_1"]['z']['dim'], nodes["volMagneticSpectrometer_1"]['z']['pos']+nodes["volMagneticSpectrometer_1"]['z']['dim'] + except KeyError: + print "Cannot find OPERA" + sys.exit() + #if nodes["volIron_12"]['z']['pos']-nodes["volIron_12"]['z']['dim'] < z < nodes["volIron_24"]['z']['pos']+nodes["volIron_24"]['z']['dim']: + #if nodes["volLayer_11"]['z']['pos']-nodes["volLayer_11"]['z']['dim'] < z < nodes["volLayer_0"]['z']['pos']+nodes["volLayer_0"]['z']['dim']: + if minz < z < maxz: + return "OPERA" + + +def inInnerLid(x,y,z, nodes): + if nodes["lidT1I_1"]['z']['pos']-nodes["lidT1I_1"]['z']['dim'] < z < nodes["lidT6I_1"]['z']['pos']+nodes["lidT6I_1"]['z']['dim']: + if z < nodes["Veto_5"]['z']['pos']-nodes["Veto_5"]['z']['dim']: + a,b = nodes['T1I_1']['x']['dim'], nodes['T1I_1']['y']['dim'] # smaller part of the tank (first 5 m) + else: + a,b = nodes['T2I_1']['x']['dim'], nodes['T2I_1']['y']['dim'] # larger part of the tank + if inEllipse(x,y,a,b): + return True + return False + + +def inEllipse(x,y,a,b): + if ( (x*x)/(a*a) + (y*y)/(b*b) ) < 1.: + return True + return False + +def pointInEllipse(point,a,b): + x = point.X() + y = point.Y() + return inEllipse(x,y,a,b) + +listOfSegments = [] + +def makeLSsegments(nodes): + zi = -2478 + zo = 3678 + zstepsize = 100 + Nz = int(math.ceil((zi-zo)/zstepsize)) + Nphi = 16 + dphi = 2*math.pi/N + phistart = dphi/2; + nBins = Nphi*Nz + SegmentIsHit = [0] * nBins + +def fillSegments(): + pass \ No newline at end of file