Timing Correction for MuonDIS generation

CHANGELOG.md

@@ -16,7 +16,9 @@ it in future.
 
 ### Fixed
 
-* Use ConstructedAt instead of remove pythonization for TClonesArray 
+* Use ConstructedAt instead of remove pythonization for TClonesArray
+* **Corrections in MuonDIS simulation**
+  The DIS interactions are now time-shifted to be consistent with the original incoming muon. Additionally, tracks from soft interactions of the original muon along with the muon's veto response are preserved (in muonDis_<>.root) and included up to the DIS interaction point. To be noted that the muon veto points are manually added using add_muonresponse.py, which modifies the simulation file. This replaces the old method of "backward-travelling muon" to generate the incoming muon's veto response. All MuonDIS simulation scripts have been updated and consolidated within FairShip/muonDIS, ensuring consistency for new productions.
 
 ### Changed
 

macro/run_simScript.py

@@ -70,8 +70,6 @@
 }
 default = '2023'
 
-inactivateMuonProcesses = False   # provisionally for making studies of various muon background sources
-
 parser = ArgumentParser()
 group = parser.add_mutually_exclusive_group()
 parser.add_argument("--evtcalc", help="Use EventCalc", action="store_true")
@@ -150,6 +148,7 @@
     const="vacuums",
     default="helium"
 )
+
 parser.add_argument("--SND", dest="SND", help="Activate SND.", action='store_true')
 parser.add_argument("--noSND", dest="SND", help="Deactivate SND. NOOP, as it currently defaults to off.", action='store_false')
 
@@ -271,7 +270,6 @@
 # import shipMuShield_only as shipDet_conf # special use case for an attempt to convert active shielding geometry for use with FLUKA
 # import shipTarget_only as shipDet_conf
 import shipDet_conf
-
 modules = shipDet_conf.configure(run,ship_geo)
 # -----Create PrimaryGenerator--------------------------------------
 primGen = ROOT.FairPrimaryGenerator()
@@ -332,7 +330,6 @@
 # P8gen.SetMom(500.*u.GeV)
 # P8gen.SetId(-211)
  primGen.AddGenerator(P8gen)
-
 if simEngine == "FixedTarget":
  P8gen = ROOT.FixedTargetGenerator()
  P8gen.SetTarget("volTarget_1",0.,0.)
@@ -378,6 +375,11 @@
  ut.checkFileExists(inputFile)
  primGen.SetTarget(0., 0.)
  DISgen = ROOT.MuDISGenerator()
+ #  add Carbon to PDG
+ pdg = ROOT.TDatabasePDG.Instance()
+ pdg.AddParticle("C12", "Carbon-12", 12.0, True, 0, 6.0, "nucleus", 1000060120)
+ pdg.AddParticle("C13", "Carbon-13", 13.003355, True, 0, 6.0, "nucleus", 1000060130)
+
  # from nu_tau detector to tracking station 2
  # mu_start, mu_end =  ship_geo.tauMudet.zMudetC,ship_geo.TrackStation2.z
  #
@@ -388,7 +390,6 @@
  DISgen.Init(inputFile,options.firstEvent)
  primGen.AddGenerator(DISgen)
  options.nEvents = min(options.nEvents,DISgen.GetNevents())
- inactivateMuonProcesses = True # avoid unwanted hadronic events of "incoming" muon flying backward
  print('Generate ',options.nEvents,' with DIS input', ' first event',options.firstEvent)
 # -----neutrino interactions from nuage------------------------
 if simEngine == "Nuage":
@@ -518,15 +519,20 @@
  sys.exit(0)
 gMC = ROOT.TVirtualMC.GetMC()
 fStack = gMC.GetStack()
+if options.mudis:
+  EnergyCut=10.*u.MeV
+else:
+  EnergyCut=100.*u.MeV
+
 if MCTracksWithHitsOnly:
  fStack.SetMinPoints(1)
  fStack.SetEnergyCut(-100.*u.MeV)
 elif MCTracksWithEnergyCutOnly:
  fStack.SetMinPoints(-1)
- fStack.SetEnergyCut(100.*u.MeV)
+ fStack.SetEnergyCut(EnergyCut)
 elif MCTracksWithHitsOrEnergyCut:
  fStack.SetMinPoints(1)
- fStack.SetEnergyCut(100.*u.MeV)
+ fStack.SetEnergyCut(EnergyCut)
 elif options.deepCopy:
  fStack.SetMinPoints(0)
  fStack.SetEnergyCut(0.*u.MeV)
@@ -560,18 +566,6 @@
 #fieldMaker.plotField(1, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-300.0, 300.0, 6.0), 'Bzx.png')
 #fieldMaker.plotField(2, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-400.0, 400.0, 6.0), 'Bzy.png')
 
-if inactivateMuonProcesses :
- ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"')
- mygMC = ROOT.TGeant4.GetMC()
- mygMC.ProcessGeantCommand("/process/inactivate muPairProd")
- mygMC.ProcessGeantCommand("/process/inactivate muBrems")
- mygMC.ProcessGeantCommand("/process/inactivate muIoni")
- mygMC.ProcessGeantCommand("/process/inactivate muonNuclear")
- mygMC.ProcessGeantCommand("/particle/select mu+")
- mygMC.ProcessGeantCommand("/particle/process/dump")
- gProcessTable = ROOT.G4ProcessTable.GetProcessTable()
- procmu = gProcessTable.FindProcess(ROOT.G4String('muIoni'),ROOT.G4String('mu+'))
- procmu.SetVerboseLevel(2)
 # -----Start run----------------------------------------------------
 run.Run(options.nEvents)
 # -----Runtime database---------------------------------------------

muonDIS/add_muonresponse.py

@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+"""Script to add the incoming muon's MC hits in the SBT to the simulation file."""
+
+import argparse
+import logging
+
+import ROOT as r
+
+parser = argparse.ArgumentParser(description=__doc__)
+parser.add_argument(
+    "-f",
+    "--inputfile",
+    required=True,
+    help="Simulation file produced from DIS, WARNING: File will be modified.",
+)
+parser.add_argument(
+    "-m",
+    "--muonfile",
+    required=True,
+    help="Muon file used for DIS production (muonDis_<>.root)",
+)
+args = parser.parse_args()
+
+
+def inspect_file(inputfile):
+    """Inspecting the produced file for successfully added muon veto points."""
+    input_file = r.TFile.Open(inputfile, "READ")
+    input_tree = input_file.cbmsim
+    input_tree.GetEntry(0)
+
+    muons = False
+    for hit in input_tree.vetoPoint:
+        if hit.GetTrackID() == 0:
+            muons = True
+    input_file.Close()
+
+    if muons:
+        print("File is already updated with incoming muon info. Nothing to do.")
+        return False
+    else:
+        print(
+            "Incoming muon's vetopoint info missing in file, proceeding with modification"
+        )
+        return True
+
+
+def modify_file(inputfile, muonfile):
+    """Add information from original muon to input simulation file."""
+    logging.warning(f"{inputfile} will be opened for modification")
+
+    input_file = r.TFile.Open(inputfile, "UPDATE")
+    try:
+        input_tree = input_file.cbmsim
+    except Exception as e:
+        print(f"Error: {e}")
+        input_file.Close()
+        exit(1)
+
+    # Open the external file with additional vetoPoints
+    muon_file = r.TFile.Open(muonfile, "READ")
+    try:
+        muon_tree = muon_file.DIS
+    except Exception as e:
+        print(f"Error: {e}")
+        muon_file.Close()
+        exit(1)
+
+    input_file.cd()
+    output_tree = input_tree.CloneTree(
+        0
+    )  # Clone the structure of the existing tree, but do not copy the entries
+
+    # Access the vetoPoint branch in the original and external trees
+    original_vetoPoint = r.TClonesArray("vetoPoint")
+    muon_vetoPoint = r.TClonesArray("vetoPoint")
+    combined_vetoPoint = r.TClonesArray("vetoPoint")
+
+    input_tree.SetBranchAddress("vetoPoint", original_vetoPoint)
+    muon_tree.SetBranchAddress("muon_vetoPoints", muon_vetoPoint)
+    output_tree.SetBranchAddress("vetoPoint", combined_vetoPoint)
+
+    for input_event, muon_event in zip(input_tree, muon_tree):
+        interaction_point = r.TVector3()
+        input_event.MCTrack[0].GetStartVertex(interaction_point)
+
+        combined_vetoPoint.Clear()
+
+        index = 0
+
+        for hit in original_vetoPoint:
+            if combined_vetoPoint.GetSize() == index:
+                combined_vetoPoint.Expand(index + 1)
+            combined_vetoPoint[index] = hit  # pending fix to support ROOT 6.32+
+            index += 1
+
+        for hit in muon_vetoPoint:
+            if hit.GetZ() < interaction_point.Z():
+                if combined_vetoPoint.GetSize() == index:
+                    combined_vetoPoint.Expand(index + 1)
+                combined_vetoPoint[index] = hit  # pending fix to support ROOT 6.32+
+                index += 1
+
+        output_tree.Fill()
+
+    # Write the updated tree back to the same file
+    input_file.cd()
+    output_tree.Write("cbmsim", r.TObject.kOverwrite)
+    input_file.Close()
+    muon_file.Close()
+
+    print(f"File updated with incoming muon info as {inputfile}")
+
+
+add_muons = inspect_file(args.inputfile)
+
+if add_muons:
+    modify_file(args.inputfile, args.muonfile)