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
160
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
1.9.8