Muon DIS with Pythia6 run through Geant4

CMakeLists.txt

@@ -257,7 +257,6 @@ IF(NOT FairRoot_FOUND)
 Message("--- FairRoot not found ...")
 ENDIF (NOT FairRoot_FOUND)
 
-add_subdirectory (sndFairTasks)
 add_subdirectory (shipLHC)
 add_subdirectory (shipdata)
 add_subdirectory (passive)
@@ -267,6 +266,8 @@ add_subdirectory (field)
 add_subdirectory (genfit)
 add_subdirectory (millepede)
 add_subdirectory (analysis)
+add_subdirectory (G4Processes)
+add_subdirectory (sndFairTasks)
 
 FIND_PATH(TEvePath NAMES TEveEventManager.h PATHS
   ${SIMPATH}/tools/root/include

G4Processes/CMakeLists.txt

@@ -0,0 +1,27 @@
+cmake_minimum_required(VERSION 3.12)
+project(G4Processes)
+
+# Find GEANT4 and ROOT
+find_package(Geant4 REQUIRED)
+find_package(ROOT REQUIRED)
+
+# Create library
+add_library(G4Processes SHARED G4MuonDISProcess.cc G4MuonDISProcess.hh ${CMAKE_SOURCE_DIR}/shipgen/MuDISGenerator.cxx)
+
+# Include headers, also from GEANT4 and ROOT
+target_include_directories(G4Processes PUBLIC
+  ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_SOURCE_DIR}/shipgen
+  ${GEANT4_INCLUDE_DIR}
+  ${ROOT_INCLUDE_DIR}
+)
+
+target_link_libraries(G4Processes PUBLIC 
+  ${GEANT4_LIBRARIES}     # Link GEANT4 libraries
+  ${ROOT_LIBRARIES}       # Link ROOT libraries
+  ${ROOT_LIBRARY_DIR}/libEGPythia6.so
+  ${CMAKE_SOURCE_DIR}/shipgen # Link shipgen
+)
+
+# Install the library
+install(TARGETS G4Processes DESTINATION lib)

G4Processes/G4MuonDISProcess.hh

@@ -0,0 +1,88 @@
+#pragma once
+
+#ifndef G4MuonDISProcess_h
+#define G4MuonDISProcess_h
+
+#include "G4VDiscreteProcess.hh"
+#include "G4Track.hh"
+#include "G4Step.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4VParticleChange.hh"
+#include "G4ThreeVector.hh"
+#include "globals.hh"
+#include <tuple>
+#include <vector>
+#include <map>
+
+class MuDISGenerator;
+class G4Navigator;
+class G4PhysicsFreeVector;
+
+class G4MuonDISProcess : public G4VDiscreteProcess {
+public:
+   /** Constructors **/
+   G4MuonDISProcess(const G4String &name = "muonDIS");
+   virtual ~G4MuonDISProcess() {}
+
+   /** Called to check if this process should apply to a given particle **/
+   virtual G4bool IsApplicable(const G4ParticleDefinition &particle) override;
+
+   virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize,
+                                                         G4ForceCondition *condition) override;
+
+   /**  Called to calculate the mean free path **/
+   virtual G4double
+   GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override;
+
+   /**  Actual interaction **/
+   virtual G4VParticleChange *PostStepDoIt(const G4Track &track, const G4Step &step) override;
+
+   /**  Set the nucleon type **/
+   void SetNucleonType(char *nucleon) { fNucleon = nucleon; }
+   /** Set the x-y-z range of the DIS interaction **/
+   void SetRange(std::vector<float> *x_range, std::vector<float> *y_range, std::vector<float> *z_range);
+   /** Set the name of the geometry volume selected for placement of the DIS interaction **/
+   void SetVolume(char* volume) { fVolumeName = volume; }
+   /**  Set the DIS cross section lookup tables for the set nucleon type  **/
+   void SetCrossSectionTables(std::shared_ptr<std::map<int, G4PhysicsFreeVector*>> xsecTables) { fXsecTables = xsecTables; }
+
+private:
+   /**  Pythia6 call to generate the DIS process **/
+   void SimulateDIS(const G4Track &track);
+
+   /** Reset private class member variables **/
+   void ResetParameters();
+
+   /** Rotation of Pythia6 outgoing particles to match physics CS **/
+   std::tuple<double, double, double>
+   rotate(double ctheta, double stheta, double cphi, double sphi, double px, double py, double pz);
+
+   /*** 3D linear track extrapolation **/
+   G4ThreeVector linear_track_extrapolation(G4ThreeVector pos, G4ThreeVector mom, G4double z_ref);
+
+   MuDISGenerator *disGen; /// class instance used to calculate material budget
+
+   G4Navigator *g4Navigator; /// Geant4 geo navigator
+
+   G4double fCrossSection;    /// DIS cross section from the external generator
+   G4double fOutTracksWeight; /// Weight of outgoing particles
+   char *fNucleon;            //! nucleon type
+   std::shared_ptr<std::map<int, G4PhysicsFreeVector*>> fXsecTables; /// DIS cross section lookup tables
+   char *fVolumeName;         //! geometry volume name to place the DIS
+   G4double fXStart{};        /// start of x range to place the DIS
+   G4double fXEnd{};          /// end of x range to place the DIS
+   G4double fYStart{};        /// start of y range to place the DIS
+   G4double fYEnd{};          /// end of y range to place the DIS
+   G4double fZStart{};        /// start of z range to place the DIS
+   G4double fZEnd{};          /// end of z range to place the DIS
+   // The following are items that update on a per-event basis and not for each step!
+   G4bool was_simulated = 0;    /// true if the DIS was already simulated, false otherwise
+   double bparam{};             /// mean density in a range
+   double mparam[8]{};          //[8] density parameters along track
+   double max_density{};        /// maximum density in the whole range
+   double prev_max_density{};   /// maximum density in the path traversed so far
+   double density_along_path{}; /// mean density in a range
+   double prev_start[3]{};      //[3] track position at previous step in the [z_start, z_end] range
+};
+
+#endif

shipLHC/Floor.cxx

@@ -109,7 +109,11 @@ void Floor::PreTrack(){
     }
     TLorentzVector  mom;
     gMC->TrackMomentum(mom);
-    if  ( (mom.E()-mom.M() )<fEmin && pos.Z()<fzPos){
+    // The preset energy cut on produced particles is applied when:
+    //    FastMuon option is OFF and pos.z < zPos
+    //    OR
+    //    FastMuon is ON. Then Ecut is valid in the whole z range.
+    if  ( (mom.E()-mom.M() )<fEmin && (fFastMuon || pos.Z()<fzPos)){
       gMC->StopTrack();
       return;
     }

shipLHC/muonDis.py

@@ -324,7 +324,7 @@ def getPythia6CrossSec(nstat,pmom=[]):
 
 Float_t fixXsec(TPythia6& g) {
     Pyint5_t* p5 = g.GetPyint5();
-    return p5->XSEC[0][3];
+    return p5->XSEC[2][0];
 }
 """)
 
@@ -2315,7 +2315,11 @@ def missing3Dproj(hist,ymin,ymax):
 elif not options.command.find("convert")<0: convertAscii2Root(options.muonIn)
 elif not options.command.find("make")<0:    makeMuDISEvents(nucleon=options.nucleon)
 elif not options.command.find("ana")<0:     analyze(options.muonIn)
-elif not options.command.find("cross")<0:   getPythiaCrossSec(options.nEvents)
+elif not options.command.find("cross")<0:   
+     if options.pythia6:
+         getPythia6CrossSec(options.nEvents)
+     else:
+         getPythia8CrossSec(options.nEvents) 
 elif not options.command.find("muonPreTransport")<0:    muonPreTransport()
 elif not options.command.find("muon")<0:    muonRateAtSND()
 

shipLHC/run_simSND.py

@@ -48,6 +48,15 @@
 parser.add_argument('--zMax', type=float, help="max distance to apply energy cut", dest='zmax', default=70000.)
 parser.add_argument("--Nuage",     dest="nuage",  help="Use Nuage, neutrino generator of OPERA", required=False, action="store_true")
 parser.add_argument("--MuDIS",     dest="mudis",  help="Use muon deep inelastic scattering generator", required=False, action="store_true")
+
+parser.add_argument("--MuDISG4",     dest="mudisg4",  help="Use external muon deep inelastic scattering generator with Geant4", required=False, action="store_true")
+parser.add_argument('--nucleon', choices=['p+', 'n'], help='Type of nucleon')
+parser.add_argument('--x_range', help='x position range placing muonDIS interaction as two values: start and end', type=float, nargs=2, default=[-1e10,1e10])
+parser.add_argument('--y_range', help='y position range placing muonDIS interaction as two values: start and end', type=float, nargs=2, default=[-1e10,1e10])
+parser.add_argument('--z_range', help='z position range placing muonDIS interaction as two values: start and end', type=float, nargs=2, default=[-1e10,1e10])
+parser.add_argument('--volume',  help='Name of geometry volume where to simulate the DIS; e.g. volTarget, volMuUpstreamDet, volMuDownstreamDet, volVetoPlane, Wall', default="")
+parser.add_argument('--dis_xsec_file', help='Name of DIS cross section file', default="/eos/experiment/sndlhc/MonteCarlo/Pythia6/MuonDIS/muDIScrossSec.root")
+
 parser.add_argument("-n", "--nEvents",dest="nEvents",  help="Number of events to generate", required=False,  default=100, type=int)
 parser.add_argument("-i", "--firstEvent",dest="firstEvent",  help="First event of input file to use", required=False,  default=0, type=int)
 parser.add_argument("-s", "--seed",dest="theSeed",  help="Seed for random number. Only for experts, see TRrandom::SetSeed documentation", required=False,  default=0, type=int)
@@ -67,6 +76,9 @@
 # user hook
 userTask = False
 
+# make log colourful
+ROOT.FairLogger.GetLogger().SetColoredLog(True)
+
 class MyTask(ROOT.FairTask):
     "user task"
 
@@ -78,7 +90,9 @@ def Exec(self,opt):
         if MCTracks.GetEntries()>100:  fMC.StopRun()
 
 checking4overlaps = False
-if options.debug: checking4overlaps = True
+if options.debug: 
+  checking4overlaps = True
+  ROOT.fair.Logger.SetConsoleSeverity("debug")
 
 if options.pythia8:       simEngine = "Pythia8"
 if options.pg:                 simEngine = "PG"
@@ -197,7 +211,8 @@ def Exec(self,opt):
         f'({options.EVx},{options.EVy},{options.EVz})[cm × cm × cm] \n'
         f'with a uniform x-y spread of (Dx,Dy)=({options.Dx},{options.Dy})[cm × cm]'
         f' and {options.nZSlices} z slices in steps of {options.zSliceStep}[cm].')
-  ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # otherwise stupid printout for each event
+  if not options.debug:
+    ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # otherwise stupid printout for each event
 # -----muon DIS Background------------------------
 if simEngine == "muonDIS":
    ut.checkFileExists(inputFile)
@@ -294,6 +309,10 @@ def Exec(self,opt):
      elif 'Floor' in modules: 
            modules['Floor'].SetFastMuon()
            modules['Floor'].SetZmax(options.zmax)
+           modules['Floor'].SetEmin(options.ecut)
+           if (options.ecut!=-1.):
+              print(f"Applying energy cut of {options.ecut} in full z range!\
+                    \nThe set zMax will be used for the FastMuon option only!")
            print('transport only-muons up to z=',options.zmax)
 # ------------------------------------------------------------------------
 #---Store the visualiztion info of the tracks, this make the output file very large!!
@@ -305,6 +324,26 @@ def Exec(self,opt):
 
 # -----Initialize simulation run------------------------------------
 run.Init()
+
+# -----muon DIS Background: add muonDIS from an external generator as a process in GEANT4
+# Only after Geant4 is initialized can one add a process to the simulation,
+# thus one adds the custom DIS after run.Init()
+if options.mudisg4:
+  # Make sure ranges run from lower to higher limit
+  options.x_range = sorted(options.x_range)
+  options.y_range = sorted(options.y_range)
+  options.z_range = sorted(options.z_range)
+  # add a process injector task for the Geant4 propagation
+  external_muDIS_injector = ROOT.MuonDISProcessInjector(options.nucleon,
+                                                        options.x_range, 
+                                                        options.y_range,
+                                                        options.z_range,
+                                                        options.volume,
+                                                        options.dis_xsec_file)
+  external_muDIS_injector.Init()
+  run.AddTask(external_muDIS_injector)
+#------------------------------------------------------------------------
+
 gMC = ROOT.TVirtualMC.GetMC()
 fStack = gMC.GetStack()
 if MCTracksWithHitsOnly:
@@ -379,7 +418,6 @@ def Exec(self,opt):
  procmu = gProcessTable.FindProcess(ROOT.G4String('muIoni'),ROOT.G4String('mu+'))
  procmu.SetVerboseLevel(2)
 
-if options.debug:  ROOT.fair.Logger.SetConsoleSeverity("debug")
 # -----Start run----------------------------------------------------
 run.Run(options.nEvents)
 # -----Runtime database---------------------------------------------

shipgen/MuDISGenerator.cxx

@@ -83,6 +83,7 @@ Double_t MuDISGenerator::MeanMaterialBudget(const Double_t *start, const Double_
   //
   Double_t bparam[6]; // total parameters
   Double_t lparam[6]; // local parameters
+  Double_t tiny_step = 1e-6;
 
   for (Int_t i=0;i<6;i++) bparam[i]=0;
 
@@ -150,6 +151,19 @@ Double_t MuDISGenerator::MeanMaterialBudget(const Double_t *start, const Double_
     currentnode = gGeoManager->GetCurrentNode();
     if (snext<2.*TGeoShape::Tolerance()) nzero++;
     else nzero = 0;
+    // If one is consistently stuck at some point,
+    // reset the navigator and move forward with a tiny_step
+    if (nzero>1){
+      gGeoManager->CdTop();
+      gGeoManager->FindNode(gGeoManager->GetCurrentPoint()[0] - dir[0] * tiny_step / dir[2],
+                            gGeoManager->GetCurrentPoint()[1] - dir[1] * tiny_step / dir[2],
+                            gGeoManager->GetCurrentPoint()[2] + tiny_step);
+      snext = std::max(snext, tiny_step);
+      nzero = 0;
+    }
+    /*
+    // The `if (nzero>1)` above makes sure one doesn't reach here.
+    // Keeping this for reference only.
     if (nzero>3) {
       // This means navigation has problems on one boundary
       // Try to cross by making a small step
@@ -164,6 +178,7 @@ Double_t MuDISGenerator::MeanMaterialBudget(const Double_t *start, const Double_
       mparam[1] = 1000000;        // and infinite rad length
       return bparam[0]/step;
     }
+    */
     mparam[6]+=1.;
     step += snext;
     bparam[1]    += snext/lparam[1];

shipgen/MuDISGenerator.h

@@ -33,7 +33,6 @@ class MuDISGenerator : public FairGenerator
     endZ        = zE; 
   }
 
- private:
   Double_t MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam);
 
 

sndFairTasks/CMakeLists.txt

@@ -10,6 +10,9 @@ ${CMAKE_SOURCE_DIR}/sndFairTasks
 ${CMAKE_SOURCE_DIR}/veto
 ${XROOTD_INCLUDE_DIR}
 ${XROOTD_INCLUDE_DIR}/..
+${GEANT4_ROOT}/include/Geant4
+${GEANT4_VMC_ROOT}/include/geant4vmc
+${CMAKE_SOURCE_DIR}/G4Processes
 )
 
 include_directories( ${INCLUDE_DIRECTORIES}  ${VMC_INCLUDE_DIRS}  ${FAIRROOT_INCLUDE_DIR} ${FairLogger_INCDIR} ${FMT_INCLUDE_DIR}
@@ -20,7 +23,7 @@ set(LINK_DIRECTORIES
 ${ROOT_LIBRARY_DIR}
 ${XROOTD_LIBRARY_DIR}
 ${FAIRROOT_LIBRARY_DIR}
-
+${G4Processes_LIBRARY_DIR}
 )
 
 link_directories( ${LINK_DIRECTORIES})
@@ -30,11 +33,12 @@ DigiTaskSND.cxx
 ConvRawData.cxx
 boardMappingParser.cxx
 MCEventBuilder.cxx
+MuonDISProcessInjector.cxx
 )
 
 Set(HEADERS)
 Set(LINKDEF sndFairTasksLinkDef.h)
 Set(LIBRARY_NAME sndFairTasks)
-Set(DEPENDENCIES Base ShipData shipLHC GeoBase ParBase Geom Core)
+Set(DEPENDENCIES Base ShipData shipLHC GeoBase ParBase Geom Core G4Processes)
 
 GENERATE_LIBRARY()