ana_ShipMuon.py

Go to the documentation of this file.

# SPDX-License-Identifier: LGPL-3.0-or-later
# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
 
# analyze muon background /media/Data/HNL/PythiaGeant4Production/pythia8_Geant4_total.root
import multiprocessing as mp
import os
from io import TextIOWrapper
 
import ROOT
from ShipGeoConfig import load_from_root_file
 
ROOT.gInterpreter.ProcessLine("typedef double Double32_t")
local = False
if not os.uname()[1].lower().find("ubuntu") < 0:
    local = True
 
parallel = True
if parallel:
    # Define an output queue
    output = mp.Queue()
    processes = []
 
 
# 11-19 with QGSP_BERT_EMV instead of QGSP_BERT_HP_PEN
# 51-59 passive shielding
# 41-49 fixed field horizontal
# 61-69  chamber with vacuum
# 71-79 without field in the wings
# 81-89 add shielding
# 91-99 switch off  muBrems, muPair
# 101-109 vacuum tube as cone, Al->Steel
 
# prefix = 'muon4' # default
# prefix = 'muon6' # Al -> vacuum
# prefix = 'muon7' # no field in wings
# prefix = 'muon8' # entry window + additional shielding
# prefix = 'muon9' # switch off  muBrems, muPair
# prefix = 'muon10' # vacuum tube as cone, Al->Steel
# prefix = 'muon14' # vacuum tube as cone, Al->Steel, slightly smaller lead shield
# prefix = 'muon12' # switch off  muIoni
# prefix = 'muon13' # switch off  muBrems, muPair AND muIoni
# prefix = 'muon15' # vacuum tube as cone, Al->Steel, cave with air
# prefix = 'muon 161-169 new geometry
# prefix = 'muon 251-259 passive shielding with concrete walls
# prefix = 'muon 171-179 new increased Bdl,  2.5m clearance from start, 9m space for tau
# 181-189 narrow tunnel around first magnet, add 2m hadron absorber, fixed problem with B-field -> G4
# 191-199 narrow tunnel bug fix, also bug fix for top/bottom, now with tungsten core in add absorber
# 201-209 test with reduced field, 1.5T
# 211-219 many bug fixes, overlaps, added Goliath magnet
# 221-229 bug fix for Goliath
# 231-239 reduce field to 1T
# 311-314 change numbering of veto volumes
# 321-324 1 Teslas
# 331-339 added TauSensitive, fixed C1, using muons from Yandex
# 341-349 should be the same as 331-339
# 361-364 4m high magnets!
# 371-374 new strawpoints
# 381-389 yandex file, overwrite with magC1 fixed
# 391-399 change to dY = 2m
# 401-409 change to dY = 2m Yandex
# 411-419 back to dY = 4m, add sensitive plance downstream det2
# 421-429 reduce size of vacuum tank before veto station, add pdgcode for muonPoint, move det2 closer
# 431-439 reduce size of vacuum tank before veto station, add pdgcode for muonPoint, move det2 closer Yandex
# 441-449 try with increased tau mu detector, and reduced height of active muon shield, dy=1m
# 451-459 try with increased tau mu detector, and reduced height of active muon shield, dy=1m  Yandex
# 461-469 dy=4m and then 1m, RPC fixed, display on
# 471-479 dy=4m and then 1m, RPC fixed, Yandex
# 481-489 dy=4m and then 1m, RPC fixed, display OFF
# 491-499 back to dy=1m display on
# 551-554 ultimate slim design
# 561-564 fixing overlaps hcal and almost overlaps in magnets
# 571-574 more undetected overlaps fixed
# 581-584 test with replacing tungsten core with iron, 10m height
# 591-599 test with replacing tungsten core with iron, 6m height, display
# 601-609 replacing tungsten core with iron, 6m height, Yandex data, display
# 610-619 replacing tungsten core with iron, 10m height, display
# 620-629 replacing tungsten core with iron, 10m height, Yandex data, display
# 630-639 10m height, change RPC width 4.5->3.6m
# 640-649 10m height, change RPC width 4.5->3.6m Yandex
# 650-669 6m height, change RPC width 4.5->3.6m
# 660-669 6m height, change RPC width 4.5->3.6m Yandex
# makeProd("muon700",10,False,False) # switch off field of active shielding    # prefix,DY,y=False,phiRandom=False,X=None
# makeProd("muon710",10,False,False) # start production with beam smeared on r=3cm disk
# makeProd("muon720",10,True,False)
# makeProd("muon711",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon721",10,True,True)
# makeProd("muon712",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon722",10,True,True)
# makeProd("muon713",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon723",10,True,True)
# makeProd("muon714",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon724",10,True,True)
# makeProd("muon715",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon725",10,True,True)
# makeProd("muon716",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon726",10,True,True)
# makeProd("muon717",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon777",10,False,True) # in case the other one does not work 717
# makeProd("muon727",10,True,True)  # ?
# makeProd("muon718",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon728",10,True,True)
# makeProd("muon719",10,False,True) # start production with beam smeared on r=3cm disk
# makeProd("muon729",10,True,True)
# makeProd("muon730",10,'Jpsi',False)
# makeProd("muon731",10,'Jpsi',True)
# makeProd("muon732",10,'Jpsi',True)
# makeProd("muon733",10,'Jpsi',True)   # back to pencil beam
# makeProd("muon630",10,False,True) # test with new muonShield code, 3cm smearing
# makeProd("muon631",10,False,True) # run with concrete wall enabled as sensitive
# makeProd("muon632",10,False,True) # run with concrete wall enabled as sensitive, active shielding polarity fixed
# but wrong geometry
# makeProd("muon810",10,False,False) # start production with latest geometry
# makeProd("muon820",10,True,False)
# makeProd("muon811",10,False,True) #
# makeProd("muon821",10,True,True)
 
 
# makeProd("muon814",10,False,True) #
# makeProd("muon824",10,True,True)
# makeProd("muon815",10,False,True)
# makeProd("muon825",10,True,True)
# makeProd("muon816",10,False,True)
# makeProd("muon826",10,True,True)
# makeProd("muon817",10,False,True)
# makeProd("muon827",10,True,True)
# makeProd("muon818",10,False,True)
# makeProd("muon828",10,True,True)
# makeProd("muon819",10,False,True)
# makeProd("muon829",10,True,True)
# makeProd("muon910",10,False,True)
# makeProd("muon920",10,True,True)
# makeProd("muon911",10,False,True)
# makeProd("muon921",10,True,True)
 
# makeProd("muon912",10,False,True)
# makeProd("muon922",10,True,True)
 
# makeProd("muon913",10,False,True)
# makeProd("muon923",10,True,True)
# makeProd("muon914",10,False,True)
# makeProd("muon924",10,True,True)
# makeProd("muon915",10,False,True)
# makeProd("muon925",10,True,True)
# makeProd("muon916",10,False,True)
# makeProd("muon926",10,True,True)
# makeProd("muon917",10,False,True)
# makeProd("muon927",10,True,True)
# makeProd("muon927",10,True,True,8)
# makeProd("muon918",10,False,True)
# makeProd("muon928",10,True,True)
# makeProd("muon919",10,False,True)
# makeProd("muon929",10,True,True)
# makeProd("muon1019",10,False,True)
# makeProd("muon1029",10,True,True)
# makeProd("muon1018",10,False,True)
# makeProd("muon1028",10,True,True)
# makeProd("muon1017",10,False,True)
# makeProd("muon1027",10,True,True)
# makeProd("muon1016",10,False,True)
# makeProd("muon1026",10,True,True)
# makeProd("muon1015",10,False,True)
# makeProd("muon1025",10,True,True)
# makeProd("muon1014",10,False,True)
# makeProd("muon1024",10,True,True)
# makeProd("muon1013",10,False,True)
# makeProd("muon1023",10,True,True)
# makeProd("muon1012",10,False,True)
# makeProd("muon1022",10,True,True)
# makeProd("muon633",10,'concrete',False) # run with concrete wall enabled as sensitive
# makeProd("muon1111",10,'False',False) # try iron for first shield
# makeProd("muon1121",10,'True',False)  # try iron for first shield
# makeProd("muon1112",10,'False',True) # try iron for first shield
# makeProd("muon1122",10,'True',True)  # try iron for first shield
# restart with also rockS sensitive and &&->||
# makeProd("muon634",10,'concrete',False) # run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon635",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon636",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon637",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon638",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon639",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon640",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon641",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon642",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon643",10,'concrete',True) #   run with concrete wall enabled as sensitive, and active shield too
# makeProd("muon650",10,'concrete',False) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon651",10,'concrete',False) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon652",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon653",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon654",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon655",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon656",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon657",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon658",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon659",10,'concrete',True) # run with concrete wall enabled as sensitive, first tunnel 10m wide
# makeProd("muon660",10,'concrete',False) # liter magnet
# makeProd("muon661",10,'concrete',True) # liter magnet
# makeProd("muon662",10,'concrete',True) # liter magnet
# makeProd("muon663",10,'concrete',True) # liter magnet
# makeProd("muon664",10,'concrete',True) # liter magnet
# makeProd("muon665",10,'concrete',True) # liter magnet
# makeProd("muon666",10,'concrete',True) # liter magnet
# makeProd("muon667",10,'concrete',True) # liter magnet
# makeProd("muon668",10,'concrete',True) # liter magnet
# makeProd("muon669",10,'concrete',True) # liter magnet
# makeProd("muon670",10,'concrete',False) # even liter magnet
# makeProd("muon671",10,'concrete',True) # even liter magnet
# makeProd("muon672",10,'concrete',True) # even liter magnet
# makeProd("muon673",10,'concrete',True) # even liter magnet
# makeProd("muon674",10,'concrete',True) # even liter magnet
# makeProd("muon675",10,'concrete',True) # even liter magnet
# makeProd("muon676",10,'concrete',True) # even liter magnet
# makeProd("muon677",10,'concrete',True) # even liter magnet
# makeProd("muon678",10,'concrete',True) # even liter magnet
# makeProd("muon679",10,'concrete',True) # even liter magnet
# new attempt, previous did not worked out
# makeProd("muon680",10,'concrete',False) # even liter magnet
# makeProd("muon681",10,'concrete',True) # even liter magnet
# makeProd("muon682",10,'concrete',True) # even liter magnet
# makeProd("muon683",10,'concrete',True) # even liter magnet
# makeProd("muon684",10,'concrete',True) # even liter magnet
# makeProd("muon685",10,'concrete',True) # even liter magnet
# makeProd("muon686",10,'concrete',True) # even liter magnet
# makeProd("muon687",10,'concrete',True) # even liter magnet
# makeProd("muon688",10,'concrete',True) # even liter magnet
# makeProd("muon689",10,'concrete',True) # even liter magnet
# new attempt, increase height of first magnet, previous did not worked out either
# makeProd("muon690",10,'concrete',False) # even liter magnet
# makeProd("muon691",10,'concrete',True) # even liter magnet
# makeProd("muon692",10,'concrete',True) # even liter magnet
# makeProd("muon693",10,'concrete',True) # even liter magnet
# makeProd("muon694",10,'concrete',True) # even liter magnet
# makeProd("muon695",10,'concrete',True) # even liter magnet
# makeProd("muon696",10,'concrete',True) # even liter magnet
# makeProd("muon697",10,'concrete',True) # even liter magnet
# makeProd("muon698",10,'concrete',True) # even liter magnet
# makeProd("muon699",10,'concrete',True) # even liter magnet
# testing height 4m->2m
# makeProd("muon700",10,'concrete',False) # even liter magnet
# makeProd("muon701",10,'concrete',True) # even liter magnet
# makeProd("muon702",10,'concrete',True) # even liter magnet
# makeProd("muon703",10,'concrete',True) # even liter magnet
# makeProd("muon704",10,'concrete',True) # even liter magnet
# makeProd("muon705",10,'concrete',True) # even liter magnet
# makeProd("muon706",10,'concrete',True) # even liter magnet
# makeProd("muon707",10,'concrete',True) # even liter magnet
# makeProd("muon708",10,'concrete',True) # even liter magnet
# makeProd("muon709",10,'concrete',True) # even liter magnet
# testing height 4m->2m + new magnets
# makeProd("muon800",10,'concrete',False) # even liter magnet
# makeProd("muon801",10,'concrete',True) # even liter magnet
# makeProd("muon802",10,'concrete',True) # even liter magnet
# makeProd("muon803",10,'concrete',True) # even liter magnet
# makeProd("muon804",10,'concrete',True) # even liter magnet
# makeProd("muon805",10,'concrete',True) # even liter magnet
# makeProd("muon806",10,'concrete',True) # even liter magnet
# makeProd("muon807",10,'concrete',True) # even liter magnet
# makeProd("muon808",10,'concrete',True) # even liter magnet
# makeProd("muon809",10,'concrete',True) # even liter magnet
# new iteration, and ship_geo.Yheight*1./10.
# makeProd("muon710",10,'concrete',False) # even liter magnet
# makeProd("muon711",10,'concrete',True) # even liter magnet
# makeProd("muon712",10,'concrete',True) # even liter magnet
# makeProd("muon713",10,'concrete',True) # even liter magnet
# makeProd("muon714",10,'concrete',True) # even liter magnet
# makeProd("muon715",10,'concrete',True) # even liter magnet
# makeProd("muon716",10,'concrete',True) # even liter magnet
# makeProd("muon717",10,'concrete',True) # even liter magnet
# makeProd("muon718",10,'concrete',True) # even liter magnet
# makeProd("muon719",10,'concrete',True) # even liter magnet
# new iteration, and back to ship_geo.Yheight*2./10.
# makeProd("muon720",10,'concrete',False) # even liter magnet
# makeProd("muon721",10,'concrete',True) # even liter magnet
# makeProd("muon722",10,'concrete',True) # even liter magnet
# makeProd("muon723",10,'concrete',True) # even liter magnet
# makeProd("muon724",10,'concrete',True) # even liter magnet
# makeProd("muon725",10,'concrete',True) # even liter magnet
# makeProd("muon726",10,'concrete',True) # even liter magnet
# makeProd("muon727",10,'concrete',True) # even liter magnet
# makeProd("muon728",10,'concrete',True) # even liter magnet
# makeProd("muon729",10,'concrete',True) # even liter magnet
# new iteration, and back to ship_geo.Yheight*2./10. but with more info in vetoPoint
# makeProd("muon730",10,'concrete',False) # even liter magnet
# makeProd("muon740",10,'concrete',False) # even liter magnet
# new iteration,  ship_geo.Yheight*1./10.
# makeProd("muon750",10,'concrete',False) #
# new iteration,  ship_geo.Yheight*1.5/10.
# makeProd("muon760",10,'concrete',False) #
# makeProd("muon761",10,'concrete',True) #
# as before but now with full simulation
# makeProd("muon1710",10,False,False)
# makeProd("muon1720",10,True,False)
# makeProd("muon1711",10,False,True)
# makeProd("muon1721",10,True,True)
# makeProd("muon1712",10,False,True)
# makeProd("muon1722",10,True,True)
# makeProd("muon1713",10,False,True)
# makeProd("muon1723",10,True,True)
# makeProd("muon1714",10,False,True)
# makeProd("muon1724",10,True,True)
# makeProd("muon1715",10,False,True)
# makeProd("muon1725",10,True,True)
# makeProd("muon1716",10,False,True)
# makeProd("muon1726",10,True,True)
# another one, hopefully better
# makeProd("muon1717",10,'concrete',False) #
# another one, increasing height to 3m
# makeProd("muon1718",10,'concrete',False) #
 
prefixes = []
withChain = 0
pref = "muon"
xx = os.path.abspath(".").lower()
if not xx.find("neutrino") < 0:
    pref = "neutrino"
if not xx.find("vdis") < 0:
    pref = "disV"
elif not xx.find("clby") < 0:
    pref = "disCLBY"
elif not xx.find("dis") < 0:
    pref = "dis"
 
if len(os.sys.argv) > 1:
    for i in range(1, len(os.sys.argv)):
        for prefix in os.sys.argv[i].split(","):
            if prefix.find(pref) < 0:
                prefix = pref + prefix
            prefixes.append(prefix)
            withChain += 1
else:
    prefixes = [""]
 
testdir = "."
path = ""
# if local: path = "/media/Data/HNL/muonBackground/"
if prefixes[0] != "":
    testdir = path + prefixes[0] + "1"
# figure out which setup
for f in os.listdir(testdir):
    if not f.find("geofile_full") < 0:
        fgeo = ROOT.TFile(testdir + "/" + f)
        sGeo = fgeo.Get("FAIRGeom")
        inputFile = f.replace("geofile_full", "ship")
        break
# try to extract from input file name
tmp = inputFile.split(".")
try:
    dy = float(tmp[1] + "." + tmp[2])
except Exception:
    dy = 10.0
 
if not inputFile.find("_D.") < 0:
    inputFile1 = inputFile.replace("_D.", ".")
    inputFile2 = inputFile
else:
    inputFile1 = inputFile
    inputFile2 = inputFile.replace(".root", "_D.root")
 
import rootUtils as ut
import shipunit as u
 
PDG = ROOT.TDatabasePDG.Instance()
import geometry_config
 
# init geometry and mag. field
if not fgeo.FindKey("ShipGeo"):
    # old geofile, missing Shipgeo dictionary
    # re-create geometry and mag. field
    ShipGeo = geometry_config.create_config(Yheight=dy)
else:
    # new geofile, load Shipgeo dictionary written by run_simScript.py
    ShipGeo = load_from_root_file(fgeo, "ShipGeo")
 
# -----Create geometry----------------------------------------------
import shipDet_conf
 
run = ROOT.FairRunSim()
modules = shipDet_conf.configure(run, ShipGeo)
 
 
rz_inter = -1.0, 0.0
 
 
def origin(sTree, it) -> None:
    at = sTree.MCTrack[it]
    im = at.GetMotherId()
    if im > 0:
        origin(sTree, im)
    if im < 0:
        # print 'does not come from muon'
        pass
    if im == 0:
        # print 'origin z',at.GetStartZ()
        ROOT.TMath.Sqrt(at.GetStartX() ** 2 + at.GetStartY() ** 2), at.GetStartZ()
 
 
otherPhysList = False
noField = False
passive = False
#
top = sGeo.GetTopVolume()
muon = top.GetNode("MuonDetector_1")
mvol = muon.GetVolume()
zmuon = muon.GetMatrix().GetTranslation()[2]
totl = (zmuon + mvol.GetShape().GetDZ()) / u.m
ztarget = -100.0
 
fchain = []
fchainRec = []
# first time trying to read a chain of Fairship files, doesn't seem to work out of the box, try some tricks.
testdir = "."
if path != "":
    testdir = path
if withChain > 0:
    for prefix in prefixes:
        for i in range(1, 10):
            if prefix + str(i) not in os.listdir(testdir):
                continue
            q1 = inputFile1 in os.listdir(path + prefix + str(i))
            q2 = inputFile2 in os.listdir(path + prefix + str(i))
            recFile1 = inputFile1.replace(".root", "_rec.root")
            recFile2 = inputFile2.replace(".root", "_rec.root")
            r1 = recFile1 in os.listdir(path + prefix + str(i))
            r2 = recFile2 in os.listdir(path + prefix + str(i))
            if q1 or r1:
                inputFile = inputFile1
            elif q2 or r2:
                inputFile = inputFile2
            else:
                continue
            fname = path + prefix + str(i) + "/" + inputFile
            recFile = inputFile.replace(".root", "_rec.root")
            if recFile not in os.listdir(path + prefix + str(i)):
                fchain.append(fname)
                continue
            fname = path + prefix + str(i) + "/" + recFile
            fchainRec.append(fname)
            fchain.append(fname)  # take rec file if exist
else:
    fchain.append(inputFile)
    prefix = ""
 
 
def makeProd() -> None:
    ntot = 736406
    ncpu = 4
    n3 = int(ntot / ncpu)
    cmd = "python $FAIRSHIP/macro/run_simScript.py --MuonBack -f $SHIPSOFT/data/pythia8_Geant4_onlyMuons.root "  # --display"
    ns = 0
    prefix = "muon18"
    for i in range(1, ncpu + 1):
        d = prefix + str(i)
        if d not in os.listdir("."):
            os.system("mkdir " + d)
    os.chdir("./" + prefix + "1")
    for i in range(1, ncpu + 1):
        if i == ncpu:
            n3 = ntot - i * n3
        os.system("cp $FAIRSHIP/macro/run_simScript.py .")
        os.system(cmd + " -n " + str(n3) + " -i " + str(ns) + " > log &")
        # print " -n "+str(n3)+" -i "+str(ns)
        ns += n3
        if i == ncpu:
            break
        os.chdir("../" + prefix + str(i + 1))
 
 
def detMap():
    sGeo = ROOT.gGeoManager
    detList = {}
    for v in sGeo.GetListOfVolumes():
        nm = v.GetName()
        i = sGeo.FindVolumeFast(nm).GetNumber()
        detList[i] = nm
    return detList
 
 
def fitSingleGauss(x, ba=None, be=None) -> None:
    name = "myGauss_" + x
    myGauss = h[x].GetListOfFunctions().FindObject(name)
    if not myGauss:
        if not ba:
            ba = h[x].GetBinCenter(1)
        if not be:
            be = h[x].GetBinCenter(h[x].GetNbinsX())
        bw = h[x].GetBinWidth(1)
        mean = h[x].GetMean()
        sigma = h[x].GetRMS()
        norm = h[x].GetEntries() * 0.3
        myGauss = ROOT.TF1(name, "[0]*" + str(bw) + "/([2]*sqrt(2*pi))*exp(-0.5*((x-[1])/[2])**2)+[3]", 4)
        myGauss.SetParameter(0, norm)
        myGauss.SetParameter(1, mean)
        myGauss.SetParameter(2, sigma)
        myGauss.SetParameter(3, 1.0)
        myGauss.SetParName(0, "Signal")
        myGauss.SetParName(1, "Mean")
        myGauss.SetParName(2, "Sigma")
    h[x].Fit(myGauss, "", "", ba, be)
 
 
h = {}
histlist = {}
txt = {}
labels = {}
logVols = detMap()
 
 
def bookHist(detName) -> None:
    for mu in ["", "_mu", "_muV0"]:
        tag = detName + mu
        if detName.find("LS") < 0:
            ut.bookHist(h, tag, "x/y " + detName, 100, -3.0, 3.0, 100, -6.0, 6.0)
        else:
            ut.bookHist(h, tag, "z phi " + detName, 100, -50.0, 50.0, 100, -1.0, 1.0)
        ut.bookHist(h, tag + "_E", "deposited energy MeV " + detName, 100, 0.0, 2.0)
        ut.bookHist(h, tag + "_P", "particle mom GeV " + detName, 500, 0.0, 50.0)
        ut.bookHist(h, tag + "_LP", "particle mom GeV " + detName, 100, 0.0, 1.0)
        ut.bookHist(h, tag + "_OP", "original particle mom GeV " + detName, 400, 0.0, 400.0)
        ut.bookHist(h, tag + "_id", "particle id " + detName, 5001, -2499.5, 2499.5)
        ut.bookHist(h, tag + "_mul", "multiplicity of hits/tracks " + detName, 100, -0.5, 99.5)
        ut.bookHist(h, tag + "_evmul", "multiplicity of hits/event " + detName, 100, -0.5, 9999.5)
        ut.bookHist(h, tag + "_origin", "r vs z", 100, ztarget, totl, 100, 0.0, 12.0)
        ut.bookHist(h, tag + "_originmu", "r vs z", 100, ztarget, totl, 100, 0.0, 12.0)
 
 
ut.bookHist(h, "origin", "r vs z", 100, ztarget, totl, 100, 0.0, 15.0)
ut.bookHist(h, "borigin", "r vs z", 100, ztarget, totl, 500, 0.0, 30.0)
ut.bookHist(h, "porigin", "r vs z secondary", 100, ztarget, totl, 500, 0.0, 30.0)
ut.bookHist(h, "mu-inter", "r vs z", 100, ztarget, totl, 50, 0.0, 3.0)
ut.bookHist(h, "muonP", "muon momentum", 100, 0.0, 400.0)
ut.bookHist(h, "weight", "muon weight", 1000, 1.0e3, 1.0e6)
ut.bookHist(h, "delx13", "x diff first and last point,  mu-", 100, -10.0, 10.0)
ut.bookHist(h, "delx-13", "x diff first and last point, mu+", 100, -10.0, 10.0)
ut.bookHist(h, "deltaElec", "energy of delta electrons", 100, 0.0, 10.0)
ut.bookHist(h, "secondaries", "energy of delta electrons", 100, 0.0, 10.0)
ut.bookHist(h, "prop_deltaElec", "energy of delta elec / muon energy", 100, 0.0, 1.0)
ut.bookHist(h, "dummy", "", 10, 0.0, 1.0)
ut.bookHist(h, "dE", "", 100, 0.0, 10.0, 100, 0.0, 10.0)
h["dummy"].SetMaximum(10.0)
 
histlistAll = {
    1: "strawstation_5",
    2: "strawstation_1",
    3: "strawstation_4",
    4: "Ecal",
    5: "muondet",
    6: "VetoTimeDet",
    7: "T1LiSc",
    8: "T2LiSc",
    9: "T3LiSc",
    10: "T5LiSc",
    11: "ShipRpc",
    12: "volHPT",
    13: "Target",
    14: "TimeDet",
    15: "Det2",
}
hLiSc = {1: {}}
for i in range(1, 7):
    hLiSc[1][i] = "T1LiSc_" + str(i)
hLiSc[2] = {}
for i in range(1, 48):
    hLiSc[2][i] = "T2LiSc_" + str(i)
hLiSc[3] = {}
for i in range(1, 3):
    hLiSc[3][i] = "T3LiSc_" + str(i)
hLiSc[5] = {}
for i in range(1, 3):
    hLiSc[5][i] = "T5LiSc_" + str(i)
hMuon = {}
for i in range(0, 4):
    hMuon[i] = "muondet" + str(i)
 
 
# start event loop
mom = ROOT.TVector3()
pos = ROOT.TVector3()
 
 
def BigEventLoop() -> None:
    pid = 1
    for fn in fchain:
        if os.path.islink(fn):
            rfn = os.path.realpath(fn).split("eos")[1]
            fn = ROOT.gSystem.Getenv("EOSSHIP") + "/eos/" + rfn
        elif not os.path.isfile(fn):
            raise RuntimeError("Don't know what to do with " + fn)
        if parallel:
            # process files parallel instead of sequential
            processes.append(mp.Process(target=executeOneFile, args=(fn, output, pid)))
            pid += 1
        else:
            processes.append(fn)
    # Run processes
    n = 0
    for p in processes:
        if parallel:
            p.start()
            n += 1
        else:
            executeOneFile(p)
    if parallel:
        # Exit the completed processes
        for p in processes:
            p.join()
            # clean histos before reading in the new ones
        for x in h:
            h[x].Reset()
        print("now, collect the output")
        pid = 1
        for p in processes:
            ut.readHists(h, "tmpHists_" + str(pid) + ".root")
            pid += 1
    # compactify liquid scintillator
    for mu in ["", "_mu", "_muV0"]:
        for x in ["", "_E", "_P", "_LP", "_OP", "_id", "_mul", "_evmul", "_origin", "_originmu"]:
            for k in [1, 2, 3, 5]:
                first = True
                for j in hLiSc[k]:
                    detName = hLiSc[k][j]
                    tag = detName + mu + x
                    newh = detName[0:2] + "LiSc" + mu + x
                    if tag not in h:
                        continue
                    if first:
                        h[newh] = h[tag].Clone(newh)
                        h[newh].SetTitle(h[tag].GetTitle().split("_")[0])
                        first = False
                    else:
                        h[newh].Add(h[tag])
    # compactify muon stations
    for mu in ["", "_mu", "_muV0"]:
        for x in ["", "_E", "_P", "_LP", "_OP", "_id", "_mul", "_evmul", "_origin", "_originmu"]:
            first = True
            for j in hMuon:
                detName = hMuon[j]
                tag = detName + mu + x
                newh = "muondet" + mu + x
                if first:
                    h[newh] = h[tag].Clone(newh)
                    h[newh].SetTitle(h[tag].GetTitle().split(" ")[0] + " " + newh)
                    first = False
                else:
                    h[newh].Add(h[tag])
 
    # make list of hists with entries
    k = 1
    for x in histlistAll:
        if histlistAll[x] in h:
            histlist[k] = histlistAll[x]
            # make cumulative histograms
            for c in ["", "_E", "_P", "_LP", "_OP", "_id", "_mul", "_evmul", "_origin", "_originmu"]:
                h[histlist[k] + "_mu" + c].Add(h[histlist[k] + "_muV0" + c])
                h[histlist[k] + c].Add(h[histlist[k] + "_mu" + c])
                h[histlist[k] + c].SetMinimum(0.0)
                h[histlist[k] + "_mu" + c].SetMinimum(0.0)
                h[histlist[k] + "_muV0" + c].SetMinimum(0.0)
            k += 1
    nstations = len(histlist)
    makePlots(nstations)
 
 
def executeOneFile(fn, output=None, pid=None) -> None:
    f = ROOT.TFile.Open(fn)
    sTree = f.Get("cbmsim")
    nEvents = sTree.GetEntries()
    if sTree.GetBranch("GeoTracks"):
        sTree.SetBranchStatus("GeoTracks", 0)
    sTree.GetEntry(0)
    hitContainers = [
        sTree.vetoPoint,
        sTree.muonPoint,
        sTree.EcalPointLite,
        sTree.strawtubesPoint,
        sTree.ShipRpcPoint,
        sTree.TargetPoint,
    ]
    ntot = 0
    for n in range(nEvents):
        sTree.GetEntry(n)
        theMuon = sTree.MCTrack[0]
        if len(sTree.MCTrack) > 1:
            w = sTree.MCTrack[1].GetWeight()  # also works for neutrinos
        else:
            print("should not happen with new files", n, fn)
            w = sTree.MCTrack[0].GetWeight()  # also works for neutrinos
        if w == 0:
            w = 1.0
        h["weight"].Fill(w)
        h["muonP"].Fill(theMuon.GetP() / u.GeV, w)
        ntot += 1
        if ntot % 10000 == 0:
            print("read event", f.GetName(), n)
        hitmult = {}
        esum = 0
        for mcp in sTree.MCTrack:
            if mcp.GetMotherId() == 0:  # mother is original muon
                E = mcp.GetEnergy()
                if mcp.GetPdgCode() == 11:
                    h["deltaElec"].Fill(E, w)
                    esum += E
                h["secondaries"].Fill(E, w)
        h["prop_deltaElec"].Fill(
            esum / sTree.MCTrack[0].GetP(), w
        )  # until energy is really made persistent GetEnergy()
        for c in hitContainers:
            for ahit in c:
                if not ahit.GetEnergyLoss() > 0:
                    continue
                detID = ahit.GetDetectorID()
                if ahit.GetName() == "strawtubesPoint":
                    detName = "strawstation_" + str(ahit.GetStationNumber())
                    x = ahit.GetX()
                    y = ahit.GetY()
                    E = ahit.GetEnergyLoss()
                else:
                    if detID not in logVols:
                        detName = c.GetName().replace("Points", "")
                        if detName not in histlistAll.values():
                            print(detID, detName, c.GetName())
                    else:
                        detName = logVols[detID]
                    x = ahit.GetX()
                    y = ahit.GetY()
                    ahit.GetZ()
                    E = ahit.GetEnergyLoss()
                if detName not in h:
                    bookHist(detName)
                mu = ""
                pdgID = -2
                if "PdgCode" in dir(ahit):
                    pdgID = ahit.PdgCode()
                trackID = ahit.GetTrackID()
                phit = -100.0
                mom = ROOT.TVector3()
                if not trackID < 0:
                    aTrack = sTree.MCTrack[trackID]
                    pdgID = aTrack.GetPdgCode()
                    aTrack.GetMomentum(mom)  # ! this is not momentum of particle at Calorimeter place
                    phit = mom.Mag() / u.GeV
                if abs(pdgID) == 13:
                    mu = "_mu"
                if phit > 3 and abs(pdgID) == 13:
                    mu = "_muV0"
                detName = detName + mu
                if detName.find("LS") < 0:
                    h[detName].Fill(x / u.m, y / u.m, w)
                else:
                    h[detName].Fill(ahit.GetZ() / u.m, ROOT.TMath.ATan2(y, x) / ROOT.TMath.Pi(), w)
                h[detName + "_E"].Fill(E / u.MeV, w)
                if detName not in hitmult:
                    hitmult[detName] = {-1: 0}
                if trackID not in hitmult[detName]:
                    hitmult[detName][trackID] = 0
                hitmult[detName][trackID] += 1
                hitmult[detName][-1] += 1
                h[detName + "_id"].Fill(pdgID, w)
                h[detName + "_P"].Fill(phit, w)
                h[detName + "_LP"].Fill(phit, w)
                if not trackID < 0:
                    r = ROOT.TMath.Sqrt(aTrack.GetStartX() ** 2 + aTrack.GetStartY() ** 2) / u.m
                    h["origin"].Fill(aTrack.GetStartZ() / u.m, r, w)
                    h[detName + "_origin"].Fill(aTrack.GetStartZ() / u.m, r, w)
                    if abs(pdgID) == 13:
                        h[detName + "_originmu"].Fill(aTrack.GetStartZ() / u.m, r, w)
                    h["borigin"].Fill(aTrack.GetStartZ() / u.m, r, w)
                    aTrack.GetMomentum(mom)
                    h[detName + "_OP"].Fill(mom.Mag() / u.GeV, w)
                    if trackID > 0:
                        origin(sTree, trackID)
                        h["porigin"].Fill(
                            aTrack.GetStartZ() / u.m,
                            ROOT.TMath.Sqrt(aTrack.GetStartX() ** 2 + aTrack.GetStartY() ** 2) / u.m,
                            w,
                        )
                    h["mu-inter"].Fill(rz_inter[1] / u.m, rz_inter[0] / u.m, w)
        for detName in hitmult:
            h[detName + "_evmul"].Fill(hitmult[detName][-1], w)
            for tr in hitmult[detName]:
                h[detName + "_mul"].Fill(hitmult[detName][tr], w)
    if output:
        ut.writeHists(h, "tmpHists_" + str(pid) + ".root")
        output.put("ok")
    f.Close()
 
 
#
def makePlots(nstations: int) -> None:
    cxcy = {
        1: [2, 1],
        2: [3, 1],
        3: [2, 2],
        4: [3, 2],
        5: [3, 2],
        6: [3, 2],
        7: [3, 3],
        8: [3, 3],
        9: [3, 3],
        10: [4, 3],
        11: [4, 3],
        12: [4, 3],
        13: [5, 3],
        14: [5, 3],
        15: [5, 3],
    }
    ut.bookCanvas(
        h,
        key="ResultsI",
        title="hit occupancies",
        nx=1100,
        ny=600 * cxcy[nstations][1],
        cx=cxcy[nstations][0],
        cy=cxcy[nstations][1],
    )
    ut.bookCanvas(
        h,
        key="ResultsImu",
        title="hit occupancies",
        nx=1100,
        ny=600 * cxcy[nstations][1],
        cx=cxcy[nstations][0],
        cy=cxcy[nstations][1],
    )
    ut.bookCanvas(
        h,
        key="ResultsImuV0",
        title="hit occupancies",
        nx=1100,
        ny=600 * cxcy[nstations][1],
        cx=cxcy[nstations][0],
        cy=cxcy[nstations][1],
    )
    ut.bookCanvas(
        h,
        key="ResultsII",
        title="deposited energies",
        nx=1100,
        ny=600 * cxcy[nstations][1],
        cx=cxcy[nstations][0],
        cy=cxcy[nstations][1],
    )
    ut.bookCanvas(
        h,
        key="ResultsIII",
        title="track info",
        nx=1100,
        ny=600 * cxcy[nstations][1],
        cx=cxcy[nstations][0],
        cy=cxcy[nstations][1],
    )
    ut.bookCanvas(
        h,
        key="ResultsIV",
        title="track info",
        nx=1100,
        ny=600 * cxcy[nstations][1],
        cx=cxcy[nstations][0],
        cy=cxcy[nstations][1],
    )
    ut.bookCanvas(h, key="ResultsV", title="origin info", nx=600, ny=400, cx=1, cy=1)
    txt[0] = "occupancy 1sec 5E13pot "
    h["dummy"].Fill(-1)
    nSpills = len(prefixes)
    for i in histlist:
        tc = h["ResultsI"].cd(i)
        hn = histlist[i]
        if hn not in h:
            continue
        h[hn].SetStats(0)
        h[hn].Draw("colz")
        txt[i] = "%5.2G" % (h[hn].GetSumOfWeights() / nSpills)
        labels[i] = ROOT.TLatex().DrawLatexNDC(0.15, 0.85, txt[i])
        #
        hn = histlist[i] + "_mu"
        tc = h["ResultsImu"].cd(i)
        h[hn].SetStats(0)
        h[hn].Draw("colz")
        txt[i + 100] = "%5.2G" % (h[hn].GetSumOfWeights() / nSpills)
        labels[i + 100] = ROOT.TLatex().DrawLatexNDC(0.15, 0.85, txt[i + 100])
        #
        hn = histlist[i] + "_muV0"
        tc = h["ResultsImuV0"].cd(i)
        h[hn].SetStats(0)
        h[hn].Draw("colz")
        txt[i + 200] = "%5.2G" % (h[hn].GetSumOfWeights() / nSpills)
        labels[i + 200] = ROOT.TLatex().DrawLatexNDC(0.15, 0.85, txt[i + 200])
    #
    for i in histlist:
        tc = h["ResultsII"].cd(i)
        h[histlist[i] + "_E"].Draw("colz")
    #
    for i in histlist:
        tc = h["ResultsIII"].cd(i)
        tc.SetLogy(1)
        hn = histlist[i]
        drawBoth("_P", hn)
        hid = h[hn + "_id"]
        print("particle statistics for " + histlist[i])
        for k in range(hid.GetNbinsX()):
            ncont = hid.GetBinContent(k + 1)
            pid = hid.GetBinCenter(k + 1)
            if ncont > 0:
                temp = int(abs(pid) + 0.5) * int(pid / abs(pid))
                nm = PDG.GetParticle(temp).GetName()
                print(f"{nm} :{ncont:5.2g}")
        hid = h[histlist[i] + "_mu_id"]
        for k in range(hid.GetNbinsX()):
            ncont = hid.GetBinContent(k + 1)
            pid = hid.GetBinCenter(k + 1)
            if ncont > 0:
                temp = int(abs(pid) + 0.5) * int(pid / abs(pid))
                nm = PDG.GetParticle(temp).GetName()
                print(f"{nm} :{ncont:5.2g}")
        #
        tc = h["ResultsV"].cd(1)
        h["origin"].SetStats(0)
        h["origin"].Draw("colz")
    #
    for i in histlist:
        tc = h["ResultsIV"].cd(i)
        tc.SetLogy(1)
        hn = histlist[i]
        drawBoth("_OP", hn)
    persistency()
 
 
#
def AnaEventLoop() -> None:
    fout = open("rareEvents.txt", "w")  # noqa: SIM115
    for fn in fchainRec:
        f = ROOT.TFile(fn)
        if not f.FindObjectAny("cbmsim"):
            print("skip file ", f.GetName())
            continue
        sTree = f.Get("cbmsim")
        sTree.GetEvent(0)
        if sTree.GetBranch("GeoTracks"):
            sTree.SetBranchStatus("GeoTracks", 0)
        nEvents = sTree.GetEntries()
        for n in range(nEvents):
            sTree.GetEntry(n)
            if n == 0:
                print("now at event ", n, f.GetName())
            if len(sTree.MCTrack) > 1:
                wg = sTree.MCTrack[1].GetWeight()  # also works for neutrinos
            else:
                wg = sTree.MCTrack[0].GetWeight()  # also works for neutrinos
            i = -1
            for atrack in sTree.FitTracks:
                i += 1
                fitStatus = atrack.getFitStatus()
                if not fitStatus.isFitConverged():
                    continue
                nmeas = atrack.getNumPoints()
                chi2 = fitStatus.getChi2() / nmeas
                fittedState = atrack.getFittedState()
                P = fittedState.getMomMag()
                fout.write("rare event with track %i, %s, %8.2F \n" % (n, f.GetName(), wg))
                originOfMuon(fout, n, f.GetName(), nEvents)
                fout.write("reco %i,%5.3F,%8.2F \n" % (nmeas, chi2, P))
                mcPartKey = sTree.fitTrack2MC[i]
                mcPart = sTree.MCTrack[mcPartKey]
                for ahit in sTree.strawtubesPoint:
                    if ahit.GetTrackID() == mcPartKey:
                        fout.write(
                            "P when making hit %i, %8.2F \n"
                            % (
                                ahit.PdgCode(),
                                ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2 + ahit.GetPz() ** 2) / u.GeV,
                            )
                        )
                        break
                if not mcPart:
                    continue
                Ptruth = mcPart.GetP()
                fout.write("Ptruth %i %8.2F \n" % (mcPart.GetPdgCode(), Ptruth / u.GeV))
 
 
#
def muDISntuple(fn) -> None:
    # take as input the rare events
    fout = ROOT.TFile("muDISVetoCounter.root", "recreate")
    h["ntuple"] = ROOT.TNtuple("muons", "muon flux VetoCounter", "id:px:py:pz:x:y:z:w")
    f = ROOT.TFile(fn)
    sTree = f.Get("cbmsim")
    nEvents = sTree.GetEntries()
    for n in range(nEvents):
        sTree.GetEntry(n)
        if len(sTree.MCTrack) > 1:
            wg = sTree.MCTrack[1].GetWeight()
        else:
            wg = sTree.MCTrack[0].GetWeight()
        for ahit in sTree.vetoPoint:
            detID = ahit.GetDetectorID()
            if logVols[detID] != "VetoTimeDet":
                continue
            pid = ahit.PdgCode()
            if abs(pid) != 13:
                continue
            P = ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2 + ahit.GetPz() ** 2)
            if 3 / u.GeV < P:
                h["ntuple"].Fill(
                    float(pid),
                    float(ahit.GetPx() / u.GeV),
                    float(ahit.GetPy() / u.GeV),
                    float(ahit.GetPz() / u.GeV),
                    float(ahit.GetX() / u.m),
                    float(ahit.GetY() / u.m),
                    float(ahit.GetZ() / u.m),
                    float(wg),
                )
    fout.cd()
    h["ntuple"].Write()
 
 
def analyzeConcrete() -> None:
    h["fout"] = ROOT.TFile("muConcrete.root", "recreate")
    h["ntuple"] = ROOT.TNtuple("muons", "muon flux concrete", "id:px:py:pz:x:y:z:w")
    for m in ["", "mu", "V0"]:
        ut.bookHist(h, "conc_hitz" + m, "concrete hit z " + m, 100, -100.0, 100.0)
        ut.bookHist(h, "conc_hitzP" + m, "concrete hit z vs P" + m, 100, -100.0, 100.0, 100, 0.0, 25.0)
        ut.bookHist(h, "conc_hity" + m, "concrete hit y " + m, 100, -15.0, 15.0)
        ut.bookHist(h, "conc_p" + m, "concrete hit p " + m, 1000, 0.0, 400.0)
        ut.bookHist(h, "conc_pt" + m, "concrete hit pt " + m, 100, 0.0, 20.0)
        ut.bookHist(h, "conc_hitzy" + m, "concrete hit zy " + m, 100, -100.0, 100.0, 100, -15.0, 15.0)
    top = sGeo.GetTopVolume()
    magn = top.GetNode("magyoke_1")
    z0 = magn.GetMatrix().GetTranslation()[2] / u.m
    for fn in fchain:
        f = ROOT.TFile(fn)
        if not f.FindObjectAny("cbmsim"):
            print("skip file ", f.GetName())
            continue
        sTree = f.Get("cbmsim")
        nEvents = sTree.GetEntries()
        ROOT.gROOT.cd()
        for n in range(nEvents):
            sTree.GetEntry(n)
            if len(sTree.MCTrack) > 1:
                wg = sTree.MCTrack[1].GetWeight()
            else:
                wg = sTree.MCTrack[0].GetWeight()
            for ahit in sTree.vetoPoint:
                detID = ahit.GetDetectorID()
                if logVols[detID] != "rockD":
                    continue
                m = ""
                pid = ahit.PdgCode()
                if abs(pid) == 13:
                    m = "mu"
                P = ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2 + ahit.GetPz() ** 2)
                if abs(pid) == 13 and 3 / u.GeV < P:
                    m = "V0"
                    h["ntuple"].Fill(
                        float(pid),
                        float(ahit.GetPx() / u.GeV),
                        float(ahit.GetPy() / u.GeV),
                        float(ahit.GetPz() / u.GeV),
                        float(ahit.GetX() / u.m),
                        float(ahit.GetY() / u.m),
                        float(ahit.GetZ() / u.m),
                        float(wg),
                    )
                h["conc_hitz" + m].Fill(ahit.GetZ() / u.m - z0, wg)
                h["conc_hity" + m].Fill(ahit.GetY() / u.m, wg)
                h["conc_p" + m].Fill(P / u.GeV, wg)
                h["conc_hitzP" + m].Fill(ahit.GetZ() / u.m, P / u.GeV, wg)
                Pt = ROOT.TMath.Sqrt(ahit.GetPx() ** 2 + ahit.GetPy() ** 2)
                h["conc_pt" + m].Fill(Pt / u.GeV, wg)
                h["conc_hitzy" + m].Fill(ahit.GetZ() / u.m - z0, ahit.GetY() / u.m, wg)
    #
    # start = [ahit.GetX()/u.m,ahit.GetY()/u.m,ahit.GetZ()/u.m]
    # direc = [-ahit.GetPx()/P,-ahit.GetPy()/P,-ahit.GetPz()/P]
    # t = - start[0]/direc[0]
 
    ut.bookCanvas(h, key="ResultsV0", title="muons hitting concrete, p>3GeV", nx=1000, ny=600, cx=2, cy=2)
    ut.bookCanvas(h, key="Resultsmu", title="muons hitting concrete", nx=1000, ny=600, cx=2, cy=2)
    ut.bookCanvas(h, key="Results", title="hitting concrete", nx=1000, ny=600, cx=2, cy=2)
    for m in ["", "mu", "V0"]:
        tc = h["Results" + m].cd(1)
        h["conc_hity" + m].Draw()
        tc = h["Results" + m].cd(2)
        h["conc_hitz" + m].Draw()
        tc = h["Results" + m].cd(3)
        tc.SetLogy(1)
        h["conc_pt" + m].Draw()
        tc = h["Results" + m].cd(4)
        tc.SetLogy(1)
        h["conc_p" + m].Draw()
        h["fout"].cd()
        h["ntuple"].Write()
 
 
def rareEventEmulsion(fname: str = "rareEmulsion.txt") -> None:
    fout = open(fname, "w")  # noqa: SIM115
    for fn in fchainRec:
        f = ROOT.TFile(fn)
        if not f.FindObjectAny("cbmsim"):
            print("skip file ", f.GetName())
            continue
        sTree = f.Get("cbmsim")
        sTree.GetEvent(0)
        if sTree.GetBranch("GeoTracks"):
            sTree.SetBranchStatus("GeoTracks", 0)
        nEvents = sTree.GetEntries()
        for n in range(nEvents):
            sTree.GetEntry(n)
            if n == 0:
                print("now at event ", n, f.GetName())
            for ahit in sTree.vetoPoint:
                detID = ahit.GetDetectorID()
                if logVols[detID] != "Emulsion":
                    continue
                ahit.GetX()
                ahit.GetY()
                ahit.GetZ()
                if len(sTree.MCTrack) > 1:
                    wg = sTree.MCTrack[1].GetWeight()  # also works for neutrinos
                else:
                    wg = sTree.MCTrack[0].GetWeight()  # also works for neutrinos
                fout.write("rare emulsion hit %i, %s, %8.3F, %i \n" % (n, f.GetName(), wg, ahit.PdgCode()))
                if ahit.GetPz() / u.GeV > 1.0:
                    fout.write(
                        "V,P when making hit %8.3F,%8.3F,%8.3F %8.3F,%8.3F,%8.3F (GeV) \n"
                        % (
                            ahit.GetX() / u.m,
                            ahit.GetY() / u.m,
                            ahit.GetZ() / u.m,
                            ahit.GetPx() / u.GeV,
                            ahit.GetPy() / u.GeV,
                            ahit.GetPz() / u.GeV,
                        )
                    )
                else:
                    fout.write(
                        "V,P when making hit %8.3F,%8.3F,%8.3F %8.3F,%8.3F,%8.3F (MeV)\n"
                        % (
                            ahit.GetX() / u.m,
                            ahit.GetY() / u.m,
                            ahit.GetZ() / u.m,
                            ahit.GetPx() / u.MeV,
                            ahit.GetPy() / u.MeV,
                            ahit.GetPz() / u.MeV,
                        )
                    )
                originOfMuon(fout, n, f.GetName(), nEvents)
 
 
#
def extractRareEvents(single=None) -> None:
    for fn in fchainRec:
        if single:
            if fn.find(str(single)) < 0:
                continue
        f = ROOT.TFile(fn)
        if not f.FindObjectAny("cbmsim"):
            print("skip file ", f.GetName())
            continue
        sTree = f.Get("cbmsim")
        nEvents = sTree.GetEntries()
        raref = ROOT.TFile(fn.replace(".root", "_rare.root"), "recreate")
        newTree = sTree.CloneTree(0)
        for n in range(nEvents):
            sTree.GetEntry(n)
            if n == 0:
                print("now at event ", n, f.GetName())
                # count fitted tracks which have converged and nDF>20:
            n = 0
            for fT in sTree.FitTracks:
                fst = fT.getFitStatus()
                if not fst.isFitConverged():
                    continue
                if fst.getNdf() < 20:
                    continue
                n += 1
            if n > 0:
                rc = newTree.Fill()
                print("filled newTree", rc)
            sTree.Clear()
        newTree.AutoSave()
        print("   --> events saved:", newTree.GetEntries())
        f.Close()
        raref.Close()
 
 
#
def extractMuCloseByEvents(single=None) -> None:
    mom = ROOT.TVector3()
    pos = ROOT.TVector3()
    # Goliath magnet has been removed from geometry
    goliath_node = top.GetNode("volGoliath_1")
    if goliath_node:
        golmx = goliath_node.GetMatrix()
        zGol = golmx.GetTranslation()[2]
    else:
        zGol = -9999.0  # No Goliath, accept all z positions
    for fn in fchainRec:
        if single:
            if fn.find(str(single)) < 0:
                continue
        f = ROOT.TFile(fn)
        if not f.FindObjectAny("cbmsim"):
            print("skip file ", f.GetName())
            continue
        sTree = f.Get("cbmsim")
        nEvents = sTree.GetEntries()
        raref = ROOT.TFile(fn.replace(".root", "_clby.root"), "recreate")
        newTree = sTree.CloneTree(0)
        for n in range(nEvents):
            sTree.GetEntry(n)
            if n == 0:
                print("now at event ", n, f.GetName())
                # look for muons p>3 hitting something
            n = 0
            for c in [sTree.vetoPoint, sTree.strawtubesPoint, sTree.ShipRpcPoint]:
                for ahit in c:
                    if abs(ahit.PdgCode()) != 13:
                        continue
                    ahit.Momentum(mom)
                    if mom.Mag() < 3.0:
                        continue
                    ahit.Position(pos)
                    if pos.z() < zGol:
                        continue
                    n += 1
            if n > 0:
                newTree.Fill()
                # print 'filled newTree',rc
            sTree.Clear()
        newTree.AutoSave()
        print("   --> events saved:", newTree.GetEntries())
        f.Close()
        raref.Close()
 
 
#
def MergeRareEvents(runs: list[str] | None = None) -> None:
    if runs is None:
        runs = ["61", "62"]
    for prefix in runs:
        cmd = "$ROOTSYS/bin/hadd rareEvents_" + prefix + ".root -f "
        for fn in fchainRec:
            if fn.find("muon" + prefix) < 0:
                continue
            fr = fn.replace(".root", "_rare.root")
            cmd = cmd + " " + fr
        os.system(cmd)
 
 
#
def persistency() -> None:
    printAndCopy(prefix)
    ut.writeHists(h, prefix + ".root", plusCanvas=True)
 
 
def reDraw(fn) -> None:
    if fn.find("root") < 0:
        fn = fn + ".root"
    if "tc" not in h:
        h["tc"] = ROOT.TFile(fn)
    for x in ["ResultsI", "ResultsII", "ResultsImu", "ResultsImuV0", "ResultsIII", "ResultsIV", "ResultsV"]:
        h[x] = h["tc"].FindObjectAny(x)
        h[x].Draw()
 
 
def printAndCopy(prefix: str | None = None) -> None:
    if not prefix:
        prefix = (h["tc"].GetName()).replace(".root", "")
    for x in ["ResultsI", "ResultsII", "ResultsImu", "ResultsImuV0", "ResultsIII", "ResultsIV", "ResultsV"]:
        h[x].Update()
    if prefix not in os.listdir("."):
        os.mkdir(prefix)
    os.chdir(prefix)
    h["ResultsI"].Print(prefix + "Back_occ.png")
    h["ResultsII"].Print(prefix + "Back_depE.png")
    h["ResultsImu"].Print(prefix + "muBack_occ.png")
    h["ResultsImuV0"].Print(prefix + "muV0Back_occ.png")
    h["ResultsIII"].Print(prefix + "Back_P.png")
    h["ResultsIV"].Print(prefix + "Back_OP.png")
    h["ResultsV"].Print(prefix + "origin.png")
    os.chdir("../")
 
 
def drawBoth(tag: str, hn) -> None:
    n1 = h[hn + "_mu" + tag].GetMaximum()
    n2 = h[hn + tag].GetMaximum()
    if n1 > n2:
        h[hn + tag].SetMaximum(n1)
    h[hn + "_mu" + tag].SetLineColor(4)
    h[hn + tag].SetLineColor(3)
    h[hn + "_mu" + tag].Draw()
    h[hn + tag].Draw("same")
 
 
def debugGeoTracks(sTree) -> None:
    for i in range(sTree.GetEntries()):
        sTree.GetEntry(i)
        n = 0
        for gt in sTree.GeoTracks:
            print(
                i, n, gt.GetFirstPoint()[2], gt.GetLastPoint()[2], gt.GetParticle().GetPdgCode(), gt.GetParticle().P()
            )
            n += 1
 
 
def eventsWithStrawPoints(i) -> None:
    sTree = fchain[i].Get("cbmsim")
    mom = ROOT.TVector3()
    for i in range(sTree.GetEntries()):
        sTree.GetEntry(i)
        nS = len(sTree.strawtubesPoint)
        if nS > 0:
            mu = sTree.MCTrack[0]
            mu.GetMomentum(mom)
            print(i, nS)
            mom.Print()
            sp = sTree.strawtubesPoint[(max(0, nS - 3))]
            sp.Momentum(mom)
            mom.Print()
            print("-----------------------")
 
 
def eventsWithEntryPoints(i) -> None:
    sTree = fchain[i].Get("cbmsim")
    mom = ROOT.TVector3()
    for i in range(sTree.GetEntries()):
        sTree.GetEntry(i)
        np = 0
        for vp in sTree.vetoPoint:
            detName = logVols[vp.GetDetectorID()]
            if detName == "VetoTimeDet":
                np += 0  # ??
            vp.Momentum(mom)
            print(i, detName, vp.PdgCode())
            mom.Print()
            print("-----------------------")
 
 
def depEnergy(sTree) -> None:
    for n in range(sTree.GetEntries()):
        sTree.GetEntry(n)
        for ahit in sTree.strawtubesPoint:
            dE = ahit.GetEnergyLoss() / u.keV
            ahit.Momentum(mom)
            pa = PDG.GetParticle(ahit.PdgCode())
            mpa = pa.Mass()
            E = ROOT.TMath.Sqrt(mom.Mag2() + mpa**2)
            ekin = E - mpa
            h["dE"].Fill(dE, ekin / u.MeV)
        h["dE"].SetXTitle("keV")
        h["dE"].SetYTitle("MeV")
 
 
def originOfMuon(fout: TextIOWrapper, n: int, fn, nEvents) -> None:
    # from fn extract Yandex or CERN/cracow
    ncpu = 9
    x = fn.find("/")
    ni = int(fn[x - 1 : x]) - 1
    if nEvents < 100000:
        fm = "$SHIPSOFT/data/pythia8_Geant4_onlyMuons.root"
    else:
        fm = "$SHIPSOFT/data/pythia8_Geant4_Yandex_onlyMuons.root"
    fmuon = ROOT.TFile(fm)
    ntupl = fmuon.Get("pythia8-Geant4")
    ntot = ntupl.GetEntries()
    n3 = int(ntot / ncpu)
    N = n3 * ni + n
    ntupl.GetEvent(N)
    P = ROOT.TMath.Sqrt(ntupl.pz * ntupl.pz + ntupl.py * ntupl.py + ntupl.px * ntupl.px)
    fout.write("original muon %i, %i, %8.2F \n" % (ntupl.parentid, ntupl.pythiaid, P))
    fmuon.Close()
 
 
#
# BigEventLoop()
# makePlots()
# AnaEventLoop()
 
 
# ShipAna
def pers() -> None:
    xdisk = "/media/Work/HNL/"
    for x in h:
        if isinstance(h[x], ROOT.TCanvas):
            h[x].Update()
            tn = h[x].GetName() + ".png"
            h[x].Print(tn)
            rpath = os.path.abspath(".").split("/HNL/")[1]
            lp = rpath.split("/")
            prefix = xdisk
            for i in range(len(lp)):
                if lp[i] not in os.listdir(prefix):
                    os.system("mkdir " + prefix + "/" + lp[i])
                prefix = prefix + "/" + lp[i]
            os.system("cp " + tn + " " + xdisk + rpath)
 
 
from operator import itemgetter
 
 
def makeNicePrintout(x: list[str] | None = None):
    if x is None:
        x = ["rareEvents_61-62.txt", "rareEvents_71-72.txt"]
    result = []
    cor = 1.0
    for fn in x:
        with open(fn) as f:
            recTrack = None
            if fn == "rareEvents_81-102.txt":
                cor = 30.0
            for lx in f.readlines():
                line = lx.replace("\n", "")
                if not line.find("rare event") < 0:
                    if recTrack:
                        result.append(recTrack)
                    tmp = line.split(",")
                    w = tmp[2].replace(" ", "")
                    ff = tmp[1].split("/")[0].replace(" ", "")
                    recTrack = {"w": w, "file": ff}
                elif not line.find("original") < 0:
                    tmp = line.split(",")
                    recTrack["origin"] = tmp[0].split(" ")[2]
                    recTrack["pytiaid"] = tmp[1].replace(" ", "")
                    recTrack["o-mom"] = tmp[2].replace(" ", "")
                elif not line.find("reco ") < 0:
                    tmp = line.split(",")
                    recTrack["nmeas"] = tmp[0].split(" ")[1]
                    recTrack["chi2"] = tmp[1]
                    recTrack["p_rec"] = tmp[2].replace(" ", "")
                elif not line.find("making") < 0:
                    tmp = line.split(",")
                    recTrack["p_hit"] = tmp[1].replace(" ", "")
                    recTrack["fp_hit"] = float(tmp[1].replace(" ", ""))
                elif not line.find("Ptruth") < 0:
                    tmp = line.split(" ")
                    recTrack["id_hit"] = tmp[1].replace(" ", "")
    # print a table
    print(
        "%4s %8s %8s %4s %8s %8s %8s %8s %8s  %8s "
        % ("nr", "origin", "pythiaID", "ID", "p_orig", "p_hit", "chi2", "weight", "file", "cor w")
    )
    # sort according to p_hit
    tmp = sorted(result, key=itemgetter("fp_hit"))
    muonrate1 = 0
    muonrate2 = 0
    muonrate3 = 0
    for i in range(len(tmp)):
        tr = tmp[i]
        corw = float(tr["w"]) / cor
        if float(tr["p_hit"]) > 1:
            muonrate1 += corw
        if float(tr["p_hit"]) > 2:
            muonrate2 += corw
        if float(tr["p_hit"]) > 3:
            muonrate3 += corw
        print(
            "%4i %8s %8s %4s %8s %8s %8s %8s %8s  %8s "
            % (
                i,
                tr["origin"],
                tr["pytiaid"],
                tr["id_hit"],
                tr["o-mom"],
                tr["p_hit"],
                tr["chi2"],
                tr["w"],
                tr["file"],
                corw,
            )
        )
    print("guestimate for muonrate above 1GeV = ", muonrate1)
    print("guestimate for muonrate above 2GeV = ", muonrate2)
    print("guestimate for muonrate above 3GeV = ", muonrate3)
    # guestimate for muonrate above 1GeV =  56025.4793333
    # guestimate for muonrate above 2GeV =  25584.9546667
    # guestimate for muonrate above 3GeV =  14792.8113333
    return tmp
 
 
#
def readAndMergeHistos(prods) -> None:
    for p in prods:
        x = p
        if p.find(".root") < 0:
            x = p + ".root"
        ut.readHists(h, x)
    # make list of hists with entries
    k = 1
    for x in histlistAll:
        if histlistAll[x] in h:
            histlist[k] = histlistAll[x]
            k += 1
    nstations = len(histlist)
    print("make plots for ", nstations)
    makePlots(nstations)
    printAndCopy(prods[0].replace(".root", ""))
 
 
# 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
# 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
# python -i $HNL/ana_ShipMuon.py 1012 1013 1014 1015 1016 1017 1018 1019 1022 1023 1024 1025 1026 1027 1028 1029
# python -i $HNL/ana_ShipMuon.py 634 635 636 637 638 639 640 641 642 643
# combine all cern-cracow
# 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
# 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
# make muonDIS ntuple: muDISntuple("/media/Data/HNL/muonBackground/rareEvents_81-102.root") -> 'muDISVetoCounter.root'
#                      second step python $FAIRSHIP/muonShieldOptimization/makeMuonDIS.py 1 10000 muDISVetoCounter.root
#                      third step run_simScript.py --MuDIS -n 10 -f  muonDis_1.root
# for concrete
# analyzeConcrete() -> muConcrete.root