SND@LHC Software
Loading...
Searching...
No Matches
mergeMbias.py
Go to the documentation of this file.
1from __future__ import print_function
2from __future__ import division
3import ROOT,os,random
4import shipunit as u
5import rootUtils as ut
6from ShipGeoConfig import ConfigRegistry
7
8from array import array
9pdg = ROOT.TDatabasePDG()
10mu = pdg.GetParticle(13)
11Mmu = mu.Mass()
12Mmu2 = Mmu * Mmu
13rnr = ROOT.TRandom()
14eospath = ROOT.gSystem.Getenv("EOSSHIP")+"/eos/experiment/ship/data/"
15ship_geo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = 10.)
16endOfHadronAbsorber = (ship_geo['hadronAbsorber'].z + ship_geo['hadronAbsorber'].length/2.) /100.
17startOfTarget = -50. # value used for Geant4 production
18
19def fillPart(t):
20 particles = {}
21 for n in range(t.GetEntries()):
22 rc = t.GetEvent(n)
23 if t.parentid not in particles :
24 particles[t.parentid] = 0
25 particles[t.parentid] +=1
26 return particles
27#
28def fillWeights():
29 weights = {}
30 for p in productions:
31 f = ROOT.TFile.Open(eospath+productions[p]["file"])
32 t = f.FindObjectAny("pythia8-Geant4")
33 weights[p]={}
34 for n in range(t.GetEntries()):
35 rc = t.GetEvent(n)
36 if t.w not in weights[p] :
37 weights[p][t.w] = [t.ecut,0]
38 weights[p][t.w][1] +=1
39 if t.ecut > weights[p][t.w][0] : weights[p][t.w][0] = t.ecut
40 return weights
41
42def TplotP(sTree):
43 ut.bookCanvas(h,key='P',title='momentum',nx=1800,ny=1200,cx=3,cy=2)
44 ut.bookCanvas(h,key='>P',title='N >P', nx=1800,ny=1200,cx=3,cy=2)
45 ut.bookCanvas(h,key='PT',title='Pt',nx=1800,ny=1200,cx=3,cy=2)
46 ut.bookCanvas(h,key='>PT',title='N >Pt', nx=1800,ny=1200,cx=3,cy=2)
47 cuts = {'mu':'abs(id)==13','nu':'abs(id)!=13','mu-':'id==13','mu+':'id==-13',
48 'nutau':'id==16','nutaubar':'id==-16','numu':'id==14','numubar':'id==-14',
49 'nue':'id==12','nuebar':'id==-12','nuesum':'abs(id)==12','numusum':'abs(id)==14','nutausum':'abs(id)==16'}
50 OpenCharm = {'':'','charm':"&(pythiaid==id & (abs(parentid) == 15 || abs(parentid) == 4112 || abs(parentid) == 4122 || abs(parentid) == 4132 \
51 || abs(parentid) == 431 || abs(parentid) == 421 || abs(parentid) == 411) )"}
52 ROOT.gROOT.cd('')
53 for x in ['','charm']:
54 for q in cuts:
55 p = q+x
56 hn = 'Tp'+p
57 hnt = 'Tpt'+p
58 ut.bookHist(h,hn, p +' ;p [GeV] ;N',400,0.0,400.0)
59 ut.bookHist(h,hnt, p +' ;pt [GeV] ;N',40,0.0,4.0)
60 sTree.Draw('sqrt(px**2+py**2+pz**2)>>'+hn,'w*('+cuts[q]+OpenCharm[x]+')','goff')
61 sTree.Draw('sqrt(px**2+py**2)>>'+hnt,'w*('+cuts[q]+OpenCharm[x]+')','goff')
62 if q=='mu+': h[hn].SetLineColor(3)
63 if q=='mu-': h[hn].SetLineColor(4)
64# integrated rates
65 for q in [ 'mu','mu-','mu+','nu','nue','nuesum','nuebar','numusum','numu','numubar','nutau','nutaubar']:
66 for x in ['','charm']:
67 for z in ['p','pt']:
68 p = z+q+x
69 hi = 'T'+p+'_>E'
70 h[hi]=h['T'+p].Clone(hi)
71 h[hi].Reset()
72 nsum = 0
73 for i in range(h[hi].GetNbinsX()+1,0,-1):
74 nsum+=h['T'+p].GetBinContent(i)
75 h[hi].SetBinContent(i,nsum)
76 for z in ['p','pt']:
77 k = 1
78 for x in [ 'mu','nu']:
79 t=z.upper()
80 p=z+x
81 cv = h[t].cd(k)
82 cv.SetLogy(1)
83 h['T'+p].Draw()
84 if h['T'+p].GetEntries()<1: continue
85 if not p.find('mu')<0:
86 h['T'+p+'+'].Draw('same')
87 h['T'+p+'-'].Draw('same')
88 cv = h['>'+t].cd(k)
89 cv.SetLogy(1)
90 h[hi].Draw()
91 k+=1
92# plot different nu species:
93 k = 3
94 cv = h[t].cd(k)
95 cv.SetLogy(1)
96 first = True
97 i = 2
98 h['tlnu'+z+str(k)] = ROOT.TLegend(0.49,0.13,0.88,0.36)
99 for p in ['numu','numubar','nue','nuebar','nutau','nutaubar']:
100 hn = 'T'+z+p
101 h['tlnu'+z+str(k)].AddEntry(h[hn],z+' '+p,'PL')
102 h[hn].SetLineColor(i)
103 h[hn+'8']=h[hn].Clone()
104 h[hn+'8'].SetName(hn+'8')
105 h[hn+'8'].Rebin(8)
106 i+=1
107 if first:
108 h[hn+'8'].Draw()
109 first = False
110 h[hn+'8'].Draw('same')
111 h['tlnu'+z+str(k)].Draw()
112 k+=1
113#
114
115productions = {}
116allProds = False
117if allProds:
118 productions["CERN-Cracow"] = {"stats":{1.:[1.1E8],10.:[1.22E9],100:[1.27E10]},
119 "file":"Mbias/pythia8_Geant4_total.root" }
120# checked, 10 variables, parentid = 8
121 productions["Yandex"] = {"stats":{5.:[2.1E9,1E9],0.5:[1E8]},
122 "file":"Mbias/pythia8_Geant4_total_Yandex.root" }
123# checked, 13 variables, parentid = 11
124 productions["Yandex2"] = {"stats":{10.:[1E10]},
125 "file":"Mbias/pythia8_Geant4_total_Yandex2.root" }
126# now with mu momentum at prodcution, NOT after hadron absorber
127productions["Yandex3"] = {"stats":{10.:[1E10]},
128 "file":"Mbias/pythia8_Geant4_total_Yandex3.root" }
129# checked, 13 variables, parentid = 11
130
131fnew = "pythia8_Geant4-noOpenCharm.root"
132noOpCharm = "!(pythiaid==id & (abs(parentid) == 15 || abs(parentid) == 4112 || abs(parentid) == 4122 || abs(parentid) == 4132 \
133 || abs(parentid) == 431 || abs(parentid) == 421 || abs(parentid) == 411) )"
134OpCharm = "(pythiaid==id & (abs(parentid) == 15 || abs(parentid) == 4112 || abs(parentid) == 4122 || abs(parentid) == 4132 \
135 || abs(parentid) == 431 || abs(parentid) == 421 || abs(parentid) == 411) )"
136cuts = {'':'abs(id)>0','_onlyMuons':'abs(id)==13','_onlyNeutrinos':'abs(id)==12||abs(id)==14||abs(id)==16'}
137
138h={}
139
140def mergeMinBias(pot,norm=5.E13,opt=''):
141 storeCharm=False
142 if opt != '':
143 storeCharm=True
144 opt=''
145 first = True
146 for p in productions:
147 f = ROOT.TFile.Open(eospath+productions[p]["file"])
148 t = f.FindObjectAny("pythia8-Geant4")
149 if first:
150 first = False
151 tuples = ''
152 for l in t.GetListOfLeaves():
153 if tuples == '': tuples += l.GetName()
154 else: tuples += ":"+l.GetName()
155 fxx = fnew.replace('.root',opt+'.root')
156 if storeCharm: fxx = fnew.replace('.root','old-charm.root')
157 h['N'] = ROOT.TFile(fxx, 'RECREATE')
158 print('new file created',fxx)
159 h['ntuple'] = ROOT.TNtuple("pythia8-Geant4","min bias flux after 3m hadron absorber "+opt,tuples)
160 ROOT.gROOT.cd()
161 t.SetEventList(0)
162 if storeCharm: t.Draw(">>temp",cuts[opt]+"&"+OpCharm)
163 else: t.Draw(">>temp",cuts[opt]+"&"+noOpCharm)
164 temp = ROOT.gROOT.FindObjectAny('temp')
165 t.SetEventList(temp)
166 nev = temp.GetN()
167 leaves = t.GetListOfLeaves()
168 nL = leaves.GetEntries()
169 for iev in range(nev) :
170 rc = t.GetEntry(temp.GetEntry(iev))
171 vlist = []
172 k=-1
173 for x in range(nL):
174 k+=1
175 if nL > 11 and (k==7 or k==8 or k==9): continue
176 vlist.append( leaves.At(x).GetValue() )
177 if len(vlist) != 11 :
178 print("this should never happen, big error",len(vlist),k,p,iev,nev)
179 1/0
180 # "id:px:py:pz:x:y:z:pythiaid:parentid:w:ecut"
181 # yandex productions have
182 # "id:px:py:pz:x:y:z:ox:oy:oz:pythiaid:parentid:w:ecut"
183 Psq = vlist[1]**2+vlist[2]**2+vlist[3]**2
184 if abs(vlist[0])==13: Ekin = ROOT.TMath.Sqrt(Mmu2+Psq)-Mmu
185 else: Ekin = ROOT.TMath.Sqrt(Psq)
186 # re calculate weights
187 # pot = 5E13/w
188 # E > 100 GeV: add all pots -> w = 5E13/pot
189 # w = norm/( pot["CERN-Cracow"][100.] + productions["Yandex"][5.] + productions["Yandex2"][10.] )
190 # E < 100 & E > 10 GeV
191 # w = norm/( pot["CERN-Cracow"][10.] + productions["Yandex"][5.] + productions["Yandex2"][10.] )
192 # E < 10 & E > 5 GeV
193 # w = norm/( pot["CERN-Cracow"][1.] + productions["Yandex"][5.] )
194 # E < 5 GeV & E > 1 GeV
195 # w = norm/( pot["CERN-Cracow"][1.] + productions["Yandex"][0.5] )
196 # E < 1 GeV & E > 0.5 GeV
197 # w = norm/( productions["Yandex"][0.5] )
198 if Ekin > 100. :
199 vlist[9] = norm/( pot["CERN-Cracow"][100.] + pot["Yandex"][5.] + pot["Yandex2"][10.] )
200 elif Ekin > 10. :
201 vlist[9] = norm/( pot["CERN-Cracow"][10.] + pot["Yandex"][5.] + pot["Yandex2"][10.] )
202 elif Ekin > 5. :
203 vlist[9] = norm/( pot["CERN-Cracow"][1.] + pot["Yandex"][5.] )
204 elif Ekin > 1. :
205 vlist[9] = norm/( pot["CERN-Cracow"][1.] + pot["Yandex"][0.5] )
206 elif Ekin > 0.5 :
207 vlist[9] = norm/( pot["Yandex"][0.5] )
208 else :
209 print("this should not happen, except some rounding errors",p,Ekin,vlist[9])
210# scoring plane, g4Ex_gap: afterHadronZ = z0Pos+targetL+absorberL+5.1*cm
211# z0Pos = -50.*m absorberL = 2*150.*cm
212# target length increased for Yandex2 production, ignore this, put all muons at current end of hadronabsorber
213 vlist[6] = endOfHadronAbsorber
214 h['ntuple'].Fill(vlist[0],vlist[1],vlist[2],vlist[3],vlist[4],vlist[5],vlist[6],
215 vlist[7],vlist[8],vlist[9],vlist[10])
216 h['N'].cd()
217 h['ntuple'].Write()
218 h['N'].Close()
219
220def runProduction(opts=''):
221 we = fillWeights()
222 pot = {}
223 norm = 5.E13
224 for p in we:
225 pot[p]={}
226 for w in we[p]:
227 pot[p][we[p][w][0]] = 5.E13/w
228 print("pots:",pot)
229 #
230 mergeMinBias(pot,norm=5.E13,opt=opts)
231
232def removeCharm(p):
233 f = ROOT.TFile.Open(eospath+productions[p]["file"])
234 t = f.FindObjectAny("pythia8-Geant4")
235 first = True
236 if first:
237 first = False
238 tuples = ''
239 for l in t.GetListOfLeaves():
240 if tuples == '': tuples += l.GetName()
241 else: tuples += ":"+l.GetName()
242 h['N'] = ROOT.TFile(fnew, 'RECREATE')
243 print('new file created',fnew)
244 h['ntuple'] = ROOT.TNtuple("pythia8-Geant4",t.GetTitle()+" no charm",tuples)
245 ROOT.gROOT.cd()
246 t.SetEventList(0)
247 t.Draw(">>temp",noOpCharm)
248 temp = ROOT.gROOT.FindObjectAny('temp')
249 t.SetEventList(temp)
250 nev = temp.GetN()
251 leaves = t.GetListOfLeaves()
252 nL = leaves.GetEntries()
253 for iev in range(nev):
254 rc = t.GetEntry(temp.GetEntry(iev))
255 vlist = array('f')
256 for x in range(leaves.GetEntries()):
257 vlist.append( leaves.At(x).GetValue() )
258 h['ntuple'].Fill(vlist)
259 h['N'].cd()
260 h['ntuple'].Write()
261 h['N'].Close()
262
263def mergeWithCharm(splitOnly=False,ramOnly=False):
264 # Ntup.Fill(par.id(),par.px(),par.py(),par.pz(),par.e(),par.m(),wspill,sTree.id,sTree.px,sTree.py,sTree.pz,sTree.E,sTree.M)
265 # i.e. the par. is for the neutrino, and the sTree. is for its mother.
266 # wspill is the weight for this file normalised/5e13.
267 if not splitOnly:
268 fcascade = ROOT.TFile.Open(eospath+"/Charm/Decay-Cascade-parp16-MSTP82-1-MSEL4-76Mpot_1.root")
269 t = fcascade.Decay
270 newFile = ROOT.TFile("pythia8_Charm.root", 'RECREATE')
271 nt = ROOT.TNtuple("pythia8-Geant4","mu/nu flux from charm","id:px:py:pz:x:y:z:opx:opy:opz:ox:oy:oz:pythiaid:parentid:w:ecut")
272 for n in range(t.GetEntries()):
273 rc = t.GetEntry(n)
274 ztarget = rnr.Exp(0.16) + startOfTarget
275 vlist = array('f')
276 x = t.id,t.px,t.py,t.pz,0.,0.,ztarget,t.px,t.py,t.pz,0.,0.,ztarget,t.id,t.mid,t.weight,0.
277 for ax in x: vlist.append(ax)
278 rc = nt.Fill(vlist)
279 newFile.cd()
280 nt.Write()
281 newFile.Close()
282 fcascade.Close()
283 os.system("hadd -f pythia8_Geant4-withCharm.root pythia8_Geant4-noOpenCharm.root pythia8_Charm.root")
284 print(" progress: minbias and charm merged")
285 ramOnly = True
286 if ramOnly:
287# put all events in memory, otherwise will take years to finish
288 event = ROOT.std.vector("float")
289 f = ROOT.TFile("pythia8_Geant4-withCharm.root")
290 t = f.FindObjectAny('pythia8-Geant4')
291 leaves = t.GetListOfLeaves()
292 L = leaves.GetEntries()
293 m=0
294 allEvents = []
295 for n in range(t.GetEntries()):
296 rc = t.GetEvent(n)
297 if m%1000000==0 : print('status read',m)
298 m+=1
299 a = event()
300 for l in range(L): a.push_back(leaves.At(l).GetValue())
301 allEvents.append(a)
302# distribute events randomly
303 evList = []
304 for n in range(len(allEvents)): evList.append(n)
305 random.shuffle(evList)
306 leaves = t.GetListOfLeaves()
307 tuples = ''
308 for l in leaves:
309 if tuples=='': tuples += l.GetName()
310 else: tuples += ":"+l.GetName()
311 newFile = ROOT.TFile("pythia8_Geant4-withCharm-ram.root", 'RECREATE')
312 randomTuple = ROOT.TNtuple(t.GetName(),t.GetTitle(),tuples)
313 m=0
314 for n in evList:
315 a = allEvents[n]
316 if m%1000000==0 : print('status write',m)
317 m+=1
318 vlist = array('f')
319 for x in a: vlist.append(x)
320 randomTuple.Fill(vlist)
321 newFile.cd()
322 randomTuple.Write()
323 newFile.Close()
324 f.Close()
325 allEvents = []
326 print(" progress: order of events randomized")
327 if 1>0:
328# split in muons and neutrinos
329 cuts = {'_onlyNeutrinos':'abs(id)==12||abs(id)==14||abs(id)==16','_onlyMuons':'abs(id)==13'}
330 fName = "pythia8_Geant4-withCharm-ram.root"
331 f = ROOT.TFile(fName)
332 t = f.FindObjectAny("pythia8-Geant4")
333 tuples = ''
334# add histograms
335 for idnu in [16,-16,14,-14,12,-12]:
336 name = pdg.GetParticle(idnu).GetName()
337 idhnu=1000+idnu
338 if idnu < 0: idhnu=2000+abs(idnu)
339 ut.bookHist(h,str(idhnu),name+' momentum (GeV)',400,0.,400.)
340 ut.bookHist(h,str(idhnu+100),name+' log10(ptot) vs log10(pt+0.01)',100,-0.3,1.7,100,-2.,1.)
341 ut.bookHist(h,str(idhnu+200),name+' log10(ptot) vs log10(pt+0.01)',25,-0.3,1.7,100,-2.,1.)
342 leaves = t.GetListOfLeaves()
343 for l in leaves:
344 if tuples == '': tuples += l.GetName()
345 else: tuples += ":"+l.GetName()
346 for opt in cuts:
347 fxx = fName.replace('-ram.root',opt+'.root')
348 N = ROOT.TFile(fxx, 'RECREATE')
349 ntuple = ROOT.TNtuple("pythia8-Geant4",opt.replace("_only","")+" flux mbias and charm"+opt,tuples)
350 ROOT.gROOT.cd()
351 t.SetEventList(0)
352 t.Draw(">>temp",cuts[opt])
353 temp = ROOT.gROOT.FindObjectAny('temp')
354 t.SetEventList(temp)
355 for iev in range(temp.GetN()) :
356 rc = t.GetEntry(temp.GetEntry(iev))
357 vlist = array('f')
358 for n in range(leaves.GetSize()):
359 vlist.append(leaves.At(n).GetValue())
360 ntuple.Fill(vlist)
361 if "Neutrinos" in opt:
362 pt2=t.px**2+t.py**2
363 ptot=ROOT.TMath.Sqrt(pt2+t.pz**2)
364 l10ptot=min(max(ROOT.TMath.Log10(ptot),-0.3),1.69999)
365 l10pt=min(ROOT.TMath.Log10(ROOT.TMath.Sqrt(pt2)+0.01),0.9999)
366 idnu = int(t.id)
367 idhnu=1000+idnu
368 if idnu < 0: idhnu=2000+abs(idnu)
369 h[str(idhnu)].Fill(ptot,t.w)
370 h[str(idhnu+100)].Fill(l10ptot,l10pt,t.w)
371 h[str(idhnu+200)].Fill(l10ptot,l10pt,t.w)
372#
373 N.cd()
374 ntuple.Write()
375 if "Neutrinos" in opt:
376 for akey in h:
377 cln = h[akey].Class().GetName()
378 if not cln.find('TH')<0: h[akey].Write()
379 N.Close()
380 print(" progress: splitted "+opt)
381def test(fname):
382 h['f'] = ROOT.TFile.Open(fname)
383 sTree = h['f'].FindObjectAny('pythia8-Geant4')
384 TplotP(sTree)
385def compare():
386 test(eospath+'pythia8_Geant4_onlyMuons.root')
387 for x in ['','_>E']:
388 for z in ['p','pt']:
389 for ahist in ['mu-','mu+','mu','mu-charm','mu+charm','mucharm']:
390 h['TP'+z+ahist+x]=h['T'+z+ahist+x].Clone('CT'+ahist+x)
391 test(eospath+'pythia8_Geant4_Yandex_onlyNeutrinos.root')
392 for x in ['','_>E']:
393 for z in ['p','pt']:
394 for ahist in ['numusum','nuesum','numusumcharm','nuesumcharm','numu','numubar','nue','nuebar','numucharm','numubarcharm','nuecharm','nuebarcharm']:
395 h['TP'+z+ahist+x]=h['T'+z+ahist+x].Clone('CT'+ahist+x)
396 test(eospath+'Mbias/pythia8_Geant4-withCharm-ram.root')
397 for x in ['','_>E']:
398 for z in ['p','pt']:
399 t = z.upper()
400 if x != '' : t='>'+t
401 t1=h[t].cd(1)
402 p = z+'mu'
403 h['T'+p+x].SetTitle('musum')
404 h['T'+p+x].Draw()
405 h['T'+p+'charm'+x].Draw('same')
406 h['TP'+p+x].Draw('same')
407 h['TP'+p+x].SetLineColor(6)
408 h['TP'+p+'charm'+x].Draw('same')
409 h['TP'+p+'charm'+x].SetLineColor(3)
410 h['T'+p+'charm'+x].SetLineColor(2)
411 h['T'+p+x].SetLineColor(4)
412 h['L'+p+x] = ROOT.TLegend(0.33,0.69,0.99,0.94)
413 h['L'+p+x].AddEntry(h['T'+p+x],'muon new with charm, cascade, k-fac','PL')
414 h['L'+p+x].AddEntry(h['T'+p+'charm'+x],'muon from charm new with charm, cascade, k-fac','PL')
415 h['L'+p+x].AddEntry(h['TP'+p+x],'muon old CERN-Cracow prod','PL')
416 h['L'+p+x].AddEntry(h['TP'+p+'charm'+x],'muon from charm old CERN-Cracow prod','PL')
417 h['L'+p+x].Draw()
418 h[t].cd(2)
419 p = z+'numusum'
420 h['T'+p+x].Draw()
421 h['T'+p+'charm'+x].Draw('same')
422 h['TP'+p+x].Draw('same')
423 h['TP'+p+x].SetLineColor(6)
424 h['TP'+p+'charm'+x].Draw('same')
425 h['TP'+p+'charm'+x].SetLineColor(3)
426 h['T'+p+'charm'+x].SetLineColor(2)
427 h['T'+p+x].SetLineColor(4)
428 h['L'+p+x] = ROOT.TLegend(0.33,0.69,0.99,0.94)
429 h['L'+p+x].AddEntry(h['T'+p+x],'nu_mu new with charm, cascade, k-fac','PL')
430 h['L'+p+x].AddEntry(h['T'+p+'charm'+x],'nu_mu from charm new with charm, cascade, k-fac','PL')
431 h['L'+p+x].AddEntry(h['TP'+p+x],'nu_mu old CERN-Cracow prod','PL')
432 h['L'+p+x].AddEntry(h['TP'+p+'charm'+x],'nu_mu from charm old CERN-Cracow prod','PL')
433 h['L'+p+x].Draw()
434 t3=h[t].cd(3)
435 t3.SetLogy(1)
436 p = z+'nuesum'
437 h['T'+p+x].Draw()
438 h['T'+p+'charm'+x].Draw('same')
439 h['TP'+p+x].Draw('same')
440 h['TP'+p+x].SetLineColor(6)
441 h['TP'+p+'charm'+x].Draw('same')
442 h['TP'+p+'charm'+x].SetLineColor(3)
443 h['T'+p+'charm'+x].SetLineColor(2)
444 h['T'+p+x].SetLineColor(4)
445 h['L'+p+x] = ROOT.TLegend(0.33,0.69,0.99,0.94)
446 h['L'+p+x].AddEntry(h['T'+p+x],'nu_e new with charm, cascade, k-fac','PL')
447 h['L'+p+x].AddEntry(h['T'+p+'charm'+x],'nu_e from charm new with charm, cascade, k-fac','PL')
448 h['L'+p+x].AddEntry(h['TP'+p+x],'nu_e old CERN-Cracow prod','PL')
449 h['L'+p+x].AddEntry(h['TP'+p+'charm'+x],'nu_e from charm old CERN-Cracow prod','PL')
450 h['L'+p+x].Draw()
451#
452 t4 = h[t].cd(4)
453 t4.SetLogy(1)
454 h['Lmuc'+z+x] = ROOT.TLegend(0.33,0.73,0.99,0.94)
455 p = z+'mu-'
456 h['T'+p+x].SetTitle('mu-/mu+')
457 h['T'+p+x].Draw()
458 h['TP'+p+x].Draw('same')
459 h['TP'+p+x].SetLineColor(6)
460 h['T'+p+x].SetLineColor(4)
461 h['Lmuc'+z+x].AddEntry(h['T'+p+x],'mu- new with charm, cascade, k-fac','PL')
462 h['Lmuc'+z+x].AddEntry(h['TP'+p+x],'mu- old Yandex prod','PL')
463 p = z+'mu+'
464 h['T'+p+x].Draw('same')
465 h['TP'+p+x].Draw('same')
466 h['TP'+p+x].SetLineColor(3)
467 h['T'+p+x].SetLineColor(2)
468 h['Lmuc'+z+x].AddEntry(h['T'+p+x],'mu+ new with charm, cascade, k-fac','PL')
469 h['Lmuc'+z+x].AddEntry(h['TP'+p+x],'mu+ old Yandex prod','PL')
470 h['Lmuc'+z+x].Draw()
471#
472 t5 = h[t].cd(5)
473 t5.SetLogy(1)
474 h['Lnumu'+z+x] = ROOT.TLegend(0.33,0.73,0.99,0.94)
475 p = z+'numu'
476 h['T'+p+x].SetTitle('numu/numubar')
477 h['T'+p+x].Draw()
478 h['TP'+p+x].Draw('same')
479 h['TP'+p+x].SetLineColor(6)
480 h['T'+p+x].SetLineColor(4)
481 h['Lnumu'+z+x].AddEntry(h['T'+p+x],'nu_mu new with charm, cascade, k-fac','PL')
482 h['Lnumu'+z+x].AddEntry(h['TP'+p+x],'nu_mu old Yandex prod','PL')
483 p = z+'numubar'
484 h['T'+p+x].Draw('same')
485 h['TP'+p+x].Draw('same')
486 h['TP'+p+x].SetLineColor(3)
487 h['T'+p+x].SetLineColor(2)
488 h['Lnumu'+z+x].AddEntry(h['T'+p+x],'anti nu_mu new with charm, cascade, k-fac','PL')
489 h['Lnumu'+z+x].AddEntry(h['TP'+p+x],'anti nu_mu old Yandex prod','PL')
490 h['Lnumu'+z+x].Draw()
491 t6 = h[t].cd(6)
492 t6.SetLogy(1)
493 p = z+'nue'
494 h['Lnue'+z+x] = ROOT.TLegend(0.33,0.73,0.99,0.94)
495 h['T'+p+x].SetTitle('nue/nuebar')
496 h['T'+p+x].Draw()
497 h['TP'+p+x].Draw('same')
498 h['TP'+p+x].SetLineColor(6)
499 h['T'+p+x].SetLineColor(4)
500 h['Lnue'+z+x].AddEntry(h['T'+p+x],'nu_e new with charm, cascade, k-fac','PL')
501 h['Lnue'+z+x].AddEntry(h['TP'+p+x],'nu_e old Yandex prod','PL')
502 p = z+'nuebar'
503 h['T'+p+x].Draw('same')
504 h['TP'+p+x].Draw('same')
505 h['TP'+p+x].SetLineColor(3)
506 h['T'+p+x].SetLineColor(2)
507 h['Lnue'+z+x].AddEntry(h['T'+p+x],'anti nu_e new with charm, cascade, k-fac','PL')
508 h['Lnue'+z+x].AddEntry(h['TP'+p+x],'anti nu_e old Yandex prod','PL')
509 h['Lnue'+z+x].Draw()
510#
511 h[t].Print('comparison'+z+x.replace('_>','')+'.png')
512 h[t].Print('comparison'+z+x.replace('_>','')+'.pdf')
513# make ratio plots
514 x = '_>E'
515 for z in ['p','pt']:
516 h[z+'muRatio'+x]=h['T'+z+'mu'+x].Clone(z+'muRatio'+x)
517 h[z+'muRatio'+x].Divide(h['TP'+z+'mu'+x])
518 h[z+'numuRatio'+x]=h['T'+z+'numusum'+x].Clone(z+'numuRatio'+x)
519 h[z+'numuRatio'+x].Divide(h['TP'+z+'numusum'+x])
520 h[z+'nueRatio'+x]=h['T'+z+'nuesum'+x].Clone(z+'nueRatio'+x)
521 h[z+'nueRatio'+x].Divide(h['TP'+z+'nuesum'+x])
522 ut.bookCanvas(h,key='ratios',title='ratios',nx=1800,ny=600,cx=2,cy=1)
523 n = 1
524 for z in ['p','pt']:
525 h['Lratio'+z+x] = ROOT.TLegend(0.21,0.74,0.71,0.85)
526 tc = h['ratios'].cd(n)
527 n+=1
528 h[z+'muRatio'+x].SetLineColor(2)
529 h[z+'muRatio'+x].SetMaximum(max(h[z+'muRatio'+x].GetMaximum(),h[z+'numuRatio'+x].GetMaximum()))
530 h[z+'muRatio'+x].SetMinimum(0)
531 h[z+'muRatio'+x].SetStats(0)
532 h[z+'muRatio'+x].Draw()
533 h[z+'numuRatio'+x].SetLineColor(3)
534 h[z+'numuRatio'+x].SetStats(0)
535 h[z+'numuRatio'+x].Draw('same')
536 #h[z+'nueRatio'+x].SetLineColor(4)
537 #h[z+'nueRatio'+x].Draw('same')
538 h['Lratio'+z+x].AddEntry(h[z+'muRatio'+x],'muon flux new / old ','PL')
539 h['Lratio'+z+x].AddEntry(h[z+'numuRatio'+x],'nu_mu flux new / old ','PL')
540 #h['Lratio'+z+x].AddEntry(h[z+'nueRatio'+x],'nu_e flux new / old ','PL')
541 h['Lratio'+z+x].Draw()
542 h['ratios'].Print('comparisonRatios.png')
543 h['ratios'].Print('comparisonRatios.pdf')
544
545print("+ merging with charm events using existing charmless Mbias file: mergeWithCharm()")
546print("+ removeCharm(p) from mbias file made with g4Ex_gap_mergeFiles.py")
547print("+ testing output: test('pythia8_Geant4-noOpenCharm.root')")
548print("+ not used anymore: to start the full production, including merging of Mbias files: runProduction()")
549
mergeMbias.mergeMinBias
mergeMinBias(pot, norm=5.E13, opt='')
Definition mergeMbias.py:140
mergeMbias.mergeWithCharm
mergeWithCharm(splitOnly=False, ramOnly=False)
Definition mergeMbias.py:263
mergeMbias.runProduction
runProduction(opts='')
Definition mergeMbias.py:220
mergeMbias.test
test(fname)
Definition mergeMbias.py:381
mergeMbias.TplotP
TplotP(sTree)
Definition mergeMbias.py:42