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

Functions

None origin (sTree, it)
 
None makeProd ()
 
def detMap ()
 
None fitSingleGauss (x, ba=None, be=None)
 
None bookHist (detName)
 
None BigEventLoop ()
 
None executeOneFile (fn, output=None, pid=None)
 
None makePlots (int nstations)
 
None AnaEventLoop ()
 
None muDISntuple (fn)
 
None analyzeConcrete ()
 
None rareEventEmulsion (str fname="rareEmulsion.txt")
 
None extractRareEvents (single=None)
 
None extractMuCloseByEvents (single=None)
 
None MergeRareEvents (list[str]|None runs=None)
 
None persistency ()
 
None reDraw (fn)
 
None printAndCopy (str|None prefix=None)
 
None drawBoth (str tag, hn)
 
None debugGeoTracks (sTree)
 
None eventsWithStrawPoints (i)
 
None eventsWithEntryPoints (i)
 
None depEnergy (sTree)
 
None originOfMuon (TextIOWrapper fout, int n, fn, nEvents)
 
None pers ()
 
def makeNicePrintout (list[str]|None x=None)
 
None readAndMergeHistos (prods)
 

Variables

bool local = False
 
bool parallel = True
 
mp output = mp.Queue()
 
list processes = []
 
list prefixes = []
 makeProd("muon812",10,False,True) # –< 831 copied back, done 16.3.2015 makeProd("muon822",10,True,True) makeProd("muon821",10,True,True) makeProd("muon822",10,True,True)
 
int withChain = 0
 
str pref = "muon"
 
os xx = os.path.abspath(".").lower()
 
str prefix = pref + prefix
 
 else :
 
str testdir = "."
 
str path = ""
 
ROOT fgeo = ROOT.TFile(testdir + "/" + f)
 
ROOT sGeo = fgeo.Get("FAIRGeom")
 
f inputFile = f.replace("geofile_full", "ship")
 
f tmp = inputFile.split(".")
 
 try :
 
float dy = float(tmp[1] + "." + tmp[2])
 
f inputFile1 = inputFile.replace("_D.", ".")
 
f inputFile2 = inputFile
 
ROOT PDG = ROOT.TDatabasePDG.Instance()
 
geometry_config ShipGeo = geometry_config.create_config(Yheight=dy)
 
ROOT run = ROOT.FairRunSim()
 
shipDet_conf modules = shipDet_conf.configure(run, ShipGeo)
 
float rz_inter = -1.0, 0.0
 
bool otherPhysList = False
 
bool noField = False
 
bool passive = False
 
ROOT top = sGeo.GetTopVolume()
 
ROOT muon = top.GetNode("MuonDetector_1")
 
ROOT mvol = muon.GetVolume()
 
ROOT zmuon = muon.GetMatrix().GetTranslation()[2]
 
tuple totl = (zmuon + mvol.GetShape().GetDZ()) / u.m
 
float ztarget = -100.0
 
list fchain = []
 
list fchainRec = []
 
f q1 = inputFile1 in os.listdir(path + prefix + str(i))
 
f q2 = inputFile2 in os.listdir(path + prefix + str(i))
 
f recFile1 = inputFile1.replace(".root", "_rec.root")
 
f recFile2 = inputFile2.replace(".root", "_rec.root")
 
f r1 = recFile1 in os.listdir(path + prefix + str(i))
 
f r2 = recFile2 in os.listdir(path + prefix + str(i))
 
str fname = path + prefix + str(i) + "/" + inputFile
 
f recFile = inputFile.replace(".root", "_rec.root")
 
dict h = {}
 
dict histlist = {}
 
dict txt = {}
 
dict labels = {}
 
def logVols = detMap()
 
dict histlistAll
 
dict hLiSc = {1: {}}
 
dict hMuon = {}
 
ROOT mom = ROOT.TVector3()
 
ROOT pos = ROOT.TVector3()
 

Function Documentation

◆ AnaEventLoop()

None ana_ShipMuon.AnaEventLoop ( )

Definition at line 933 of file ana_ShipMuon.py.

933def AnaEventLoop() -> None:
934 fout = open("rareEvents.txt", "w") # noqa: SIM115
935 for fn in fchainRec:
936 f = ROOT.TFile(fn)
937 if not f.FindObjectAny("cbmsim"):
938 print("skip file ", f.GetName())
939 continue
940 sTree = f.Get("cbmsim")
941 sTree.GetEvent(0)
942 if sTree.GetBranch("GeoTracks"):
943 sTree.SetBranchStatus("GeoTracks", 0)
944 nEvents = sTree.GetEntries()
945 for n in range(nEvents):
946 sTree.GetEntry(n)
947 if n == 0:
948 print("now at event ", n, f.GetName())
949 if len(sTree.MCTrack) > 1:
950 wg = sTree.MCTrack[1].GetWeight() # also works for neutrinos
951 else:
952 wg = sTree.MCTrack[0].GetWeight() # also works for neutrinos
953 i = -1
954 for atrack in sTree.FitTracks:
955 i += 1
956 fitStatus = atrack.getFitStatus()
957 if not fitStatus.isFitConverged():
958 continue
959 nmeas = atrack.getNumPoints()
960 chi2 = fitStatus.getChi2() / nmeas
961 fittedState = atrack.getFittedState()
962 P = fittedState.getMomMag()
963 fout.write("rare event with track %i, %s, %8.2F \n" % (n, f.GetName(), wg))
964 originOfMuon(fout, n, f.GetName(), nEvents)
965 fout.write("reco %i,%5.3F,%8.2F \n" % (nmeas, chi2, P))
966 mcPartKey = sTree.fitTrack2MC[i]
967 mcPart = sTree.MCTrack[mcPartKey]
968 for ahit in sTree.strawtubesPoint:
969 if ahit.GetTrackID() == mcPartKey:
970 fout.write(
971 "P when making hit %i, %8.2F \n"
972 % (
973 ahit.PdgCode(),
974 ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2 + ahit.GetPz() ** 2) / u.GeV,
975 )
976 )
977 break
978 if not mcPart:
979 continue
980 Ptruth = mcPart.GetP()
981 fout.write("Ptruth %i %8.2F \n" % (mcPart.GetPdgCode(), Ptruth / u.GeV))
982
983
984#
open
int open(const char *, int)
Opens a file descriptor.

◆ analyzeConcrete()

None ana_ShipMuon.analyzeConcrete ( )

Definition at line 1021 of file ana_ShipMuon.py.

1021def analyzeConcrete() -> None:
1022 h["fout"] = ROOT.TFile("muConcrete.root", "recreate")
1023 h["ntuple"] = ROOT.TNtuple("muons", "muon flux concrete", "id:px:py:pz:x:y:z:w")
1024 for m in ["", "mu", "V0"]:
1025 ut.bookHist(h, "conc_hitz" + m, "concrete hit z " + m, 100, -100.0, 100.0)
1026 ut.bookHist(h, "conc_hitzP" + m, "concrete hit z vs P" + m, 100, -100.0, 100.0, 100, 0.0, 25.0)
1027 ut.bookHist(h, "conc_hity" + m, "concrete hit y " + m, 100, -15.0, 15.0)
1028 ut.bookHist(h, "conc_p" + m, "concrete hit p " + m, 1000, 0.0, 400.0)
1029 ut.bookHist(h, "conc_pt" + m, "concrete hit pt " + m, 100, 0.0, 20.0)
1030 ut.bookHist(h, "conc_hitzy" + m, "concrete hit zy " + m, 100, -100.0, 100.0, 100, -15.0, 15.0)
1031 top = sGeo.GetTopVolume()
1032 magn = top.GetNode("magyoke_1")
1033 z0 = magn.GetMatrix().GetTranslation()[2] / u.m
1034 for fn in fchain:
1035 f = ROOT.TFile(fn)
1036 if not f.FindObjectAny("cbmsim"):
1037 print("skip file ", f.GetName())
1038 continue
1039 sTree = f.Get("cbmsim")
1040 nEvents = sTree.GetEntries()
1041 ROOT.gROOT.cd()
1042 for n in range(nEvents):
1043 sTree.GetEntry(n)
1044 if len(sTree.MCTrack) > 1:
1045 wg = sTree.MCTrack[1].GetWeight()
1046 else:
1047 wg = sTree.MCTrack[0].GetWeight()
1048 for ahit in sTree.vetoPoint:
1049 detID = ahit.GetDetectorID()
1050 if logVols[detID] != "rockD":
1051 continue
1052 m = ""
1053 pid = ahit.PdgCode()
1054 if abs(pid) == 13:
1055 m = "mu"
1056 P = ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2 + ahit.GetPz() ** 2)
1057 if abs(pid) == 13 and 3 / u.GeV < P:
1058 m = "V0"
1059 h["ntuple"].Fill(
1060 float(pid),
1061 float(ahit.GetPx() / u.GeV),
1062 float(ahit.GetPy() / u.GeV),
1063 float(ahit.GetPz() / u.GeV),
1064 float(ahit.GetX() / u.m),
1065 float(ahit.GetY() / u.m),
1066 float(ahit.GetZ() / u.m),
1067 float(wg),
1068 )
1069 h["conc_hitz" + m].Fill(ahit.GetZ() / u.m - z0, wg)
1070 h["conc_hity" + m].Fill(ahit.GetY() / u.m, wg)
1071 h["conc_p" + m].Fill(P / u.GeV, wg)
1072 h["conc_hitzP" + m].Fill(ahit.GetZ() / u.m, P / u.GeV, wg)
1073 Pt = ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2)
1074 h["conc_pt" + m].Fill(Pt / u.GeV, wg)
1075 h["conc_hitzy" + m].Fill(ahit.GetZ() / u.m - z0, ahit.GetY() / u.m, wg)
1076 #
1077 # start = [ahit.GetX()/u.m,ahit.GetY()/u.m,ahit.GetZ()/u.m]
1078 # direc = [-ahit.GetPx()/P,-ahit.GetPy()/P,-ahit.GetPz()/P]
1079 # t = - start[0]/direc[0]
1080
1081 ut.bookCanvas(h, key="ResultsV0", title="muons hitting concrete, p>3GeV", nx=1000, ny=600, cx=2, cy=2)
1082 ut.bookCanvas(h, key="Resultsmu", title="muons hitting concrete", nx=1000, ny=600, cx=2, cy=2)
1083 ut.bookCanvas(h, key="Results", title="hitting concrete", nx=1000, ny=600, cx=2, cy=2)
1084 for m in ["", "mu", "V0"]:
1085 tc = h["Results" + m].cd(1)
1086 h["conc_hity" + m].Draw()
1087 tc = h["Results" + m].cd(2)
1088 h["conc_hitz" + m].Draw()
1089 tc = h["Results" + m].cd(3)
1090 tc.SetLogy(1)
1091 h["conc_pt" + m].Draw()
1092 tc = h["Results" + m].cd(4)
1093 tc.SetLogy(1)
1094 h["conc_p" + m].Draw()
1095 h["fout"].cd()
1096 h["ntuple"].Write()
1097
1098

◆ BigEventLoop()

None ana_ShipMuon.BigEventLoop ( )

Definition at line 585 of file ana_ShipMuon.py.

585def BigEventLoop() -> None:
586 pid = 1
587 for fn in fchain:
588 if os.path.islink(fn):
589 rfn = os.path.realpath(fn).split("eos")[1]
590 fn = ROOT.gSystem.Getenv("EOSSHIP") + "/eos/" + rfn
591 elif not os.path.isfile(fn):
592 raise RuntimeError("Don't know what to do with " + fn)
593 if parallel:
594 # process files parallel instead of sequential
595 processes.append(mp.Process(target=executeOneFile, args=(fn, output, pid)))
596 pid += 1
597 else:
598 processes.append(fn)
599 # Run processes
600 n = 0
601 for p in processes:
602 if parallel:
603 p.start()
604 n += 1
605 else:
606 executeOneFile(p)
607 if parallel:
608 # Exit the completed processes
609 for p in processes:
610 p.join()
611 # clean histos before reading in the new ones
612 for x in h:
613 h[x].Reset()
614 print("now, collect the output")
615 pid = 1
616 for p in processes:
617 ut.readHists(h, "tmpHists_" + str(pid) + ".root")
618 pid += 1
619 # compactify liquid scintillator
620 for mu in ["", "_mu", "_muV0"]:
621 for x in ["", "_E", "_P", "_LP", "_OP", "_id", "_mul", "_evmul", "_origin", "_originmu"]:
622 for k in [1, 2, 3, 5]:
623 first = True
624 for j in hLiSc[k]:
625 detName = hLiSc[k][j]
626 tag = detName + mu + x
627 newh = detName[0:2] + "LiSc" + mu + x
628 if tag not in h:
629 continue
630 if first:
631 h[newh] = h[tag].Clone(newh)
632 h[newh].SetTitle(h[tag].GetTitle().split("_")[0])
633 first = False
634 else:
635 h[newh].Add(h[tag])
636 # compactify muon stations
637 for mu in ["", "_mu", "_muV0"]:
638 for x in ["", "_E", "_P", "_LP", "_OP", "_id", "_mul", "_evmul", "_origin", "_originmu"]:
639 first = True
640 for j in hMuon:
641 detName = hMuon[j]
642 tag = detName + mu + x
643 newh = "muondet" + mu + x
644 if first:
645 h[newh] = h[tag].Clone(newh)
646 h[newh].SetTitle(h[tag].GetTitle().split(" ")[0] + " " + newh)
647 first = False
648 else:
649 h[newh].Add(h[tag])
650
651 # make list of hists with entries
652 k = 1
653 for x in histlistAll:
654 if histlistAll[x] in h:
655 histlist[k] = histlistAll[x]
656 # make cumulative histograms
657 for c in ["", "_E", "_P", "_LP", "_OP", "_id", "_mul", "_evmul", "_origin", "_originmu"]:
658 h[histlist[k] + "_mu" + c].Add(h[histlist[k] + "_muV0" + c])
659 h[histlist[k] + c].Add(h[histlist[k] + "_mu" + c])
660 h[histlist[k] + c].SetMinimum(0.0)
661 h[histlist[k] + "_mu" + c].SetMinimum(0.0)
662 h[histlist[k] + "_muV0" + c].SetMinimum(0.0)
663 k += 1
664 nstations = len(histlist)
665 makePlots(nstations)
666
667

◆ bookHist()

None ana_ShipMuon.bookHist (   detName)

Definition at line 513 of file ana_ShipMuon.py.

513def bookHist(detName) -> None:
514 for mu in ["", "_mu", "_muV0"]:
515 tag = detName + mu
516 if detName.find("LS") < 0:
517 ut.bookHist(h, tag, "x/y " + detName, 100, -3.0, 3.0, 100, -6.0, 6.0)
518 else:
519 ut.bookHist(h, tag, "z phi " + detName, 100, -50.0, 50.0, 100, -1.0, 1.0)
520 ut.bookHist(h, tag + "_E", "deposited energy MeV " + detName, 100, 0.0, 2.0)
521 ut.bookHist(h, tag + "_P", "particle mom GeV " + detName, 500, 0.0, 50.0)
522 ut.bookHist(h, tag + "_LP", "particle mom GeV " + detName, 100, 0.0, 1.0)
523 ut.bookHist(h, tag + "_OP", "original particle mom GeV " + detName, 400, 0.0, 400.0)
524 ut.bookHist(h, tag + "_id", "particle id " + detName, 5001, -2499.5, 2499.5)
525 ut.bookHist(h, tag + "_mul", "multiplicity of hits/tracks " + detName, 100, -0.5, 99.5)
526 ut.bookHist(h, tag + "_evmul", "multiplicity of hits/event " + detName, 100, -0.5, 9999.5)
527 ut.bookHist(h, tag + "_origin", "r vs z", 100, ztarget, totl, 100, 0.0, 12.0)
528 ut.bookHist(h, tag + "_originmu", "r vs z", 100, ztarget, totl, 100, 0.0, 12.0)
529
530

◆ debugGeoTracks()

None ana_ShipMuon.debugGeoTracks (   sTree)

Definition at line 1300 of file ana_ShipMuon.py.

1300def debugGeoTracks(sTree) -> None:
1301 for i in range(sTree.GetEntries()):
1302 sTree.GetEntry(i)
1303 n = 0
1304 for gt in sTree.GeoTracks:
1305 print(
1306 i, n, gt.GetFirstPoint()[2], gt.GetLastPoint()[2], gt.GetParticle().GetPdgCode(), gt.GetParticle().P()
1307 )
1308 n += 1
1309
1310

◆ depEnergy()

None ana_ShipMuon.depEnergy (   sTree)

Definition at line 1344 of file ana_ShipMuon.py.

1344def depEnergy(sTree) -> None:
1345 for n in range(sTree.GetEntries()):
1346 sTree.GetEntry(n)
1347 for ahit in sTree.strawtubesPoint:
1348 dE = ahit.GetEnergyLoss() / u.keV
1349 ahit.Momentum(mom)
1350 pa = PDG.GetParticle(ahit.PdgCode())
1351 mpa = pa.Mass()
1352 E = ROOT.TMath.Sqrt(mom.Mag2() + mpa**2)
1353 ekin = E - mpa
1354 h["dE"].Fill(dE, ekin / u.MeV)
1355 h["dE"].SetXTitle("keV")
1356 h["dE"].SetYTitle("MeV")
1357
1358

◆ detMap()

def ana_ShipMuon.detMap ( )

Definition at line 473 of file ana_ShipMuon.py.

473def detMap():
474 sGeo = ROOT.gGeoManager
475 detList = {}
476 for v in sGeo.GetListOfVolumes():
477 nm = v.GetName()
478 i = sGeo.FindVolumeFast(nm).GetNumber()
479 detList[i] = nm
480 return detList
481
482

◆ drawBoth()

None ana_ShipMuon.drawBoth ( str  tag,
  hn 
)

Definition at line 1289 of file ana_ShipMuon.py.

1289def drawBoth(tag: str, hn) -> None:
1290 n1 = h[hn + "_mu" + tag].GetMaximum()
1291 n2 = h[hn + tag].GetMaximum()
1292 if n1 > n2:
1293 h[hn + tag].SetMaximum(n1)
1294 h[hn + "_mu" + tag].SetLineColor(4)
1295 h[hn + tag].SetLineColor(3)
1296 h[hn + "_mu" + tag].Draw()
1297 h[hn + tag].Draw("same")
1298
1299

◆ eventsWithEntryPoints()

None ana_ShipMuon.eventsWithEntryPoints (   i)

Definition at line 1328 of file ana_ShipMuon.py.

1328def eventsWithEntryPoints(i) -> None:
1329 sTree = fchain[i].Get("cbmsim")
1330 mom = ROOT.TVector3()
1331 for i in range(sTree.GetEntries()):
1332 sTree.GetEntry(i)
1333 np = 0
1334 for vp in sTree.vetoPoint:
1335 detName = logVols[vp.GetDetectorID()]
1336 if detName == "VetoTimeDet":
1337 np += 0 # ??
1338 vp.Momentum(mom)
1339 print(i, detName, vp.PdgCode())
1340 mom.Print()
1341 print("-----------------------")
1342
1343

◆ eventsWithStrawPoints()

None ana_ShipMuon.eventsWithStrawPoints (   i)

Definition at line 1311 of file ana_ShipMuon.py.

1311def eventsWithStrawPoints(i) -> None:
1312 sTree = fchain[i].Get("cbmsim")
1313 mom = ROOT.TVector3()
1314 for i in range(sTree.GetEntries()):
1315 sTree.GetEntry(i)
1316 nS = len(sTree.strawtubesPoint)
1317 if nS > 0:
1318 mu = sTree.MCTrack[0]
1319 mu.GetMomentum(mom)
1320 print(i, nS)
1321 mom.Print()
1322 sp = sTree.strawtubesPoint[(max(0, nS - 3))]
1323 sp.Momentum(mom)
1324 mom.Print()
1325 print("-----------------------")
1326
1327

◆ executeOneFile()

None ana_ShipMuon.executeOneFile (   fn,
  output = None,
  pid = None 
)

Definition at line 668 of file ana_ShipMuon.py.

668def executeOneFile(fn, output=None, pid=None) -> None:
669 f = ROOT.TFile.Open(fn)
670 sTree = f.Get("cbmsim")
671 nEvents = sTree.GetEntries()
672 if sTree.GetBranch("GeoTracks"):
673 sTree.SetBranchStatus("GeoTracks", 0)
674 sTree.GetEntry(0)
675 hitContainers = [
676 sTree.vetoPoint,
677 sTree.muonPoint,
678 sTree.EcalPointLite,
679 sTree.strawtubesPoint,
680 sTree.ShipRpcPoint,
681 sTree.TargetPoint,
682 ]
683 ntot = 0
684 for n in range(nEvents):
685 sTree.GetEntry(n)
686 theMuon = sTree.MCTrack[0]
687 if len(sTree.MCTrack) > 1:
688 w = sTree.MCTrack[1].GetWeight() # also works for neutrinos
689 else:
690 print("should not happen with new files", n, fn)
691 w = sTree.MCTrack[0].GetWeight() # also works for neutrinos
692 if w == 0:
693 w = 1.0
694 h["weight"].Fill(w)
695 h["muonP"].Fill(theMuon.GetP() / u.GeV, w)
696 ntot += 1
697 if ntot % 10000 == 0:
698 print("read event", f.GetName(), n)
699 hitmult = {}
700 esum = 0
701 for mcp in sTree.MCTrack:
702 if mcp.GetMotherId() == 0: # mother is original muon
703 E = mcp.GetEnergy()
704 if mcp.GetPdgCode() == 11:
705 h["deltaElec"].Fill(E, w)
706 esum += E
707 h["secondaries"].Fill(E, w)
708 h["prop_deltaElec"].Fill(
709 esum / sTree.MCTrack[0].GetP(), w
710 ) # until energy is really made persistent GetEnergy()
711 for c in hitContainers:
712 for ahit in c:
713 if not ahit.GetEnergyLoss() > 0:
714 continue
715 detID = ahit.GetDetectorID()
716 if ahit.GetName() == "strawtubesPoint":
717 detName = "strawstation_" + str(ahit.GetStationNumber())
718 x = ahit.GetX()
719 y = ahit.GetY()
720 E = ahit.GetEnergyLoss()
721 else:
722 if detID not in logVols:
723 detName = c.GetName().replace("Points", "")
724 if detName not in histlistAll.values():
725 print(detID, detName, c.GetName())
726 else:
727 detName = logVols[detID]
728 x = ahit.GetX()
729 y = ahit.GetY()
730 ahit.GetZ()
731 E = ahit.GetEnergyLoss()
732 if detName not in h:
733 bookHist(detName)
734 mu = ""
735 pdgID = -2
736 if "PdgCode" in dir(ahit):
737 pdgID = ahit.PdgCode()
738 trackID = ahit.GetTrackID()
739 phit = -100.0
740 mom = ROOT.TVector3()
741 if not trackID < 0:
742 aTrack = sTree.MCTrack[trackID]
743 pdgID = aTrack.GetPdgCode()
744 aTrack.GetMomentum(mom) # ! this is not momentum of particle at Calorimeter place
745 phit = mom.Mag() / u.GeV
746 if abs(pdgID) == 13:
747 mu = "_mu"
748 if phit > 3 and abs(pdgID) == 13:
749 mu = "_muV0"
750 detName = detName + mu
751 if detName.find("LS") < 0:
752 h[detName].Fill(x / u.m, y / u.m, w)
753 else:
754 h[detName].Fill(ahit.GetZ() / u.m, ROOT.TMath.ATan2(y, x) / ROOT.TMath.Pi(), w)
755 h[detName + "_E"].Fill(E / u.MeV, w)
756 if detName not in hitmult:
757 hitmult[detName] = {-1: 0}
758 if trackID not in hitmult[detName]:
759 hitmult[detName][trackID] = 0
760 hitmult[detName][trackID] += 1
761 hitmult[detName][-1] += 1
762 h[detName + "_id"].Fill(pdgID, w)
763 h[detName + "_P"].Fill(phit, w)
764 h[detName + "_LP"].Fill(phit, w)
765 if not trackID < 0:
766 r = ROOT.TMath.Sqrt(aTrack.GetStartX() ** 2 + aTrack.GetStartY() ** 2) / u.m
767 h["origin"].Fill(aTrack.GetStartZ() / u.m, r, w)
768 h[detName + "_origin"].Fill(aTrack.GetStartZ() / u.m, r, w)
769 if abs(pdgID) == 13:
770 h[detName + "_originmu"].Fill(aTrack.GetStartZ() / u.m, r, w)
771 h["borigin"].Fill(aTrack.GetStartZ() / u.m, r, w)
772 aTrack.GetMomentum(mom)
773 h[detName + "_OP"].Fill(mom.Mag() / u.GeV, w)
774 if trackID > 0:
775 origin(sTree, trackID)
776 h["porigin"].Fill(
777 aTrack.GetStartZ() / u.m,
778 ROOT.TMath.Sqrt(aTrack.GetStartX() ** 2 + aTrack.GetStartY() ** 2) / u.m,
779 w,
780 )
781 h["mu-inter"].Fill(rz_inter[1] / u.m, rz_inter[0] / u.m, w)
782 for detName in hitmult:
783 h[detName + "_evmul"].Fill(hitmult[detName][-1], w)
784 for tr in hitmult[detName]:
785 h[detName + "_mul"].Fill(hitmult[detName][tr], w)
786 if output:
787 ut.writeHists(h, "tmpHists_" + str(pid) + ".root")
788 output.put("ok")
789 f.Close()
790
791
792#

◆ extractMuCloseByEvents()

None ana_ShipMuon.extractMuCloseByEvents (   single = None)

Definition at line 1192 of file ana_ShipMuon.py.

1192def extractMuCloseByEvents(single=None) -> None:
1193 mom = ROOT.TVector3()
1194 pos = ROOT.TVector3()
1195 # Goliath magnet has been removed from geometry
1196 goliath_node = top.GetNode("volGoliath_1")
1197 if goliath_node:
1198 golmx = goliath_node.GetMatrix()
1199 zGol = golmx.GetTranslation()[2]
1200 else:
1201 zGol = -9999.0 # No Goliath, accept all z positions
1202 for fn in fchainRec:
1203 if single:
1204 if fn.find(str(single)) < 0:
1205 continue
1206 f = ROOT.TFile(fn)
1207 if not f.FindObjectAny("cbmsim"):
1208 print("skip file ", f.GetName())
1209 continue
1210 sTree = f.Get("cbmsim")
1211 nEvents = sTree.GetEntries()
1212 raref = ROOT.TFile(fn.replace(".root", "_clby.root"), "recreate")
1213 newTree = sTree.CloneTree(0)
1214 for n in range(nEvents):
1215 sTree.GetEntry(n)
1216 if n == 0:
1217 print("now at event ", n, f.GetName())
1218 # look for muons p>3 hitting something
1219 n = 0
1220 for c in [sTree.vetoPoint, sTree.strawtubesPoint, sTree.ShipRpcPoint]:
1221 for ahit in c:
1222 if abs(ahit.PdgCode()) != 13:
1223 continue
1224 ahit.Momentum(mom)
1225 if mom.Mag() < 3.0:
1226 continue
1227 ahit.Position(pos)
1228 if pos.z() < zGol:
1229 continue
1230 n += 1
1231 if n > 0:
1232 newTree.Fill()
1233 # print 'filled newTree',rc
1234 sTree.Clear()
1235 newTree.AutoSave()
1236 print(" --> events saved:", newTree.GetEntries())
1237 f.Close()
1238 raref.Close()
1239
1240
1241#

◆ extractRareEvents()

None ana_ShipMuon.extractRareEvents (   single = None)

Definition at line 1155 of file ana_ShipMuon.py.

1155def extractRareEvents(single=None) -> None:
1156 for fn in fchainRec:
1157 if single:
1158 if fn.find(str(single)) < 0:
1159 continue
1160 f = ROOT.TFile(fn)
1161 if not f.FindObjectAny("cbmsim"):
1162 print("skip file ", f.GetName())
1163 continue
1164 sTree = f.Get("cbmsim")
1165 nEvents = sTree.GetEntries()
1166 raref = ROOT.TFile(fn.replace(".root", "_rare.root"), "recreate")
1167 newTree = sTree.CloneTree(0)
1168 for n in range(nEvents):
1169 sTree.GetEntry(n)
1170 if n == 0:
1171 print("now at event ", n, f.GetName())
1172 # count fitted tracks which have converged and nDF>20:
1173 n = 0
1174 for fT in sTree.FitTracks:
1175 fst = fT.getFitStatus()
1176 if not fst.isFitConverged():
1177 continue
1178 if fst.getNdf() < 20:
1179 continue
1180 n += 1
1181 if n > 0:
1182 rc = newTree.Fill()
1183 print("filled newTree", rc)
1184 sTree.Clear()
1185 newTree.AutoSave()
1186 print(" --> events saved:", newTree.GetEntries())
1187 f.Close()
1188 raref.Close()
1189
1190
1191#

◆ fitSingleGauss()

None ana_ShipMuon.fitSingleGauss (   x,
  ba = None,
  be = None 
)

Definition at line 483 of file ana_ShipMuon.py.

483def fitSingleGauss(x, ba=None, be=None) -> None:
484 name = "myGauss_" + x
485 myGauss = h[x].GetListOfFunctions().FindObject(name)
486 if not myGauss:
487 if not ba:
488 ba = h[x].GetBinCenter(1)
489 if not be:
490 be = h[x].GetBinCenter(h[x].GetNbinsX())
491 bw = h[x].GetBinWidth(1)
492 mean = h[x].GetMean()
493 sigma = h[x].GetRMS()
494 norm = h[x].GetEntries() * 0.3
495 myGauss = ROOT.TF1(name, "[0]*" + str(bw) + "/([2]*sqrt(2*pi))*exp(-0.5*((x-[1])/[2])**2)+[3]", 4)
496 myGauss.SetParameter(0, norm)
497 myGauss.SetParameter(1, mean)
498 myGauss.SetParameter(2, sigma)
499 myGauss.SetParameter(3, 1.0)
500 myGauss.SetParName(0, "Signal")
501 myGauss.SetParName(1, "Mean")
502 myGauss.SetParName(2, "Sigma")
503 h[x].Fit(myGauss, "", "", ba, be)
504
505

◆ makeNicePrintout()

def ana_ShipMuon.makeNicePrintout ( list[str] | None   x = None)

Definition at line 1406 of file ana_ShipMuon.py.

1406def makeNicePrintout(x: list[str] | None = None):
1407 if x is None:
1408 x = ["rareEvents_61-62.txt", "rareEvents_71-72.txt"]
1409 result = []
1410 cor = 1.0
1411 for fn in x:
1412 with open(fn) as f:
1413 recTrack = None
1414 if fn == "rareEvents_81-102.txt":
1415 cor = 30.0
1416 for lx in f.readlines():
1417 line = lx.replace("\n", "")
1418 if not line.find("rare event") < 0:
1419 if recTrack:
1420 result.append(recTrack)
1421 tmp = line.split(",")
1422 w = tmp[2].replace(" ", "")
1423 ff = tmp[1].split("/")[0].replace(" ", "")
1424 recTrack = {"w": w, "file": ff}
1425 elif not line.find("original") < 0:
1426 tmp = line.split(",")
1427 recTrack["origin"] = tmp[0].split(" ")[2]
1428 recTrack["pytiaid"] = tmp[1].replace(" ", "")
1429 recTrack["o-mom"] = tmp[2].replace(" ", "")
1430 elif not line.find("reco ") < 0:
1431 tmp = line.split(",")
1432 recTrack["nmeas"] = tmp[0].split(" ")[1]
1433 recTrack["chi2"] = tmp[1]
1434 recTrack["p_rec"] = tmp[2].replace(" ", "")
1435 elif not line.find("making") < 0:
1436 tmp = line.split(",")
1437 recTrack["p_hit"] = tmp[1].replace(" ", "")
1438 recTrack["fp_hit"] = float(tmp[1].replace(" ", ""))
1439 elif not line.find("Ptruth") < 0:
1440 tmp = line.split(" ")
1441 recTrack["id_hit"] = tmp[1].replace(" ", "")
1442 # print a table
1443 print(
1444 "%4s %8s %8s %4s %8s %8s %8s %8s %8s %8s "
1445 % ("nr", "origin", "pythiaID", "ID", "p_orig", "p_hit", "chi2", "weight", "file", "cor w")
1446 )
1447 # sort according to p_hit
1448 tmp = sorted(result, key=itemgetter("fp_hit"))
1449 muonrate1 = 0
1450 muonrate2 = 0
1451 muonrate3 = 0
1452 for i in range(len(tmp)):
1453 tr = tmp[i]
1454 corw = float(tr["w"]) / cor
1455 if float(tr["p_hit"]) > 1:
1456 muonrate1 += corw
1457 if float(tr["p_hit"]) > 2:
1458 muonrate2 += corw
1459 if float(tr["p_hit"]) > 3:
1460 muonrate3 += corw
1461 print(
1462 "%4i %8s %8s %4s %8s %8s %8s %8s %8s %8s "
1463 % (
1464 i,
1465 tr["origin"],
1466 tr["pytiaid"],
1467 tr["id_hit"],
1468 tr["o-mom"],
1469 tr["p_hit"],
1470 tr["chi2"],
1471 tr["w"],
1472 tr["file"],
1473 corw,
1474 )
1475 )
1476 print("guestimate for muonrate above 1GeV = ", muonrate1)
1477 print("guestimate for muonrate above 2GeV = ", muonrate2)
1478 print("guestimate for muonrate above 3GeV = ", muonrate3)
1479 # guestimate for muonrate above 1GeV = 56025.4793333
1480 # guestimate for muonrate above 2GeV = 25584.9546667
1481 # guestimate for muonrate above 3GeV = 14792.8113333
1482 return tmp
1483
1484
1485#

◆ makePlots()

None ana_ShipMuon.makePlots ( int  nstations)

Definition at line 793 of file ana_ShipMuon.py.

793def makePlots(nstations: int) -> None:
794 cxcy = {
795 1: [2, 1],
796 2: [3, 1],
797 3: [2, 2],
798 4: [3, 2],
799 5: [3, 2],
800 6: [3, 2],
801 7: [3, 3],
802 8: [3, 3],
803 9: [3, 3],
804 10: [4, 3],
805 11: [4, 3],
806 12: [4, 3],
807 13: [5, 3],
808 14: [5, 3],
809 15: [5, 3],
810 }
811 ut.bookCanvas(
812 h,
813 key="ResultsI",
814 title="hit occupancies",
815 nx=1100,
816 ny=600 * cxcy[nstations][1],
817 cx=cxcy[nstations][0],
818 cy=cxcy[nstations][1],
819 )
820 ut.bookCanvas(
821 h,
822 key="ResultsImu",
823 title="hit occupancies",
824 nx=1100,
825 ny=600 * cxcy[nstations][1],
826 cx=cxcy[nstations][0],
827 cy=cxcy[nstations][1],
828 )
829 ut.bookCanvas(
830 h,
831 key="ResultsImuV0",
832 title="hit occupancies",
833 nx=1100,
834 ny=600 * cxcy[nstations][1],
835 cx=cxcy[nstations][0],
836 cy=cxcy[nstations][1],
837 )
838 ut.bookCanvas(
839 h,
840 key="ResultsII",
841 title="deposited energies",
842 nx=1100,
843 ny=600 * cxcy[nstations][1],
844 cx=cxcy[nstations][0],
845 cy=cxcy[nstations][1],
846 )
847 ut.bookCanvas(
848 h,
849 key="ResultsIII",
850 title="track info",
851 nx=1100,
852 ny=600 * cxcy[nstations][1],
853 cx=cxcy[nstations][0],
854 cy=cxcy[nstations][1],
855 )
856 ut.bookCanvas(
857 h,
858 key="ResultsIV",
859 title="track info",
860 nx=1100,
861 ny=600 * cxcy[nstations][1],
862 cx=cxcy[nstations][0],
863 cy=cxcy[nstations][1],
864 )
865 ut.bookCanvas(h, key="ResultsV", title="origin info", nx=600, ny=400, cx=1, cy=1)
866 txt[0] = "occupancy 1sec 5E13pot "
867 h["dummy"].Fill(-1)
868 nSpills = len(prefixes)
869 for i in histlist:
870 tc = h["ResultsI"].cd(i)
871 hn = histlist[i]
872 if hn not in h:
873 continue
874 h[hn].SetStats(0)
875 h[hn].Draw("colz")
876 txt[i] = "%5.2G" % (h[hn].GetSumOfWeights() / nSpills)
877 labels[i] = ROOT.TLatex().DrawLatexNDC(0.15, 0.85, txt[i])
878 #
879 hn = histlist[i] + "_mu"
880 tc = h["ResultsImu"].cd(i)
881 h[hn].SetStats(0)
882 h[hn].Draw("colz")
883 txt[i + 100] = "%5.2G" % (h[hn].GetSumOfWeights() / nSpills)
884 labels[i + 100] = ROOT.TLatex().DrawLatexNDC(0.15, 0.85, txt[i + 100])
885 #
886 hn = histlist[i] + "_muV0"
887 tc = h["ResultsImuV0"].cd(i)
888 h[hn].SetStats(0)
889 h[hn].Draw("colz")
890 txt[i + 200] = "%5.2G" % (h[hn].GetSumOfWeights() / nSpills)
891 labels[i + 200] = ROOT.TLatex().DrawLatexNDC(0.15, 0.85, txt[i + 200])
892 #
893 for i in histlist:
894 tc = h["ResultsII"].cd(i)
895 h[histlist[i] + "_E"].Draw("colz")
896 #
897 for i in histlist:
898 tc = h["ResultsIII"].cd(i)
899 tc.SetLogy(1)
900 hn = histlist[i]
901 drawBoth("_P", hn)
902 hid = h[hn + "_id"]
903 print("particle statistics for " + histlist[i])
904 for k in range(hid.GetNbinsX()):
905 ncont = hid.GetBinContent(k + 1)
906 pid = hid.GetBinCenter(k + 1)
907 if ncont > 0:
908 temp = int(abs(pid) + 0.5) * int(pid / abs(pid))
909 nm = PDG.GetParticle(temp).GetName()
910 print(f"{nm} :{ncont:5.2g}")
911 hid = h[histlist[i] + "_mu_id"]
912 for k in range(hid.GetNbinsX()):
913 ncont = hid.GetBinContent(k + 1)
914 pid = hid.GetBinCenter(k + 1)
915 if ncont > 0:
916 temp = int(abs(pid) + 0.5) * int(pid / abs(pid))
917 nm = PDG.GetParticle(temp).GetName()
918 print(f"{nm} :{ncont:5.2g}")
919 #
920 tc = h["ResultsV"].cd(1)
921 h["origin"].SetStats(0)
922 h["origin"].Draw("colz")
923 #
924 for i in histlist:
925 tc = h["ResultsIV"].cd(i)
926 tc.SetLogy(1)
927 hn = histlist[i]
928 drawBoth("_OP", hn)
929 persistency()
930
931
932#

◆ makeProd()

None ana_ShipMuon.makeProd ( )

Definition at line 449 of file ana_ShipMuon.py.

449def makeProd() -> None:
450 ntot = 736406
451 ncpu = 4
452 n3 = int(ntot / ncpu)
453 cmd = "python $FAIRSHIP/macro/run_simScript.py --MuonBack -f $SHIPSOFT/data/pythia8_Geant4_onlyMuons.root " # --display"
454 ns = 0
455 prefix = "muon18"
456 for i in range(1, ncpu + 1):
457 d = prefix + str(i)
458 if d not in os.listdir("."):
459 os.system("mkdir " + d)
460 os.chdir("./" + prefix + "1")
461 for i in range(1, ncpu + 1):
462 if i == ncpu:
463 n3 = ntot - i * n3
464 os.system("cp $FAIRSHIP/macro/run_simScript.py .")
465 os.system(cmd + " -n " + str(n3) + " -i " + str(ns) + " > log &")
466 # print " -n "+str(n3)+" -i "+str(ns)
467 ns += n3
468 if i == ncpu:
469 break
470 os.chdir("../" + prefix + str(i + 1))
471
472

◆ MergeRareEvents()

None ana_ShipMuon.MergeRareEvents ( list[str] | None   runs = None)

Definition at line 1242 of file ana_ShipMuon.py.

1242def MergeRareEvents(runs: list[str] | None = None) -> None:
1243 if runs is None:
1244 runs = ["61", "62"]
1245 for prefix in runs:
1246 cmd = "$ROOTSYS/bin/hadd rareEvents_" + prefix + ".root -f "
1247 for fn in fchainRec:
1248 if fn.find("muon" + prefix) < 0:
1249 continue
1250 fr = fn.replace(".root", "_rare.root")
1251 cmd = cmd + " " + fr
1252 os.system(cmd)
1253
1254
1255#

◆ muDISntuple()

None ana_ShipMuon.muDISntuple (   fn)

Definition at line 985 of file ana_ShipMuon.py.

985def muDISntuple(fn) -> None:
986 # take as input the rare events
987 fout = ROOT.TFile("muDISVetoCounter.root", "recreate")
988 h["ntuple"] = ROOT.TNtuple("muons", "muon flux VetoCounter", "id:px:py:pz:x:y:z:w")
989 f = ROOT.TFile(fn)
990 sTree = f.Get("cbmsim")
991 nEvents = sTree.GetEntries()
992 for n in range(nEvents):
993 sTree.GetEntry(n)
994 if len(sTree.MCTrack) > 1:
995 wg = sTree.MCTrack[1].GetWeight()
996 else:
997 wg = sTree.MCTrack[0].GetWeight()
998 for ahit in sTree.vetoPoint:
999 detID = ahit.GetDetectorID()
1000 if logVols[detID] != "VetoTimeDet":
1001 continue
1002 pid = ahit.PdgCode()
1003 if abs(pid) != 13:
1004 continue
1005 P = ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2 + ahit.GetPz() ** 2)
1006 if 3 / u.GeV < P:
1007 h["ntuple"].Fill(
1008 float(pid),
1009 float(ahit.GetPx() / u.GeV),
1010 float(ahit.GetPy() / u.GeV),
1011 float(ahit.GetPz() / u.GeV),
1012 float(ahit.GetX() / u.m),
1013 float(ahit.GetY() / u.m),
1014 float(ahit.GetZ() / u.m),
1015 float(wg),
1016 )
1017 fout.cd()
1018 h["ntuple"].Write()
1019
1020

◆ origin()

None ana_ShipMuon.origin (   sTree,
  it 
)

Definition at line 389 of file ana_ShipMuon.py.

389def origin(sTree, it) -> None:
390 at = sTree.MCTrack[it]
391 im = at.GetMotherId()
392 if im > 0:
393 origin(sTree, im)
394 if im < 0:
395 # print 'does not come from muon'
396 pass
397 if im == 0:
398 # print 'origin z',at.GetStartZ()
399 ROOT.TMath.Sqrt(at.GetStartX() ** 2 + at.GetStartY() ** 2), at.GetStartZ()
400
401

◆ originOfMuon()

None ana_ShipMuon.originOfMuon ( TextIOWrapper  fout,
int  n,
  fn,
  nEvents 
)

Definition at line 1359 of file ana_ShipMuon.py.

1359def originOfMuon(fout: TextIOWrapper, n: int, fn, nEvents) -> None:
1360 # from fn extract Yandex or CERN/cracow
1361 ncpu = 9
1362 x = fn.find("/")
1363 ni = int(fn[x - 1 : x]) - 1
1364 if nEvents < 100000:
1365 fm = "$SHIPSOFT/data/pythia8_Geant4_onlyMuons.root"
1366 else:
1367 fm = "$SHIPSOFT/data/pythia8_Geant4_Yandex_onlyMuons.root"
1368 fmuon = ROOT.TFile(fm)
1369 ntupl = fmuon.Get("pythia8-Geant4")
1370 ntot = ntupl.GetEntries()
1371 n3 = int(ntot / ncpu)
1372 N = n3 * ni + n
1373 ntupl.GetEvent(N)
1374 P = ROOT.TMath.Sqrt(ntupl.pz * ntupl.pz + ntupl.py * ntupl.py + ntupl.px * ntupl.px)
1375 fout.write("original muon %i, %i, %8.2F \n" % (ntupl.parentid, ntupl.pythiaid, P))
1376 fmuon.Close()
1377
1378
1379#
1380# BigEventLoop()
1381# makePlots()
1382# AnaEventLoop()
1383
1384
1385# ShipAna

◆ pers()

None ana_ShipMuon.pers ( )

Definition at line 1386 of file ana_ShipMuon.py.

1386def pers() -> None:
1387 xdisk = "/media/Work/HNL/"
1388 for x in h:
1389 if isinstance(h[x], ROOT.TCanvas):
1390 h[x].Update()
1391 tn = h[x].GetName() + ".png"
1392 h[x].Print(tn)
1393 rpath = os.path.abspath(".").split("/HNL/")[1]
1394 lp = rpath.split("/")
1395 prefix = xdisk
1396 for i in range(len(lp)):
1397 if lp[i] not in os.listdir(prefix):
1398 os.system("mkdir " + prefix + "/" + lp[i])
1399 prefix = prefix + "/" + lp[i]
1400 os.system("cp " + tn + " " + xdisk + rpath)
1401
1402

◆ persistency()

None ana_ShipMuon.persistency ( )

Definition at line 1256 of file ana_ShipMuon.py.

1256def persistency() -> None:
1257 printAndCopy(prefix)
1258 ut.writeHists(h, prefix + ".root", plusCanvas=True)
1259
1260

◆ printAndCopy()

None ana_ShipMuon.printAndCopy ( str | None   prefix = None)

Definition at line 1271 of file ana_ShipMuon.py.

1271def printAndCopy(prefix: str | None = None) -> None:
1272 if not prefix:
1273 prefix = (h["tc"].GetName()).replace(".root", "")
1274 for x in ["ResultsI", "ResultsII", "ResultsImu", "ResultsImuV0", "ResultsIII", "ResultsIV", "ResultsV"]:
1275 h[x].Update()
1276 if prefix not in os.listdir("."):
1277 os.mkdir(prefix)
1278 os.chdir(prefix)
1279 h["ResultsI"].Print(prefix + "Back_occ.png")
1280 h["ResultsII"].Print(prefix + "Back_depE.png")
1281 h["ResultsImu"].Print(prefix + "muBack_occ.png")
1282 h["ResultsImuV0"].Print(prefix + "muV0Back_occ.png")
1283 h["ResultsIII"].Print(prefix + "Back_P.png")
1284 h["ResultsIV"].Print(prefix + "Back_OP.png")
1285 h["ResultsV"].Print(prefix + "origin.png")
1286 os.chdir("../")
1287
1288

◆ rareEventEmulsion()

None ana_ShipMuon.rareEventEmulsion ( str   fname = "rareEmulsion.txt")

Definition at line 1099 of file ana_ShipMuon.py.

1099def rareEventEmulsion(fname: str = "rareEmulsion.txt") -> None:
1100 fout = open(fname, "w") # noqa: SIM115
1101 for fn in fchainRec:
1102 f = ROOT.TFile(fn)
1103 if not f.FindObjectAny("cbmsim"):
1104 print("skip file ", f.GetName())
1105 continue
1106 sTree = f.Get("cbmsim")
1107 sTree.GetEvent(0)
1108 if sTree.GetBranch("GeoTracks"):
1109 sTree.SetBranchStatus("GeoTracks", 0)
1110 nEvents = sTree.GetEntries()
1111 for n in range(nEvents):
1112 sTree.GetEntry(n)
1113 if n == 0:
1114 print("now at event ", n, f.GetName())
1115 for ahit in sTree.vetoPoint:
1116 detID = ahit.GetDetectorID()
1117 if logVols[detID] != "Emulsion":
1118 continue
1119 ahit.GetX()
1120 ahit.GetY()
1121 ahit.GetZ()
1122 if len(sTree.MCTrack) > 1:
1123 wg = sTree.MCTrack[1].GetWeight() # also works for neutrinos
1124 else:
1125 wg = sTree.MCTrack[0].GetWeight() # also works for neutrinos
1126 fout.write("rare emulsion hit %i, %s, %8.3F, %i \n" % (n, f.GetName(), wg, ahit.PdgCode()))
1127 if ahit.GetPz() / u.GeV > 1.0:
1128 fout.write(
1129 "V,P when making hit %8.3F,%8.3F,%8.3F %8.3F,%8.3F,%8.3F (GeV) \n"
1130 % (
1131 ahit.GetX() / u.m,
1132 ahit.GetY() / u.m,
1133 ahit.GetZ() / u.m,
1134 ahit.GetPx() / u.GeV,
1135 ahit.GetPy() / u.GeV,
1136 ahit.GetPz() / u.GeV,
1137 )
1138 )
1139 else:
1140 fout.write(
1141 "V,P when making hit %8.3F,%8.3F,%8.3F %8.3F,%8.3F,%8.3F (MeV)\n"
1142 % (
1143 ahit.GetX() / u.m,
1144 ahit.GetY() / u.m,
1145 ahit.GetZ() / u.m,
1146 ahit.GetPx() / u.MeV,
1147 ahit.GetPy() / u.MeV,
1148 ahit.GetPz() / u.MeV,
1149 )
1150 )
1151 originOfMuon(fout, n, f.GetName(), nEvents)
1152
1153
1154#

◆ readAndMergeHistos()

None ana_ShipMuon.readAndMergeHistos (   prods)

Definition at line 1486 of file ana_ShipMuon.py.

1486def readAndMergeHistos(prods) -> None:
1487 for p in prods:
1488 x = p
1489 if p.find(".root") < 0:
1490 x = p + ".root"
1491 ut.readHists(h, x)
1492 # make list of hists with entries
1493 k = 1
1494 for x in histlistAll:
1495 if histlistAll[x] in h:
1496 histlist[k] = histlistAll[x]
1497 k += 1
1498 nstations = len(histlist)
1499 print("make plots for ", nstations)
1500 makePlots(nstations)
1501 printAndCopy(prods[0].replace(".root", ""))
1502
1503
1504# python -i $HNL/ana_ShipMuon.py 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829
1505# python -i $HNL/ana_ShipMuon.py 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929
1506# python -i $HNL/ana_ShipMuon.py 1012 1013 1014 1015 1016 1017 1018 1019 1022 1023 1024 1025 1026 1027 1028 1029
1507# python -i $HNL/ana_ShipMuon.py 634 635 636 637 638 639 640 641 642 643
1508# combine all cern-cracow
1509# python -i $HNL/ana_ShipMuon.py 810 811 812 813 814 815 816 817 818 819 910 911 912 913 914 915 916 917 918 919 1012 1013 1014 1015 1016 1017 1018 1019
1510# python -i $HNL/ana_ShipMuon.py 820 821 822 823 824 825 826 827 828 829 920 921 922 923 924 925 926 927 928 929 1022 1023 1024 1025 1026 1027 1028 1029
1511# make muonDIS ntuple: muDISntuple("/media/Data/HNL/muonBackground/rareEvents_81-102.root") -> 'muDISVetoCounter.root'
1512# second step python $FAIRSHIP/muonShieldOptimization/makeMuonDIS.py 1 10000 muDISVetoCounter.root
1513# third step run_simScript.py --MuDIS -n 10 -f muonDis_1.root
1514# for concrete
1515# analyzeConcrete() -> muConcrete.root

◆ reDraw()

None ana_ShipMuon.reDraw (   fn)

Definition at line 1261 of file ana_ShipMuon.py.

1261def reDraw(fn) -> None:
1262 if fn.find("root") < 0:
1263 fn = fn + ".root"
1264 if "tc" not in h:
1265 h["tc"] = ROOT.TFile(fn)
1266 for x in ["ResultsI", "ResultsII", "ResultsImu", "ResultsImuV0", "ResultsIII", "ResultsIV", "ResultsV"]:
1267 h[x] = h["tc"].FindObjectAny(x)
1268 h[x].Draw()
1269
1270

Variable Documentation

◆ dy

float ana_ShipMuon.dy = float(tmp[1] + "." + tmp[2])

Definition at line 353 of file ana_ShipMuon.py.

◆ else

ana_ShipMuon.else :

Definition at line 335 of file ana_ShipMuon.py.

◆ fchain

list ana_ShipMuon.fchain = []

Definition at line 413 of file ana_ShipMuon.py.

◆ fchainRec

list ana_ShipMuon.fchainRec = []

Definition at line 414 of file ana_ShipMuon.py.

◆ fgeo

ROOT ana_ShipMuon.fgeo = ROOT.TFile(testdir + "/" + f)

Definition at line 346 of file ana_ShipMuon.py.

◆ fname

str ana_ShipMuon.fname = path + prefix + str(i) + "/" + inputFile

Definition at line 436 of file ana_ShipMuon.py.

◆ h

dict ana_ShipMuon.h = {}

Definition at line 506 of file ana_ShipMuon.py.

◆ histlist

dict ana_ShipMuon.histlist = {}

Definition at line 507 of file ana_ShipMuon.py.

◆ histlistAll

dict ana_ShipMuon.histlistAll
Initial value:
1= {
2 1: "strawstation_5",
3 2: "strawstation_1",
4 3: "strawstation_4",
5 4: "Ecal",
6 5: "muondet",
7 6: "VetoTimeDet",
8 7: "T1LiSc",
9 8: "T2LiSc",
10 9: "T3LiSc",
11 10: "T5LiSc",
12 11: "ShipRpc",
13 12: "volHPT",
14 13: "Target",
15 14: "TimeDet",
16 15: "Det2",
17}

Definition at line 546 of file ana_ShipMuon.py.

◆ hLiSc

dict ana_ShipMuon.hLiSc = {1: {}}

Definition at line 563 of file ana_ShipMuon.py.

◆ hMuon

dict ana_ShipMuon.hMuon = {}

Definition at line 575 of file ana_ShipMuon.py.

◆ inputFile

f ana_ShipMuon.inputFile = f.replace("geofile_full", "ship")

Definition at line 348 of file ana_ShipMuon.py.

◆ inputFile1

f ana_ShipMuon.inputFile1 = inputFile.replace("_D.", ".")

Definition at line 358 of file ana_ShipMuon.py.

◆ inputFile2

f ana_ShipMuon.inputFile2 = inputFile

Definition at line 359 of file ana_ShipMuon.py.

◆ labels

dict ana_ShipMuon.labels = {}

Definition at line 509 of file ana_ShipMuon.py.

◆ local

bool ana_ShipMuon.local = False

Definition at line 13 of file ana_ShipMuon.py.

◆ logVols

def ana_ShipMuon.logVols = detMap()

Definition at line 510 of file ana_ShipMuon.py.

◆ modules

shipDet_conf ana_ShipMuon.modules = shipDet_conf.configure(run, ShipGeo)

Definition at line 383 of file ana_ShipMuon.py.

◆ mom

ROOT ana_ShipMuon.mom = ROOT.TVector3()

Definition at line 581 of file ana_ShipMuon.py.

◆ muon

ROOT ana_ShipMuon.muon = top.GetNode("MuonDetector_1")

Definition at line 407 of file ana_ShipMuon.py.

◆ mvol

ROOT ana_ShipMuon.mvol = muon.GetVolume()

Definition at line 408 of file ana_ShipMuon.py.

◆ noField

bool ana_ShipMuon.noField = False

Definition at line 403 of file ana_ShipMuon.py.

◆ otherPhysList

bool ana_ShipMuon.otherPhysList = False

Definition at line 402 of file ana_ShipMuon.py.

◆ output

mp ana_ShipMuon.output = mp.Queue()

Definition at line 20 of file ana_ShipMuon.py.

◆ parallel

bool ana_ShipMuon.parallel = True

Definition at line 17 of file ana_ShipMuon.py.

◆ passive

bool ana_ShipMuon.passive = False

Definition at line 404 of file ana_ShipMuon.py.

◆ path

str ana_ShipMuon.path = ""

Definition at line 339 of file ana_ShipMuon.py.

◆ PDG

ROOT ana_ShipMuon.PDG = ROOT.TDatabasePDG.Instance()

Definition at line 367 of file ana_ShipMuon.py.

◆ pos

ROOT ana_ShipMuon.pos = ROOT.TVector3()

Definition at line 582 of file ana_ShipMuon.py.

◆ pref

str ana_ShipMuon.pref = "muon"

Definition at line 317 of file ana_ShipMuon.py.

◆ prefix

str ana_ShipMuon.prefix = pref + prefix

Definition at line 332 of file ana_ShipMuon.py.

◆ prefixes

list ana_ShipMuon.prefixes = []

makeProd("muon812",10,False,True) # –< 831 copied back, done 16.3.2015 makeProd("muon822",10,True,True) makeProd("muon821",10,True,True) makeProd("muon822",10,True,True)

makeProd("muon813",10,False,True) # makeProd("muon823",10,True,True)

Definition at line 315 of file ana_ShipMuon.py.

◆ processes

list ana_ShipMuon.processes = []

Definition at line 21 of file ana_ShipMuon.py.

◆ q1

f ana_ShipMuon.q1 = inputFile1 in os.listdir(path + prefix + str(i))

Definition at line 424 of file ana_ShipMuon.py.

◆ q2

f ana_ShipMuon.q2 = inputFile2 in os.listdir(path + prefix + str(i))

Definition at line 425 of file ana_ShipMuon.py.

◆ r1

f ana_ShipMuon.r1 = recFile1 in os.listdir(path + prefix + str(i))

Definition at line 428 of file ana_ShipMuon.py.

◆ r2

f ana_ShipMuon.r2 = recFile2 in os.listdir(path + prefix + str(i))

Definition at line 429 of file ana_ShipMuon.py.

◆ recFile

f ana_ShipMuon.recFile = inputFile.replace(".root", "_rec.root")

Definition at line 437 of file ana_ShipMuon.py.

◆ recFile1

f ana_ShipMuon.recFile1 = inputFile1.replace(".root", "_rec.root")

Definition at line 426 of file ana_ShipMuon.py.

◆ recFile2

f ana_ShipMuon.recFile2 = inputFile2.replace(".root", "_rec.root")

Definition at line 427 of file ana_ShipMuon.py.

◆ run

ROOT ana_ShipMuon.run = ROOT.FairRunSim()

Definition at line 382 of file ana_ShipMuon.py.

◆ rz_inter

float ana_ShipMuon.rz_inter = -1.0, 0.0

Definition at line 386 of file ana_ShipMuon.py.

◆ sGeo

ROOT ana_ShipMuon.sGeo = fgeo.Get("FAIRGeom")

Definition at line 347 of file ana_ShipMuon.py.

◆ ShipGeo

load_from_root_file ana_ShipMuon.ShipGeo = geometry_config.create_config(Yheight=dy)

Definition at line 374 of file ana_ShipMuon.py.

◆ testdir

str ana_ShipMuon.testdir = "."

Definition at line 338 of file ana_ShipMuon.py.

◆ tmp

f ana_ShipMuon.tmp = inputFile.split(".")

Definition at line 351 of file ana_ShipMuon.py.

◆ top

ROOT ana_ShipMuon.top = sGeo.GetTopVolume()

Definition at line 406 of file ana_ShipMuon.py.

◆ totl

tuple ana_ShipMuon.totl = (zmuon + mvol.GetShape().GetDZ()) / u.m

Definition at line 410 of file ana_ShipMuon.py.

◆ try

ana_ShipMuon.try :

Definition at line 352 of file ana_ShipMuon.py.

◆ txt

dict ana_ShipMuon.txt = {}

Definition at line 508 of file ana_ShipMuon.py.

◆ withChain

int ana_ShipMuon.withChain = 0

Definition at line 316 of file ana_ShipMuon.py.

◆ xx

os ana_ShipMuon.xx = os.path.abspath(".").lower()

Definition at line 318 of file ana_ShipMuon.py.

◆ zmuon

ROOT ana_ShipMuon.zmuon = muon.GetMatrix().GetTranslation()[2]

Definition at line 409 of file ana_ShipMuon.py.

◆ ztarget

float ana_ShipMuon.ztarget = -100.0

Definition at line 411 of file ana_ShipMuon.py.