runCharmHadProd.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
 
import os
import sys
 
import ROOT
 
ncpus = 20
msel = "4"
# msel = "5"
 
if msel == "4":
    nev = 2000000  # 0.1s / event
    path = "/afs/cern.ch/project/lbcern/vol1/truf/charm/"
    chicc = 1.7e-3  # prob to produce primary ccbar pair/pot
else:
    nev = 1000000
    path = "/afs/cern.ch/project/lbcern/vol1/truf/beauty/"
    chicc = 1.6e-7  # prob to produce primary bbbar pair/pot
 
# path = "/home/truf/charm/"
 
# store only charm
 
# story about weights, run_fixedTarget puts into file header: pot = nrpotspill / wspill
# with wspill = nrpotspill*chicc/nrcpot*nEvents/nev -> pot =  nev * nrcpot / (chicc * nEvents)
# nEvents nr of events in input file, number of events requested for job, ratio should be 1
# pot =  nrcpot / chicc
# nrcpot=((TH1F*)fin->Get("2"))->GetBinContent(1)/2.;  // pot are counted double, i.e. for each signal, i.e. pot/2.
 
 
def makeHadrons(run) -> None:
    for n in range(ncpus):
        x.Rndm() * 1000000000.0
        os.system("mkdir run" + str(run))
        os.chdir("run" + str(run))
        cmd = (
            "python $FAIRSHIP/macro/makeCascade.py -m "
            + msel
            + " -n "
            + str(nev)
            + " -t  Cascade-run"
            + str(run)
            + "-parp16-MSTP82-1-MSEL"
            + msel
            + ".root"
        )
        # if not run in runList:
        os.system(cmd + " >log" + str(run) + " &")
        os.chdir("../")
        run += 1
 
 
def makeBackground(run, cycle: int = 0) -> None:
    for n in range(ncpus):
        orun = run + cycle * 1000
        os.chdir("run" + str(run))
        inputFile = path + "run" + str(run) + "/Cascade-run" + str(run) + "-parp16-MSTP82-1-MSEL" + msel + ".root"
        f = ROOT.TFile(inputFile)
        nt = f.pythia6
        N = nt.GetEntries()
        f.Close()
        if msel == "4":
            cmd = (
                "python $FAIRSHIP/muonShieldOptimization/run_fixedTarget.py --force --charm -V -e 10 -P -n "
                + str(N)
                + " -r "
                + str(orun)
                + " -b 100 -X 100 -I "
                + inputFile
            )
        else:
            cmd = (
                "python $FAIRSHIP/muonShieldOptimization/run_fixedTarget.py --force --beauty -V -e 10 -P -n "
                + str(N)
                + " -r "
                + str(orun)
                + " -b 100 -X 100 -I "
                + inputFile
            )
        os.system(cmd + " >logFT" + str(orun) + " &")
        os.chdir("../")
        run += 1
 
 
cycle = 2
runList = [20, 28, 30, 33]
 
 
def makeBackgroundX(runList, cycle: int = 0) -> None:
    for run in runList:
        orun = run + cycle * 1000
        os.chdir("run" + str(run))
        inputFile = path + "run" + str(run) + "/Cascade-run" + str(run) + "-parp16-MSTP82-1-MSEL" + msel + ".root"
        f = ROOT.TFile(inputFile)
        nt = f.pythia6
        N = nt.GetEntries()
        f.Close()
        cmd = (
            "python $FAIRSHIP/muonShieldOptimization/run_fixedTarget.py --force --charm -V -e 10 -P -n "
            + str(N)
            + " -r "
            + str(orun)
            + " -b 100 -X 100 -I "
            + inputFile
        )
        if run in runList:
            os.system(cmd + " >logFT" + str(orun) + " &")
        os.chdir("../")
 
 
def merge(run, cycle: int = 0) -> None:
    fname = "pythia8_Geant4_XX_10.0.root"
    cmd = " "
    for n in range(ncpus):
        for x in os.listdir(path + "/run" + str(run + n)):
            orun = run + cycle * 1000
            if not x.find("run_fixedTarget_" + str(orun + n)) < 0:
                if cycle == 0 and run == 0 and not x.find("1001") < 0:
                    continue
                if cycle == 0 and run == 0 and not x.find("1010") < 0:
                    continue
                cmd += path + "/run" + str(run + n) + "/" + x + "/" + fname.replace("XX", str(orun + n)) + " "
    if msel == "4":
        outFile = fname.replace("XX", "charm_" + str(orun) + "-" + str(orun + ncpus - 1))
    else:
        outFile = fname.replace("XX", "beauty_" + str(orun) + "-" + str(orun + ncpus - 1))
    rc = os.system("hadd -O " + outFile + " " + cmd)
    if rc != 0:
        print("hadd failed, stop", outFile)
    else:
        rc = os.system("xrdcp " + outFile + " $EOSSHIP/eos/experiment/ship/data/Mbias/background-prod-2018/" + outFile)
        if rc != 0:
            print("copy to EOS failed, stop", outFile)
        else:
            rc = os.system("rm " + outFile)
 
 
def mergeAll() -> None:
    cmd = "hadd pythia8_Geant4_charm_153.3B_10.0_mu.root "
    tmp = "/eos/experiment/ship/data/Mbias/background-prod-2018/pythia8_Geant4_charm_XX_10.0_mu.root"
    for x in [
        "0-19",
        "20-39",
        "40-59",
        "60-79",
        "80-99",
        "1000-1019",
        "1020-1039",
        "1040-1059",
        "1060-1079",
        "1080-1099",
        "2000-2019",
        "2020-2039",
        "2040-2059",
        "2060-2079",
        "2080-2099",
    ]:
        cmd += tmp.replace("XX", x) + " "
    os.system(cmd)
 
 
def compactifyCascade(run) -> None:
    ncpus = 20
    cmd = ""
    Ntot = 0
    NperJob = nev
    for i in range(run, +ncpus):
        fName = path + "run" + str(i) + "/Cascade-run" + str(i) + "-parp16-MSTP82-1-MSEL" + msel + ".root"
        with open(path + "run" + str(i) + "/log" + str(i)) as f:
            success = False
            for line in f.readlines():
                if not line.find("Macro finished successfully") < 0:
                    success = True
        if not success:
            print("job not finished properly", fName)
            continue
        cmd += fName + " "
        Ntot += NperJob
    if cmd.find("root") < 0:
        print("no file found, exit")
    else:
        stat = str(int(Ntot / 1e6)) + "Mpot"
        outFile = (
            "Cascade-run"
            + str(run)
            + "-"
            + str(run + ncpus - 1)
            + "-parp16-MSTP82-1-MSEL"
            + msel
            + "-"
            + stat
            + ".root"
        )
        rc = os.system("hadd -O " + outFile + " " + cmd)
        f = ROOT.TFile(outFile)
        Npot = f.Get("2").GetBinContent(1) / 2.0 / chicc
        f.Close()
        stat = str(int(Npot / 1e9)) + "Bpot"
        oldOutFile = outFile
        outFile = (
            "Cascade-run"
            + str(run)
            + "-"
            + str(run + ncpus - 1)
            + "-parp16-MSTP82-1-MSEL"
            + msel
            + "-"
            + stat
            + ".root"
        )
        os.system("mv " + oldOutFile + " " + outFile)
        rc = os.system("xrdcp " + outFile + " $EOSSHIP/eos/experiment/ship/data/Mbias/background-prod-2018/" + outFile)
        if rc != 0:
            print("copy to EOS failed, stop", outFile)
        else:
            rc = os.system("rm " + outFile)
 
 
def statistics() -> None:
    path = os.environ["EOSSHIP"] + "/eos/experiment/ship/data/Mbias/background-prod-2018/"
    fname = "Cascade-runAA-BB-parp16-MSTP82-1-MSEL4-40Mpot.root"
    nPot = 0
    nhadrons = 0
    for x in [0, 20, 40, 60, 80]:
        fn = fname.replace("AA", str(x)).replace("BB", str(x + 19))
        f = ROOT.TFile.Open(path + fn)
        nPot += f.Get("2").GetBinContent(1) / 2.0
        nhadrons += f.Get("pythia6").GetEntries()
    print("total nr of hadrons:", nhadrons, nPot / chicc / 1.0e9, "Billion")
 
 
def potFromFileHeader(f) -> None:
    pot = 0
    for x in f.GetListOfKeys():
        if x.GetName() == "FileHeader":
            pot += float(x.GetTitle().split(" ")[3])
    print("PoT = ", pot)
 
 
x = ROOT.TRandom3()
x.SetSeed(0)
 
run = int(sys.argv[1])
 
print("following functions exist")
print(" - makeHadrons(run): will run makeCascade")
print(" - makeBackground(run): will run fixedTarget generator")