SND@LHC Software
Loading...
Searching...
No Matches
g4Ex_gap_mergeFiles.py
Go to the documentation of this file.
1from __future__ import print_function
2# result_1Bn_ecut_5.root 1E9 with Ecut > 5 GeV
3# result_0.1Bn_ecut_0.5.root 1E8 with Ecut > 0.5 GeV
4Yandex = False # False
5Yandex2 = False # summer 2015 production, 10B, Ecut>10GeV, Mo/W/ target
6Yandex3 = True # spring 2016 production, 10B, Ecut>10GeV, Mo/W/ target, record start of track
7afterHadronAbsorber = False
8JPsi = False # True
9MoTarget = False
10Tau = False
11
12import ROOT,os
13from array import array
14from ROOT import TDatabasePDG,TMath,gDirectory
15from rootUtils import *
16pdg = TDatabasePDG()
17mu = pdg.GetParticle(13)
18Mmu = mu.Mass()
19Mmu2 = Mmu * Mmu
20
21if Yandex:
22 stats = {5.:[1E9,1E9],0.5:[1E8]}
23 files = {5.:['$SHIPSOFT/data/result_1Bn_ecut_5.root','$SHIPSOFT/data/result_1Bn_ecut_5-v02.root'],0.5:['$SHIPSOFT/data/result_0.1Bn_ecut_0.5.root']}
24 fnew = 'pythia8_Geant4_total_Yandex.root'
25if Yandex2:
26 stats = {10.:[1E10]}
27 files = {10.:['$SHIPSOFT/data/result_10Bn_ecut_10.root']}
28 fnew = 'pythia8_Geant4_total_Yandex2.root'
29if Yandex3:
30 stats = {10.:[10.05E9]}
31 files = {10.:['/media/truf/Elements/Data/HNL/ShipSoft/data/Yandex-April2016_10.05B.root']}
32 fnew = 'pythia8_Geant4_total_Yandex3.root'
33if JPsi:
34 BR = 0.05961
35 stats = {10.:[]}
36 files = {10.:[]}
37 jobs = {'b6276c506a_g4Ex_gap_':['53','54','55','56'],'b602f257b0_g4Ex_gap_':['52']}
38 nevts = 0
39 ntot = 0
40 tag = '#events='
41 for job in jobs:
42 for run in jobs[job]:
43 for i in range(10):
44 path = job+run+str(i)+'/'
45 fl = open(path+'log'+run+str(i))
46 for l in fl.readlines():
47 k = l.find(tag)
48 if k<0: continue
49 nevts = int( l[k+len(tag):].replace(')',''))
50 fl.close()
51 files[10.].append(path+'pythia8_Geant4_'+run+str(i)+'_10.0.root')
52 stats[10.].append(nevts/BR)
53 ntot += nevts/BR
54 print(run+str(i),' --> nevts = ',nevts)
55 fnew = 'pythia8_Geant4_total_Jpsi.root'
56 print('total statistics ',ntot/1.E9,' *1E9')
57#
58if Tau:
59 BR = 0.0554
60 stats = {0.:[]}
61 files = {0.:[]}
62 jobs = {'b63edb1b12_g4Ex_gap_':['13']}
63 nevts = 0
64 ntot = 0
65 tag = '#events='
66 for job in jobs:
67 for run in jobs[job]:
68 for i in range(10):
69 path = job+run+str(i)+'/'
70 if not 'log'+run+str(i) in os.listdir(path): continue
71 fl = open(path+'log'+run+str(i))
72 for l in fl.readlines():
73 k = l.find(tag)
74 if k<0: continue
75 nevts = int( l[k+len(tag):].replace(')',''))
76 fl.close()
77 files[0.].append(path+'pythia8_Geant4_'+run+str(i)+'_0.0.root')
78 stats[0.].append(nevts/BR)
79 ntot += nevts/BR
80 print(run+str(i),' --> nevts = ',nevts)
81 fnew = 'pythia8_Geant4_total_tauOnly_MoTarget_E0.root'
82 print('total statistics ',ntot/1.E9,' *1E9')
83#
84if MoTarget:
85 stats = {}
86 files = {}
87 jobs = {'b64a4b817c.cern.ch_g4Ex_gap_':['60','61','62','65','66'],'b63edb1b12_g4Ex_gap_':['63','75','77','78','81','83','85','87','90','93','95'],\
88 'b65dfad864_g4Ex_gap_':['70'],'b6276c506a_g4Ex_gap_':['71'],'b6e72959f6_g4Ex_gap_':['72'],\
89 'b6429e4f38.cern.ch_g4Ex_gap_':['73','74','76','79','80','82','84','86','88','89','91','94','96'] }
90 nevts = 0
91 ntot = 0
92 tag = '#events='
93 for job in jobs:
94 for run in jobs[job]:
95 for i in range(10):
96 path = job+run+str(i)+'/'
97 fl = open(path+'log'+run+str(i))
98 nevts = 0
99 for l in fl.readlines():
100 k = l.find(tag)
101 if k<0: continue
102 nevts = int( l[k+len(tag):].replace(')',''))
103 fl.close()
104 if nevts==0: continue
105 for r in os.listdir(path):
106 if r.find('.root')<0: continue
107 ecut = float(r.split('.root')[0].split('_')[3])
108 if ecut not in stats:
109 stats[ecut] = []
110 files[ecut] = []
111 files[ecut].append(path+r)
112 stats[ecut].append(nevts)
113 ntot += nevts
114 print(run+str(i),' --> nevts = ',nevts)
115 fnew = 'pythia8_Geant4_total_MoTarget.root'
116 print('total statistics ',ntot/1.E9,' *1E9')
117
118ntot = {}
119for ecut in stats:
120 ntot[ecut] = 0
121 for s in stats[ecut]: ntot[ecut]+=s
122print(ntot)
123
124h={}
125def makeFinalNtuples(norm=5.E13,opt=''):
126 cuts = {'':'','_onlyMuons':'abs(id)==13','_onlyNeutrinos':'abs(id)==12||abs(id)==14||abs(id)==16'}
127 first = True
128 tuples = ''
129 fn = 1
130 for ecut in stats:
131 for i in range(len(stats[ecut])):
132 h[fn] = ROOT.TFile(files[ecut][i])
133 t = h[fn].FindObjectAny("pythia8-Geant4")
134 fn+=1
135 if first:
136 first = False
137 for l in t.GetListOfLeaves(): tuples += l.GetName()+':'
138 tuples+='w:ecut'
139 fxx = fnew.replace('.root',opt+'.root')
140 if opt!='': fxx = fxx.replace('_total','')
141 h['N'] = ROOT.TFile(fxx, 'RECREATE')
142 print('new file created',fxx)
143 if afterHadronAbsorber:
144 h['ntuple'] = ROOT.TNtuple("pythia8-Geant4","flux after 3m hadron absorber",tuples)
145 else:
146 h['ntuple'] = ROOT.TNtuple("pythia8-Geant4","flux after 3m hadron absorber, position&momentum before",tuples)
147 gROOT.cd()
148 t.SetEventList(0)
149 t.Draw(">>temp",cuts[opt])
150 temp = gROOT.FindObjectAny('temp')
151 t.SetEventList(temp)
152 nev = temp.GetN()
153 for iev in range(nev) :
154 rc = t.GetEntry(temp.GetEntry(iev))
155 leaves = t.GetListOfLeaves()
156 vlist = array('f')
157 for x in range(leaves.GetEntries()):
158 vlist.append( leaves.At(x).GetValue() )
159 # get kinetic energy "id:px:py:pz:x:y:z:ox:oy:oz:pythiaid:parentid"
160 # since 2016 "id:px:py:pz:x:y:z:opx:opy:opz:ox:oy:oz:pythiaid:parentid"
161 Psq = vlist[1]**2+vlist[2]**2+vlist[3]**2
162 if abs(vlist[0])==13: Ekin = ROOT.TMath.Sqrt(Mmu2+Psq)-Mmu
163 else: Ekin = ROOT.TMath.Sqrt(Psq)
164 if Yandex:
165 if Ekin < ecut : continue
166 if Ekin > 5. : weight = norm/(ntot[5.] + ntot[0.5])
167 else : weight = norm/(ntot[0.5])
168 if Yandex2 or Yandex3:
169 if Ekin < ecut : continue
170 weight = norm/(ntot[10.])
171 if JPsi : weight = norm/(ntot[10.])
172 if Tau : weight = norm/(ntot[0.])
173 if MoTarget:
174 if Ekin < ecut : continue
175 if Ekin > 50. : weight = norm/(ntot[0.5] + ntot[10.] + ntot[20.]+ ntot[50.])
176 elif Ekin > 20. : weight = norm/(ntot[0.5] + ntot[10.] + ntot[20.])
177 elif Ekin > 10. : weight = norm/(ntot[0.5] + ntot[10.])
178 else : weight = norm/(ntot[0.5])
179 vlist.append(weight)
180 vlist.append( float(ecut) )
181 h['ntuple'].Fill(vlist)
182 h['N'].cd()
183 h['ntuple'].Write()
184
185def interactionRegion():
186 import rootUtils as ut
187 import shipunit as u
188 f = ROOT.TFile('pythia8_Geant4_13_350.0.root')
189 sTree=f.FindObjectAny('pythia8-Geant4')
190 ut.bookHist(h,'originz','z',100,-50.5,-49.)
191 ut.bookHist(h,'originzr','r vs z',100,-50.5,-49.,100,0.,0.5)
192 ut.bookHist(h,'originxy','x vs y',100,-0.5,0.5,100,-0.5,0.5)
193 ROOT.gROOT.cd()
194 sTree.Draw('z>>originz')
195 sTree.Draw('1000*sqrt(x*x+y*y):z>>originzr')
196 sTree.Draw('1000*x:1000*y>>originxy')
197
198
199makeFinalNtuples(norm=5.E13,opt='')
200makeFinalNtuples(norm=5.E13,opt='_onlyMuons')
201makeFinalNtuples(norm=5.E13,opt='_onlyNeutrinos')
g4Ex_gap_mergeFiles.interactionRegion
interactionRegion()
Definition g4Ex_gap_mergeFiles.py:185
g4Ex_gap_mergeFiles.makeFinalNtuples
makeFinalNtuples(norm=5.E13, opt='')
Definition g4Ex_gap_mergeFiles.py:125