study_thinTarget.py

Go to the documentation of this file.

#!/usr/bin/env python
# SPDX-License-Identifier: LGPL-3.0-or-later
# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
 
import time
 
import ROOT
import shipRoot_conf
 
ROOT.gROOT.ProcessLine('#include "FairModule.h"')
time.sleep(20)
 
import geometry_config
import shipunit as u
 
mcEngine = "TGeant4"
runnr = 1
nev = 1000000
 
setup = {}
setup["TLV"] = {
    "thickness": 0.1 * u.cm,
    "material": "tungsten",
    "min momentum": 400 * u.GeV,
    "max momentum": 400 * u.GeV,
}
 
s = "TLV"
thickness = setup[s]["thickness"]
material = setup[s]["material"]
minmomentum = setup[s]["min momentum"]
maxmomentum = setup[s]["max momentum"]
 
checkOverlap = True
 
outFile = "TLV.root"
theSeed = 0
ecut = 0.0
 
# -------------------------------------------------------------------
# PYTHIA8 requires Random:seed to be in range [0, 900000000]
# When theSeed=0, ROOT generates a seed from system time which can exceed this limit
if theSeed == 0:
    ROOT.gRandom.SetSeed(0)  # Generate time-based seed
    theSeed = ROOT.gRandom.GetSeed()
    # Clamp to PYTHIA8's maximum allowed seed value
if theSeed > 900000000:
    theSeed = theSeed % 900000000
ROOT.gRandom.SetSeed(theSeed)
shipRoot_conf.configure()  # load basic libraries, prepare atexit for python
ship_geo = geometry_config.create_config(Yheight=10, shieldName="warm_opt")
 
# -----Timer--------------------------------------------------------
timer = ROOT.TStopwatch()
timer.Start()
 
# -----Create simulation run----------------------------------------
gFairBaseContFact = ROOT.FairBaseContFact()  # required by change to FairBaseContFact to avoid TList::Clear errors
run = ROOT.FairRunSim()
run.SetName(mcEngine)  # Transport engine
run.SetSink(ROOT.FairRootFileSink(outFile))  # Output file
run.SetUserConfig("g4Config.C")  # user configuration file default g4Config.C
rtdb = run.GetRuntimeDb()
 
# -----Materials----------------------------------------------
run.SetMaterials("media.geo")
# -----Create geometry----------------------------------------------
cave = ROOT.ShipCave("CAVE")
cave.SetGeometryFileName("cave.geo")
run.AddModule(cave)
 
 
class Block(ROOT.pyFairModule):
    "block of material"
 
    def __init__(self) -> None:
        ROOT.pyFairModule.__init__(self, self)
 
    def ConstructGeometry(self) -> None:
        print("Construct Block")
        top = ROOT.gGeoManager.GetTopVolume()
        geoLoad = ROOT.FairGeoLoader.Instance()
        geoFace = geoLoad.getGeoInterface()
        media = geoFace.getMedia()
        geoBuild = geoLoad.getGeoBuilder()
        ShipMedium = media.getMedium(material)
        W = ROOT.gGeoManager.GetMedium(material)
        if not W:
            geoBuild.createMedium(ShipMedium)
            W = ROOT.gGeoManager.GetMedium(material)
        aBox = ROOT.gGeoManager.MakeBox("target", W, 100.0 * u.cm, 100.0 * u.cm, thickness)
        top.AddNode(aBox, 1, ROOT.TGeoTranslation(0, 0, 0))
 
    def InitParContainers() -> None:
        print("not implemented!")
 
 
sensPlane = ROOT.exitHadronAbsorber()
sensPlane.SetEnergyCut(ecut * u.GeV)
sensPlane.SetZposition(thickness + 10 * u.cm)
run.AddModule(sensPlane)
target = Block()
run.AddModule(target)
 
primGen = ROOT.FairPrimaryGenerator()
myPgun = ROOT.FairBoxGenerator(2212, 1)  # pdg id  and multiplicity
myPgun.SetPRange(minmomentum, maxmomentum)
myPgun.SetPhiRange(0, 0)  # // Azimuth angle range [degree]
myPgun.SetThetaRange(0, 0)  # // Polar angle in lab system range [degree]
myPgun.SetXYZ(0.0 * u.cm, 0.0 * u.cm, -10.0 * u.cm - (thickness))
primGen.AddGenerator(myPgun)
run.SetGenerator(primGen)
#
run.SetGenerator(primGen)
# -----Initialize simulation run------------------------------------
run.Init()
gMC = ROOT.TVirtualMC.GetMC()
 
fStack = gMC.GetStack()
fStack.SetMinPoints(1)
fStack.SetEnergyCut(-1.0)
 
# -----Start run----------------------------------------------------
run.Run(nev)
 
# -----Start Analysis---------------
 
import rootUtils as ut
 
f = ROOT.TFile("TLV.root")
pdg = ROOT.TDatabasePDG.Instance()
h = {}
sTree = f.Get("cbmsim")
ut.bookHist(h, "Ekin", "Ekin of particles in sens plane", 400000, 0.0, 400)
ut.bookHist(h, "EkinLow", "Ekin of particles in sens plane", 1000, 0.0, 0.001)
for n in range(sTree.GetEntries()):
    rc = sTree.GetEvent(n)
    for aHit in sTree.vetoPoint:
        oTrack = sTree.MCTrack[aHit.GetTrackID()]
        M = pdg.GetParticle(oTrack.GetPdgCode()).Mass()
        Ekin = ROOT.TMath.Sqrt(aHit.GetPx() ** 2 + aHit.GetPy() ** 2 + aHit.GetPz() ** 2 + M**2) - M
        rc = h["Ekin"].Fill(Ekin)
        rc = h["EkinLow"].Fill(Ekin)
ut.bookCanvas(h, key=s, title=s, nx=900, ny=600, cx=1, cy=1)
tc = h[s].cd(1)
tc.SetLogy(1)
h["Ekin"].Draw()
 
# tungsten dedex mev cm**2/g 1.145 * rho g/cm-3 19.3 * 0.1
#          interactionLength 191.9 / 19.3 = 9.94cm, 0.1/9.94 = 1%
 
# -----Finish-------------------------------------------------------
timer.Stop()
rtime = timer.RealTime()
ctime = timer.CpuTime()
print(" ")
print("Macro finished successfully.")
print("Output file is ", outFile)
print("Real time ", rtime, " s, CPU time ", ctime, "s")