From a6fec0b06bd91918ba6e8674f79f5e9a036d0e34 Mon Sep 17 00:00:00 2001 From: siilieva Date: Thu, 10 Jul 2025 14:35:49 +0200 Subject: [PATCH 1/8] Add MuonDISProcessInjector task It adds to the Geant4 process manager a custom process that is using Pythia6 to generate muonDIS. --- sndFairTasks/CMakeLists.txt | 8 ++- sndFairTasks/MuonDISProcessInjector.cxx | 86 +++++++++++++++++++++++++ sndFairTasks/MuonDISProcessInjector.h | 33 ++++++++++ sndFairTasks/sndFairTasksLinkDef.h | 5 +- 4 files changed, 126 insertions(+), 6 deletions(-) create mode 100644 sndFairTasks/MuonDISProcessInjector.cxx create mode 100644 sndFairTasks/MuonDISProcessInjector.h diff --git a/sndFairTasks/CMakeLists.txt b/sndFairTasks/CMakeLists.txt index ef4d3f8bd8..d66ce569e2 100644 --- a/sndFairTasks/CMakeLists.txt +++ b/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() diff --git a/sndFairTasks/MuonDISProcessInjector.cxx b/sndFairTasks/MuonDISProcessInjector.cxx new file mode 100644 index 0000000000..3cc763ef57 --- /dev/null +++ b/sndFairTasks/MuonDISProcessInjector.cxx @@ -0,0 +1,86 @@ +#include "MuonDISProcessInjector.h" +#include "G4MuonMinus.hh" +#include "G4MuonPlus.hh" +#include "G4ProcessManager.hh" +#include "G4MuonDISProcess.hh" +#include "FairLogger.h" +#include "TFile.h" +#include "TGraph.h" +#include +#include + +using std::vector; + +// Convert ROOT's TGraph to G4PhysicsFreeVector +G4PhysicsFreeVector* GraphToFreeVector(TGraph* graph) { + auto xsec_data = new G4PhysicsFreeVector(); + for (int i = 0; i < graph->GetN(); ++i) { + double x, y; + graph->GetPoint(i, x, y); + xsec_data->InsertValues(x, y); + } + return xsec_data; +} + +MuonDISProcessInjector::MuonDISProcessInjector(char *nucleon, vector x_range, vector y_range, + vector z_range, char *volume, char *xsec_filename) + : FairTask("MuonDISProcessInjector") +{ + fNucleon = nucleon; + fVolumeName = volume; + fXRange = x_range; + fYRange = y_range; + fZRange = z_range; + + TFile* xsec_file = new TFile(xsec_filename); + if(!xsec_file || xsec_file->IsZombie()){ + LOG(FATAL) << "DIS cross section file not found"; + exit(0); + } + TGraph* crsec_muminus = static_cast (xsec_file->Get(Form("g_13_%s", fNucleon))); + TGraph* crsec_muplus = static_cast (xsec_file->Get(Form("g_-13_%s", fNucleon))); + if (!crsec_muminus || !crsec_muplus){ + LOG(FATAL) << "DIS cross section graphs per " << fNucleon <<" are missing in "<< xsec_filename; + exit(0); + } + + fXsecTables = std::make_shared>(); + fXsecTables->insert({13, GraphToFreeVector(crsec_muminus)}); + fXsecTables->insert({-13, GraphToFreeVector(crsec_muplus)}); + + LOG(WARNING) << "MuonDISProcessInjector: Setting xyz ranges[cm] for muon DIS\nx: " + << fXRange[0] << ", " << fXRange[1] << "\ny: " << fYRange[0] << ", " << fYRange[1] << "\nz: " + << fZRange[0] << ", " << fZRange[1] << "\nand volume name '" << fVolumeName + << "' and nucleon type " << fNucleon; +} + +InitStatus MuonDISProcessInjector::Init() +{ + muMinus = G4MuonMinus::Definition(); + muPlus = G4MuonPlus::Definition(); + + DISProcess = new G4MuonDISProcess(); + + // Set the nucleon type + DISProcess->SetNucleonType(fNucleon); + // Set DIS xsec lookup tables per set nucleon type + DISProcess->SetCrossSectionTables(fXsecTables); + // Set the z range of the DIS interaction + DISProcess->SetRange(&fXRange, &fYRange, &fZRange); + // Set the name of the geo volume for the DIS interaction + DISProcess->SetVolume(fVolumeName); + + auto process_man_muminus = muMinus->GetProcessManager(); + if (process_man_muminus) { + G4int id_muminus = process_man_muminus->AddDiscreteProcess(DISProcess); + } + auto process_man_muplus = muPlus->GetProcessManager(); + if (process_man_muplus) { + G4int id_muplus = process_man_muplus->AddDiscreteProcess(DISProcess); + } + + LOG(WARNING) << "Adding the external generation for muon DIS using Pythia6"; + return kSUCCESS; +} + +ClassImp(MuonDISProcessInjector) diff --git a/sndFairTasks/MuonDISProcessInjector.h b/sndFairTasks/MuonDISProcessInjector.h new file mode 100644 index 0000000000..733223aed3 --- /dev/null +++ b/sndFairTasks/MuonDISProcessInjector.h @@ -0,0 +1,33 @@ +#ifndef MUONDISPROCESSINJECTOR_H +#define MUONDISPROCESSINJECTOR_H + +#include "FairTask.h" +#include "G4ParticleDefinition.hh" +#include "G4MuonDISProcess.hh" +#include +#include + +class G4PhysicsFreeVector; + +class MuonDISProcessInjector : public FairTask { +public: + MuonDISProcessInjector(char *nucleon, std::vector x_range, std::vector y_range, + std::vector z_range, char *volume, char *xsec_filename); + virtual ~MuonDISProcessInjector() {} + + virtual InitStatus Init(); + + char *fNucleon; //! nucleon type + char *fVolumeName; //! geometry volume name to place the DIS + std::vector fXRange; /// x range to place the DIS + std::vector fYRange; /// y range to place the DIS + std::vector fZRange; /// z range to place the DIS + std::shared_ptr> fXsecTables; /// DIS cross section lookup tables + const G4ParticleDefinition *muPlus; + const G4ParticleDefinition *muMinus; + G4MuonDISProcess *DISProcess; + + ClassDef(MuonDISProcessInjector, 1); +}; + +#endif diff --git a/sndFairTasks/sndFairTasksLinkDef.h b/sndFairTasks/sndFairTasksLinkDef.h index 7ea592291b..4b494fe1db 100755 --- a/sndFairTasks/sndFairTasksLinkDef.h +++ b/sndFairTasks/sndFairTasksLinkDef.h @@ -7,8 +7,5 @@ #pragma link C++ class DigiTaskSND; #pragma link C++ class ConvRawData; #pragma link C++ class MCEventBuilder; +#pragma link C++ class MuonDISProcessInjector; #endif - - - - From 3d32a5305856b698228c5c82502e6a4f93cceb52 Mon Sep 17 00:00:00 2001 From: siilieva Date: Thu, 10 Jul 2025 14:44:39 +0200 Subject: [PATCH 2/8] Make the MeanMaterialBudget calculation public to reuse the function globally in sndsw. --- shipgen/MuDISGenerator.h | 1 - 1 file changed, 1 deletion(-) diff --git a/shipgen/MuDISGenerator.h b/shipgen/MuDISGenerator.h index 6d6dd4606e..1e3a5edaa4 100644 --- a/shipgen/MuDISGenerator.h +++ b/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); From f128a66edeeffbe01c3823903011b465d8e26fe6 Mon Sep 17 00:00:00 2001 From: siilieva Date: Thu, 10 Jul 2025 15:48:15 +0200 Subject: [PATCH 3/8] run_simSND: Add options to run muonDIS with Geant4 --- shipLHC/run_simSND.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 85d5e19778..60b78420d1 100644 --- a/shipLHC/run_simSND.py +++ b/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) @@ -305,6 +314,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: From 877e6bda067e73ec7f73403ac012fe2f40ef63af Mon Sep 17 00:00:00 2001 From: siilieva Date: Thu, 10 Jul 2025 15:50:31 +0200 Subject: [PATCH 4/8] Add muon DIS generated my Pythia6 as a Geant4 process Using standalone CMake and don't rely on the FairRoot CMake macros as we usually do. --- CMakeLists.txt | 3 +- G4Processes/CMakeLists.txt | 27 +++ G4Processes/G4MuonDISProcess.cc | 386 ++++++++++++++++++++++++++++++++ G4Processes/G4MuonDISProcess.hh | 88 ++++++++ 4 files changed, 503 insertions(+), 1 deletion(-) create mode 100644 G4Processes/CMakeLists.txt create mode 100644 G4Processes/G4MuonDISProcess.cc create mode 100644 G4Processes/G4MuonDISProcess.hh diff --git a/CMakeLists.txt b/CMakeLists.txt index b54c3eed16..fdaa51ab88 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/G4Processes/CMakeLists.txt b/G4Processes/CMakeLists.txt new file mode 100644 index 0000000000..1cf355d7b4 --- /dev/null +++ b/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) diff --git a/G4Processes/G4MuonDISProcess.cc b/G4Processes/G4MuonDISProcess.cc new file mode 100644 index 0000000000..b6391ed6f6 --- /dev/null +++ b/G4Processes/G4MuonDISProcess.cc @@ -0,0 +1,386 @@ +#include "G4MuonDISProcess.hh" +#include "G4SystemOfUnits.hh" +#include "G4PhysicalConstants.hh" +#include "G4ParticleTable.hh" +#include "G4MuonMinus.hh" +#include "G4MuonPlus.hh" +#include "G4StepPoint.hh" +#include "G4TrackStatus.hh" +#include "G4DynamicParticle.hh" +#include "G4ParticleChange.hh" +#include "MuDISGenerator.h" +#include "G4RandomTools.hh" +#include "G4EventManager.hh" +#include "G4ThreeVector.hh" +#include "G4RunManager.hh" +#include "G4TransportationManager.hh" +#include "G4Navigator.hh" +#include "G4ProcessType.hh" +#include "G4PhysicsFreeVector.hh" +#include "globals.hh" + +#include "TGeoManager.h" +#include "TGeoNode.h" +#include "TPythia6.h" +#include "TRandom.h" +#include "TMath.h" +#include "TROOT.h" + +#include +#include +#include +#include +#include +#include +#include + +using std::vector; +using std::string; + +G4MuonDISProcess::G4MuonDISProcess(const G4String &name) : G4VDiscreteProcess(name) +{ + if (verboseLevel > 0) { + std::cout << "G4MuonDISProcess created: " << name << " " << fCrossSection << std::endl; + } +} + +void G4MuonDISProcess::ResetParameters(){ + was_simulated = 0; + bparam = 0.; + std::memset(mparam, 0, sizeof(mparam)); + max_density = 0.; + prev_max_density = 0.; + density_along_path = 0.; + std::memset(prev_start, 0, sizeof(prev_start)); +} + +void G4MuonDISProcess::SetRange(vector *x_range, vector *y_range, vector *z_range) +{ + // The base length unit in sndsw is cm, but mm in Pythia 6 and GEANT4 + fXStart = x_range->at(0) * cm; + fXEnd = x_range->at(1) * cm; + fYStart = y_range->at(0) * cm; + fYEnd = y_range->at(1) * cm; + fZStart = z_range->at(0) * cm; + fZEnd = z_range->at(1) * cm; + LOG(INFO) << "G4MuonDISProcess: Setting xyz ranges[cm] for muon DIS\nx: " + << fXStart/cm << ", " << fXEnd/cm << "\ny: " << fYStart/cm << ", " << fYEnd/cm << "\nz: " + << fZStart/cm << ", " << fZEnd/cm; +} + +G4bool G4MuonDISProcess::IsApplicable(const G4ParticleDefinition &particle) +{ + return (&particle == G4MuonMinus::Definition() || &particle == G4MuonPlus::Definition()); +} + +G4VParticleChange *G4MuonDISProcess::PostStepDoIt(const G4Track &track, const G4Step &step) +{ + aParticleChange.Initialize(track); + G4ThreeVector pos = track.GetPosition(); + if (track.GetParentID() == 0 && was_simulated == 0 && (fZEnd - pos.z()) > 1e-3) { + if ((fZStart - pos.z()) > 1e-3) { + // do nothing + } else { + string geoVolumePath =gGeoManager->GetPath(); + LOG(DEBUG) << "Proposing DIS in volume path: " << geoVolumePath; + // Check if the track position is in the specified x-y range, z is surely in the z_range, + // and in the desired geometry volume (e.g. the target), if such is set + if (pos.x() >= fXStart && pos.x() <= fXEnd && pos.y() >= fYStart && pos.y() <= fYEnd &&\ + geoVolumePath.find(fVolumeName) != string::npos) { + LOG(DEBUG) << "DIS vertex in the x-y-z range, proceeding to DIS generation at "<< pos/cm + << "\nin pre-set geometry volume "<< fVolumeName; + // Simulate the DIS interaction - call Pythia6 + SimulateDIS(track); + } + // To save on computing, drop the whole event shall the selected DIS interaction point be + // outside the desired x-y range! + else { + LOG(DEBUG) << "DIS vertex out of x-y-z range and/or volume, aborting the event at " << pos/cm; + ResetParameters(); + G4RunManager::GetRunManager()->AbortEvent(); + } + } + } + return &aParticleChange; +} + +G4double G4MuonDISProcess::GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) +{ + return DBL_MAX; +} + +G4double G4MuonDISProcess::PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, + G4ForceCondition *condition) +{ + if ((previousStepSize < 0.0) || (theNumberOfInteractionLengthLeft <= 0.0)) { + // beggining of tracking (or just after DoIt of this process) + ResetNumberOfInteractionLengthLeft(); + } else if (previousStepSize > 0.0) { + // subtract NumberOfInteractionLengthLeft + SubtractNumberOfInteractionLengthLeft(previousStepSize); + } else { + // zero step + } + + G4ThreeVector pos = track.GetPosition(); + G4ThreeVector mom = track.GetMomentum(); + g4Navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4ThreeVector extrap_pos; + if (track.GetParentID() == 0){ + LOG(DEBUG) << "primary muons's position[cm]: " << pos/cm << " and momentum[GeV] " << mom.mag()/GeV; + } + // don't bother simulating DIS if the energy of the primary is below 5 GeV + if (track.GetParentID() != 0 || + (track.GetParentID() == 0 && (pos.z() > fZEnd || track.GetTotalEnergy() < 5. * GeV))) { + *condition = NotForced; + currentInteractionLength = DBL_MAX; // make sure this process won't be picked + // reset parameters once primary particle passes z_end and/or we're done simulating DIS for the event + if (max_density != 0 || density_along_path != 0 || was_simulated == 1) { + ResetParameters(); + } + } else if ((fZStart - pos.z()) > 2e-3) { + // Make sure the muon eventually lands at fZStart, use tiny tolerances. + currentInteractionLength = fZStart - pos.z() - 1e-3; + LOG(DEBUG) << "pre-range currentInteractionLength[cm]: " << currentInteractionLength/cm; + theNumberOfInteractionLengthLeft = 0; // to force this process + *condition = NotForced; // does nothing at this point + } else { // primary muon that didn't undergo DIS yet + *condition = NotForced; + + // In regions having multiple tiny volumes(e.g. ECC, SciFi), particle propagation is stopped + // numerous times at boundaries. This leads to sampling z for DIS too frequently, and due to + // competition with 'border crossing' leads to preference for small steps for the DIS. + // As a result, all DIS vertices are placed close to the set z_start. To avoid that, if the + // previous step was limited by boundary crossing (fTransportation), subtract its length + // from the proposed currentInteractionLength for DIS, and continue. + if (currentInteractionLength >= previousStepSize &&\ + track.GetStep()->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessType()==fTransportation ) { + currentInteractionLength -=previousStepSize; + return currentInteractionLength; + } + + // Calculate start/end positions along this muon track, and amount of material in between + // Use the track's position as start since position of DIS cannot be before current pos + // Convert to cm as needed in the MuDISGenerator and sndsw! + Double_t start[3] = {pos.x() / 10., pos.y() / 10., pos.z() / 10.}; + extrap_pos = linear_track_extrapolation(pos, mom, fZEnd); + extrap_pos /= 10.; + Double_t end[3] = {extrap_pos.x(), extrap_pos.y(), extrap_pos.z()}; + + Double_t current_dir[3] = {track.GetMomentumDirection().x(), track.GetMomentumDirection().y(), + track.GetMomentumDirection().z()}; + + // Use the MeanMaterialBudget code from the MuDISGenerator class + // Keep track of mean densities in the (likely) multiple steps + // from z_start to z_end. Update the maxmimum density encountered so far. + if (bparam == 0) { + memcpy(prev_start, start, sizeof(prev_start)); + } + if (bparam != 0) { + // Subtract density of previous range i.e. prev_pos to z_end + density_along_path -= (mparam[0] * mparam[4]); + // Recalculate the density of material traversed in previous step + bparam = disGen->MeanMaterialBudget(prev_start, start, mparam); + // Update with the actual mean density between prev_pos to current pos + density_along_path += (mparam[0] * mparam[4]); + // Keep track of the max density encountered so far: from z_start to current pos + prev_max_density = std::max(prev_max_density, mparam[7]); + memcpy(prev_start, start, sizeof(prev_start)); // update to current posistion + // Reset navigator + gGeoManager->CdTop(); // go to top volume + gGeoManager->InitTrack(start, current_dir); + gGeoManager->FindNode(); + g4Navigator->LocateGlobalPointAndSetup(pos, nullptr, false); + } + + // Calculate the density of material predicted to be traversed + bparam = disGen->MeanMaterialBudget(start, end, mparam); + + // Reset navigator + gGeoManager->CdTop(); // go to top volume + gGeoManager->InitTrack(start, current_dir); + gGeoManager->FindNode(); + g4Navigator->LocateGlobalPointAndSetup(pos, nullptr, false); + + // revert back to mm unit + for (int i = 0; i < 3; i++) { + start[i] *= 10.; + end[i] *= 10.; + } + + // Set the max density to the larger value between max density encountered from z_start to current pos + // and max density along the predicted path from current pos to z_end + max_density = std::max(prev_max_density, mparam[7]); + // Density along trajectory is the weight, g/cm^2 + density_along_path += (mparam[0] * mparam[4]); + fOutTracksWeight = density_along_path; + + LOG(DEBUG) << "in-range: density_along_path[g/cm2] " << density_along_path + << " max_density[g/cm3] " << max_density; + + G4double prob2int{}, zmu; + G4int count{}; + + // Loop over trajectory between start and end to pick an interaction point + while (prob2int < gRandom->Uniform(0., 1.)) { + zmu = gRandom->Uniform(start[2], end[2]); + extrap_pos = linear_track_extrapolation(pos, mom, zmu); + // get local material at this point + TGeoNode *node = gGeoManager->FindNode(extrap_pos.x() / 10., extrap_pos.y() / 10., zmu / 10.); + TGeoMaterial *material = 0; + if (node && !gGeoManager->IsOutside()){ + material = node->GetVolume()->GetMaterial(); + } + LOG(DEBUG) << "#trials= " << count << ", material name is " << material->GetName() + << ", density is " << material->GetDensity(); + // density relative to Prob largest density along this trajectory, i.e. use rho(Pt) + prob2int = material->GetDensity() / max_density; + if (prob2int > 1.) { + LOG(WARNING) << "MuDISGenerator: prob2int > Maximum density???? " << prob2int + << " maxrho along path in [z_start, z_end]: " << max_density + << " material: " << material->GetName() + << " full path: " << gGeoManager->GetPath() + << "\ncurrent pos " << pos + << "\nstart pos " << start[0] << " " << start[1] << " " << start[2] + << "\nend pos " << end[0] << " " << end[1] << " " << end[2]; + } + // Reset navigator again + gGeoManager->CdTop(); // go to top volume + gGeoManager->InitTrack(start, current_dir); + gGeoManager->FindNode(); + g4Navigator->LocateGlobalPointAndSetup(pos, nullptr, false); + + count += 1; + } + + LOG(DEBUG) << "prob2int " << prob2int << ", " << count; + LOG(DEBUG) << "put position[cm] " << extrap_pos/cm; + + // brilliant solution! Set the next to be the distance btw current pos and the chosen z for the DIS + currentInteractionLength = (extrap_pos - pos).mag(); + LOG(DEBUG) << "in-range currentInteractionLength[cm]: " << currentInteractionLength/cm; + theNumberOfInteractionLengthLeft = 0; // to force the interaction + } + +#ifdef G4VERBOSE + if (verboseLevel > 1) { + G4cout << "G4VDiscreteProcess::PostStepGetPhysicalInteractionLength "; + G4cout << "[ " << GetProcessName() << "]" << G4endl; + track.GetDynamicParticle()->DumpInfo(); + G4cout << " in Material " << track.GetMaterial()->GetName() << G4endl; + G4cout << "InteractionLength= " << currentInteractionLength / cm << "[cm] " << G4endl; + } +#endif + return currentInteractionLength; +} + +void G4MuonDISProcess::SimulateDIS(const G4Track &track) +{ + was_simulated = 1; + G4ThreeVector pos = track.GetPosition(); + + G4ThreeVector mom = track.GetMomentum(); + G4int pid = track.GetDefinition()->GetPDGEncoding(); + G4double w = track.GetWeight(); // it is the weight from the input file + G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable(); + + G4cout << "[DIS] Simulating muon DIS on "<< fNucleon <<" at: " + << pos/cm << "cm, E=" + << track.GetTotalEnergy() / GeV << " GeV\n"; + + double theta = mom.theta(); + double phi = mom.phi(); + double ctheta = cos(theta); + double stheta = sin(theta); + double cphi = cos(phi); + double sphi = sin(phi); + + TPythia6 *myPythia = new TPythia6(); + myPythia->SetMSEL(2); // minimum bias + // PARP(2) is the lowest c.m. energy for the event as a whole that the program + // will accept to simulate, default is 10 GeV. + // Set it to 2 GeV instead. + myPythia->SetPARP(2, 2); + + for (int kf : {211, 321, 130, 310, 3112, 3122, 3222, 3312, 3322, 3334}) { + int kc = myPythia->Pycomp(kf); + myPythia->SetMDCY(kc, 1, 0); // make them stable for pythia6 part + } + + int R = static_cast(std::time(nullptr) % 900000000); + myPythia->SetMRPY(1, R); // set the generator seed + + std::map mutype = {{-13, "gamma/mu+"}, {13, "gamma/mu-"}}; + + myPythia->SetMSTU(11, 11); // stop pythia printout during loop + myPythia->Initialize("FIXT", mutype[pid], fNucleon, mom.mag()/GeV); // using ’FIXT’, so 4th argument is beam mom in GeV + + myPythia->GenerateEvent(); + myPythia->Pyedit(2); // remove decayed particles + // Integrated cross section for subprocess ISUB, based on the statistics accumulated so far + // is directly accessible, see below. However it suffers from large uncertainty and Pythia6 + // manual reads that a proper xsec estimate can be obtained only by generating large statistics. + // This is why we use an external file with tabulated muon DIS xsec, where 1M events per reaction + // were used. + fCrossSection = (fXsecTables->at(pid))->Value(track.GetTotalEnergy()/GeV); // in mb + G4cout << "Total xsec for this 1 event: " << myPythia->GetPyint5()->XSEC[2][0] + << ", total xsec using 1M events: " << fCrossSection <<" mb\n"; + + // Kill the incoming muon after the interaction + // One can only propose the particle state change as the track is const! + aParticleChange.ProposeTrackStatus(fStopAndKill); + + // outgoing tracks + for (int itrk = 1; itrk <= myPythia->GetN(); ++itrk) { + int did = myPythia->GetK(itrk, 2); + double dpx, dpy, dpz; + std::tie(dpx, dpy, dpz) = + rotate(ctheta, stheta, cphi, sphi, myPythia->GetP(itrk, 1)*GeV, myPythia->GetP(itrk, 2)*GeV, myPythia->GetP(itrk, 3)*GeV); + double mass = myPythia->GetP(itrk, 5) * GeV; // mass in MeV: Pythia6 uses GeV/c2, so convert to Geant4's unit of MeV! + double p = sqrt(dpx * dpx + dpy * dpy + dpz * dpz); // mom in MeV + G4double kineticEnergy = sqrt(p * p + mass * mass) - mass; // Ekin in MeV + // Retrieve the particle definition using the PDG code + G4ParticleDefinition *particleDef = particleTable->FindParticle(did); + G4ThreeVector momentumDirection(dpx / p, dpy / p, dpz / p); + + G4DynamicParticle *dynamicParticle = new G4DynamicParticle(particleDef, momentumDirection, kineticEnergy); + G4ThreeVector position(pos.x(), pos.y(), pos.z()); // use position of mother + G4double time = track.GetGlobalTime(); // use age of mother + G4Track *newTrack = new G4Track(dynamicParticle, time, position); + + // Add to particle change + aParticleChange.AddSecondary(newTrack); + + newTrack->SetParentID(track.GetTrackID()); + newTrack->SetTouchableHandle(track.GetTouchableHandle()); // Important for geometry context + // weight = in_muon_weight * density along trajectory * cross_section + // On top, one needs to make two adjustments: + // 1. normalize material density to nucleon mass to get the number of 'targets'. + // We use proton_mass = 1.67x1E-24 grams. + // 2. convert xsec from mb to cm^2: 1mb = 1E-27 cm^2 + newTrack->SetWeight(w*fOutTracksWeight/1.67E-24*fCrossSection*1E-27); + } + //myPythia->Pystat(1); // extensive printout of parameters: xsec, BR, etc. + //myPythia->Pylist(1); // list table of properties for all particles + myPythia->SetMSTU(11, 6); // file number to which all program output is directed + + delete myPythia; +} + +std::tuple +G4MuonDISProcess::rotate(double ctheta, double stheta, double cphi, double sphi, double px, double py, double pz) +{ + // rotate around y-axis + double px1 = ctheta * px + stheta * pz; + double pzr = -stheta * px + ctheta * pz; + // rotate around z-axis + double pxr = cphi * px1 - sphi * py; + double pyr = sphi * px1 + cphi * py; + return std::make_tuple(pxr, pyr, pzr); +} + +G4ThreeVector G4MuonDISProcess::linear_track_extrapolation(G4ThreeVector pos, G4ThreeVector mom, G4double z_ref) +{ + double lam = (z_ref - pos.z()) / mom.z(); + return G4ThreeVector(pos.x() + lam * mom.x(), pos.y() + lam * mom.y(), z_ref); +} diff --git a/G4Processes/G4MuonDISProcess.hh b/G4Processes/G4MuonDISProcess.hh new file mode 100644 index 0000000000..cecac91f54 --- /dev/null +++ b/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 +#include +#include + +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 *x_range, std::vector *y_range, std::vector *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> 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 + 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> 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 From 36b336c0127d37175f0c38c89a18d0f9535024c9 Mon Sep 17 00:00:00 2001 From: siilieva Date: Tue, 15 Jul 2025 18:08:52 +0200 Subject: [PATCH 5/8] run_simSND: update logger settings Better control of severity level and make the log print in colour. --- shipLHC/run_simSND.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 60b78420d1..2b6a933027 100644 --- a/shipLHC/run_simSND.py +++ b/shipLHC/run_simSND.py @@ -76,6 +76,9 @@ # user hook userTask = False +# make log colourful +ROOT.FairLogger.GetLogger().SetColoredLog(True) + class MyTask(ROOT.FairTask): "user task" @@ -87,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" @@ -206,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) @@ -408,7 +414,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--------------------------------------------- From 52c44a1b57a7468d75d30962a81ed61c45a6e61c Mon Sep 17 00:00:00 2001 From: siilieva Date: Fri, 18 Jul 2025 14:40:09 +0200 Subject: [PATCH 6/8] MeanMaterialBudget: Unblock navigation along volume borders In rare cases, a track is moving right along the border of two volumes. This makes the geo navigation stuck. To unblock it, one now moves along the particle trajectory with a step of 1e-6cm, smaller than any position measurement precision in SND. This way, one allows the navigator to move away from any border at a reasonable computing cost. --- shipgen/MuDISGenerator.cxx | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index 211bfae8f9..437893ed3e 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/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]; From 2c9636c27b4cdd5816f9fb58580d54121aef5f5f Mon Sep 17 00:00:00 2001 From: siilieva Date: Wed, 20 Aug 2025 16:15:13 +0200 Subject: [PATCH 7/8] Fix the direct extraction of pythia xsec from the code When going from original p6 code(Fortran) to ROOT's TPyhia6 some arrays, like the XSEC one, have their indices swapped, also counting starts from 0 in TPythia6 case. With this fix one can read the P6 xsec directly from the code, no need to look into the produced txt file. Another part of the fix is to allow using P6 or P8 when creating the DIS xsec file using existing CL options. --- shipLHC/muonDis.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/shipLHC/muonDis.py b/shipLHC/muonDis.py index 25401882fe..6b235fb89b 100755 --- a/shipLHC/muonDis.py +++ b/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() From 590ce52a75fa014cc33d03fc7e30e77feb6a2489 Mon Sep 17 00:00:00 2001 From: siilieva Date: Wed, 3 Sep 2025 21:10:17 +0200 Subject: [PATCH 8/8] Allow setting Ekin cut when using FastMuon option This is essential when running muonDIS with activated emulsions. --- shipLHC/Floor.cxx | 6 +++++- shipLHC/run_simSND.py | 4 ++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/shipLHC/Floor.cxx b/shipLHC/Floor.cxx index ecc70bc024..d39d33bb81 100644 --- a/shipLHC/Floor.cxx +++ b/shipLHC/Floor.cxx @@ -109,7 +109,11 @@ void Floor::PreTrack(){ } TLorentzVector mom; gMC->TrackMomentum(mom); - if ( (mom.E()-mom.M() )StopTrack(); return; } diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 2b6a933027..e5d708fb12 100644 --- a/shipLHC/run_simSND.py +++ b/shipLHC/run_simSND.py @@ -309,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!!