study_GammaConv.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 os
import sys
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 = 10  # 000000
 
setup = {}
setup["10"] = {"thickness": 2 * u.cm, "material": "aluminium", "momentum": 10 * u.GeV}
setup["100"] = {"thickness": 2 * u.cm, "material": "aluminium", "momentum": 100 * u.GeV}
setup["200"] = {"thickness": 2 * u.cm, "material": "aluminium", "momentum": 200 * u.GeV}
 
# 8cm = 0.9X0
 
 
s = sys.argv[1]
thickness = setup[s]["thickness"]
material = setup[s]["material"]
momentum = setup[s]["momentum"]
 
checkOverlap = True
 
outFile = "gconv" + s + ".root"
theSeed = int(10000 * time.time() % 10000000)
ecut = 0.0
 
# -------------------------------------------------------------------
ROOT.gRandom.SetSeed(theSeed)  # this should be propagated via ROOT to Pythia8 and Geant4VMC
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----------------------------------------
run = ROOT.FairRunSim()
run.SetName(mcEngine)  # Transport engine
run.SetSink(ROOT.FairRootFileSink(outFile))  # Output file
 
boostFactor = 100.0
if boostFactor > 1:
    # Turn off UseGeneralProcess to access GammaToMuons directly when cross-sections need to be changed
    os.environ["SET_GENERAL_PROCESS_TO_FALSE"] = "1"
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)
        self.sensPlane = None
 
    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))
        if self.sensPlane:
            self.sensPlane.AddSensitiveVolume(aBox)
 
    def InitParContainers() -> None:
        print("not implemented!")
 
    def makeSensitive(self, sensPlane) -> None:
        self.sensPlane = sensPlane
 
 
sensPlane = ROOT.exitHadronAbsorber()
sensPlane.SetEnergyCut(ecut * u.GeV)
sensPlane.SetZposition(thickness + 10 * u.cm)
run.AddModule(sensPlane)
target = Block()
target.makeSensitive(sensPlane)
run.AddModule(target)
 
primGen = ROOT.FairPrimaryGenerator()
myPgun = ROOT.FairBoxGenerator(22, 1)  # pdg id  and multiplicity
myPgun.SetPRange(momentum - 0.01, momentum + 0.01)
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)
 
if boostFactor > 1:
    ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"')
    ROOT.gROOT.ProcessLine('#include "Geant4/G4AnnihiToMuPair.hh"')
    ROOT.gROOT.ProcessLine('#include "Geant4/G4GammaConversionToMuons.hh"')
    gProcessTable = ROOT.G4ProcessTable.GetProcessTable()
    procAnnihil = gProcessTable.FindProcess(ROOT.G4String("AnnihiToMuPair"), ROOT.G4String("e+"))
    procGMuPair = gProcessTable.FindProcess(ROOT.G4String("GammaToMuPair"), ROOT.G4String("gamma"))
    procGMuPair.SetCrossSecFactor(boostFactor)
    procAnnihil.SetCrossSecFactor(boostFactor)
 
# -----Start run----------------------------------------------------
run.Run(nev)
 
# -----Start Analysis---------------
ROOT.gROOT.ProcessLine('#include "Geant4/G4EmParameters.hh"')
emP = ROOT.G4EmParameters.Instance()
emP.Dump()
 
import rootUtils as ut
 
f = ROOT.gROOT.GetListOfFiles()[0]
h = {}
ut.bookHist(h, "muons", "muon mult", 10, -0.5, 9.5)
ut.bookHist(h, "electrons", "e mult", 10, -0.5, 9.5)
sTree = f.Get("cbmsim")
for n in range(sTree.GetEntries()):
    rc = sTree.GetEvent(n)
    nMu = 0
    nEl = 0
    for aHit in sTree.vetoPoint:
        if abs(sTree.MCTrack[aHit.GetTrackID()].GetPdgCode()) == 13:
            nMu += 1
        if abs(sTree.MCTrack[aHit.GetTrackID()].GetPdgCode()) == 11:
            nEl += 1
    rc = h["muons"].Fill(nMu)
    rc = h["electrons"].Fill(nEl)
h["muons"].Draw()
 
# -----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")