FairShip
Loading...
Searching...
No Matches
Functions | Variables
mergeMbias Namespace Reference

Functions

def fillPart (t)
 
def fillWeights ()
 
None TplotP (sTree)
 
None mergeMinBias (pot, norm=5.0e13, opt="")
 
None runProduction (str opts="")
 
None removeCharm (p)
 
None mergeWithCharm (bool splitOnly=False, bool ramOnly=False)
 
None test (fname)
 
None compare ()
 

Variables

ROOT pdg = ROOT.TDatabasePDG()
 
ROOT mu = pdg.GetParticle(13)
 
ROOT Mmu = mu.Mass()
 
ROOT Mmu2 = Mmu * Mmu
 
ROOT rnr = ROOT.TRandom()
 
ROOT eospath = ROOT.gSystem.Getenv("EOSSHIP") + "/eos/experiment/ship/data/"
 
geometry_config ship_geo = geometry_config.create_config(Yheight=10.0)
 
tuple endOfHadronAbsorber = (ship_geo["hadronAbsorber"].z + ship_geo["hadronAbsorber"].length / 2.0) / 100.0
 
float startOfTarget = -50.0
 
dict productions = {}
 
bool allProds = False
 
str fnew = "pythia8_Geant4-noOpenCharm.root"
 
tuple noOpCharm
 
tuple OpCharm
 
dict cuts = {"": "abs(id)>0", "_onlyMuons": "abs(id)==13", "_onlyNeutrinos": "abs(id)==12||abs(id)==14||abs(id)==16"}
 
dict h = {}
 

Function Documentation

◆ compare()

None mergeMbias.compare ( )

Definition at line 474 of file mergeMbias.py.

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

◆ fillPart()

def mergeMbias.fillPart (   t)

Definition at line 23 of file mergeMbias.py.

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#

◆ fillWeights()

def mergeMbias.fillWeights ( )

Definition at line 34 of file mergeMbias.py.

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

◆ mergeMinBias()

None mergeMbias.mergeMinBias (   pot,
  norm = 5.0e13,
  opt = "" 
)

Definition at line 172 of file mergeMbias.py.

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

◆ mergeWithCharm()

None mergeMbias.mergeWithCharm ( bool   splitOnly = False,
bool   ramOnly = False 
)

Definition at line 315 of file mergeMbias.py.

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

◆ removeCharm()

None mergeMbias.removeCharm (   p)

Definition at line 281 of file mergeMbias.py.

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

◆ runProduction()

None mergeMbias.runProduction ( str   opts = "")

Definition at line 269 of file mergeMbias.py.

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

◆ test()

None mergeMbias.test (   fname)

Definition at line 468 of file mergeMbias.py.

468def test(fname) -> None:
469 h["f"] = ROOT.TFile.Open(fname)
470 sTree = h["f"].FindObjectAny("pythia8-Geant4")
471 TplotP(sTree)
472
473

◆ TplotP()

None mergeMbias.TplotP (   sTree)

Definition at line 50 of file mergeMbias.py.

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

Variable Documentation

◆ allProds

bool mergeMbias.allProds = False

Definition at line 144 of file mergeMbias.py.

◆ cuts

dict mergeMbias.cuts = {"": "abs(id)>0", "_onlyMuons": "abs(id)==13", "_onlyNeutrinos": "abs(id)==12||abs(id)==14||abs(id)==16"}

Definition at line 167 of file mergeMbias.py.

◆ endOfHadronAbsorber

tuple mergeMbias.endOfHadronAbsorber = (ship_geo["hadronAbsorber"].z + ship_geo["hadronAbsorber"].length / 2.0) / 100.0

Definition at line 19 of file mergeMbias.py.

◆ eospath

ROOT mergeMbias.eospath = ROOT.gSystem.Getenv("EOSSHIP") + "/eos/experiment/ship/data/"

Definition at line 17 of file mergeMbias.py.

◆ fnew

str mergeMbias.fnew = "pythia8_Geant4-noOpenCharm.root"

Definition at line 158 of file mergeMbias.py.

◆ h

dict mergeMbias.h = {}

Definition at line 169 of file mergeMbias.py.

◆ Mmu

ROOT mergeMbias.Mmu = mu.Mass()

Definition at line 14 of file mergeMbias.py.

◆ Mmu2

ROOT mergeMbias.Mmu2 = Mmu * Mmu

Definition at line 15 of file mergeMbias.py.

◆ mu

ROOT mergeMbias.mu = pdg.GetParticle(13)

Definition at line 13 of file mergeMbias.py.

◆ noOpCharm

tuple mergeMbias.noOpCharm
Initial value:
1= (
2 "!(pythiaid==id & (abs(parentid) == 15 || abs(parentid) == 4112 || abs(parentid) == 4122 || abs(parentid) == 4132 \
3 || abs(parentid) == 431 || abs(parentid) == 421 || abs(parentid) == 411) )"
4)

Definition at line 159 of file mergeMbias.py.

◆ OpCharm

tuple mergeMbias.OpCharm
Initial value:
1= (
2 "(pythiaid==id & (abs(parentid) == 15 || abs(parentid) == 4112 || abs(parentid) == 4122 || abs(parentid) == 4132 \
3 || abs(parentid) == 431 || abs(parentid) == 421 || abs(parentid) == 411) )"
4)

Definition at line 163 of file mergeMbias.py.

◆ pdg

ROOT mergeMbias.pdg = ROOT.TDatabasePDG()

Definition at line 12 of file mergeMbias.py.

◆ productions

dict mergeMbias.productions = {}

Definition at line 143 of file mergeMbias.py.

◆ rnr

ROOT mergeMbias.rnr = ROOT.TRandom()

Definition at line 16 of file mergeMbias.py.

◆ ship_geo

geometry_config mergeMbias.ship_geo = geometry_config.create_config(Yheight=10.0)

Definition at line 18 of file mergeMbias.py.

◆ startOfTarget

float mergeMbias.startOfTarget = -50.0

Definition at line 20 of file mergeMbias.py.