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