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 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/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() diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 85d5e19778..e5d708fb12 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) @@ -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--------------------------------------------- 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]; 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); 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 - - - -