Merge pull request #113 from mnovak42/addGNuclear

Adding gamma-nuclear interaction.

G4HepEm/G4HepEm/include/G4HepEmTrackingManager.hh

@@ -11,6 +11,7 @@ class G4HepEmTLData;
 class G4SafetyHelper;
 class G4Step;
 class G4VProcess;
+class G4VParticleChange;
 class G4Region;
 
 #include <vector>
@@ -51,6 +52,14 @@ private:
   double StackSecondaries(G4HepEmTLData* aTLData, G4Track* aG4PrimaryTrack,
                           const G4VProcess* aG4CreatorProcess, int aG4IMC);
 
+  // Stacks secondaries created by Geant4 physics (if any) and returns with the
+  // energy deposit while stacking due to applying secondary production cuts
+  double StackG4Secondaries(G4VParticleChange* particleChange,
+                            G4Track* aG4PrimaryTrack,
+                            const G4VProcess* aG4CreatorProcess, int aG4IMC);
+
+  void InitNuclearProcesses(int particleID);
+
   // Checks if the particles has fast simulation maanger process attached and
   // stores in the local `fFastSimProcess` array (indexed by HepEm particle ID)
   void InitFastSimRelated(int particleID);
@@ -82,6 +91,9 @@ private:
   std::vector<G4HepEmNoProcess *> fGammaNoProcessVector;
   G4HepEmNoProcess *fTransportNoProcess;
 
+  // Pointers to the Gamma-nuclear process (if any)
+  G4VProcess* fGNucProcess;
+
   // Pointers to the fast simulation manager processes of the 3 particles if any
   // [0] e-; [1] e+; [2] gamma; nullptr: no fast sim manager process attached
   G4VProcess* fFastSimProcess[3];

G4HepEm/G4HepEm/src/G4HepEmTrackingManager.cc

@@ -32,6 +32,7 @@
 
 #include "G4VProcess.hh"
 #include "G4EmProcessSubType.hh"
+#include "G4HadronicProcessType.hh"
 #include "G4ProcessType.hh"
 #include "G4TransportationProcessType.hh"
 
@@ -41,6 +42,7 @@
 #include "G4Gamma.hh"
 #include "G4Positron.hh"
 
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4HepEmTrackingManager::G4HepEmTrackingManager() {
@@ -75,9 +77,16 @@ G4HepEmTrackingManager::G4HepEmTrackingManager() {
   fGammaNoProcessVector.push_back(
       new G4HepEmNoProcess("phot", G4ProcessType::fElectromagnetic,
                            G4EmProcessSubType::fPhotoElectricEffect));
+   fGammaNoProcessVector.push_back(
+       new G4HepEmNoProcess("photonNuclear", G4ProcessType::fHadronic,
+                            G4HadronicProcessType::fHadronInelastic));
+
   fTransportNoProcess = new G4HepEmNoProcess(
       "Transportation", G4ProcessType::fTransportation, TRANSPORTATION);
 
+  // Init the gamma-nuclear process
+  fGNucProcess = nullptr;
+
   // Init the fast sim mamanger process ptrs of the 3 particles
   fFastSimProcess[0] = nullptr;
   fFastSimProcess[1] = nullptr;
@@ -119,6 +128,8 @@ void G4HepEmTrackingManager::BuildPhysicsTable(const G4ParticleDefinition &part)
   } else if (&part == G4Gamma::Definition()) {
     int particleID = 2;
     fRunManager->Initialize(fRandomEngine, particleID);
+    // Find the gamma-nuclear process if has been attached
+    InitNuclearProcesses(particleID);
     // Find the fast simulation manager process for gamma (if has been attached)
     InitFastSimRelated(particleID);
   } else {
@@ -751,6 +762,10 @@ void G4HepEmTrackingManager::TrackGamma(G4Track *aTrack) {
   if (fFastSimProc != nullptr) {
     fFastSimProc->StartTracking(aTrack);
   }
+
+  if (fGNucProcess != nullptr) {
+    fGNucProcess->StartTracking(aTrack);
+  }
   // === StartTracking ===
 
   while (aTrack->GetTrackStatus() == fAlive) {
@@ -850,25 +865,67 @@ void G4HepEmTrackingManager::TrackGamma(G4Track *aTrack) {
         thePrimaryTrack->SetGStepLength(aTrack->GetStepLength());
         G4HepEmGammaManager::UpdateNumIALeft(thePrimaryTrack);
       } else {
-        // (NOTE: Ekin, MC-index, step-length, onBoundary have all set)
-        G4HepEmGammaManager::Perform(theHepEmData, theHepEmPars, theTLData);
-        proc = fGammaNoProcessVector[thePrimaryTrack->GetWinnerProcessIndex()];
-
-        // energy, e-depo, momentum direction and status
-        const double ekin = thePrimaryTrack->GetEKin();
-        double edep = thePrimaryTrack->GetEnergyDeposit();
-        postStepPoint.SetKineticEnergy(ekin);
-        if (ekin <= 0.0) {
-          aTrack->SetTrackStatus(fStopAndKill);
-        }
-        const double *pdir = thePrimaryTrack->GetDirection();
-        postStepPoint.SetMomentumDirection( G4ThreeVector(pdir[0], pdir[1], pdir[2]) );
-
-        step.UpdateTrack();
+        double edep = 0.0;
+        // NOTE: gamma-nuclear interaction needs to be done here while others in HepEm
+        const int iDProc = thePrimaryTrack->GetWinnerProcessIndex();
+        if (iDProc != 3) {
+          // Conversion, Compton or photoelectric --> use HepEm for the interaction
+          // (NOTE: Ekin, MC-index, step-length, onBoundary have all set)
+          G4HepEmGammaManager::Perform(theHepEmData, theHepEmPars, theTLData);
+          // energy, e-depo, momentum direction and status
+          const double ekin = thePrimaryTrack->GetEKin();
+          edep = thePrimaryTrack->GetEnergyDeposit();
+          postStepPoint.SetKineticEnergy(ekin);
+          if (ekin <= 0.0) {
+            aTrack->SetTrackStatus(fStopAndKill);
+          }
+          const double *pdir = thePrimaryTrack->GetDirection();
+          postStepPoint.SetMomentumDirection( G4ThreeVector(pdir[0], pdir[1], pdir[2]) );
 
-        // Stack secondaries created by the HepEm physics above
-        edep += StackSecondaries(theTLData, aTrack, proc, g4IMC);
+          step.UpdateTrack();
 
+          // Stack secondaries created by the HepEm physics above
+          edep += StackSecondaries(theTLData, aTrack, proc, g4IMC);
+
+        } else {
+          // Gamma-nuclear: --> use Geant4 for the interaction:
+          // NOTE: it's destructive i.e. stopps and kills the gammma when the
+          //    interaction happens so we do not use anymore the primary HepEm
+          //    track thus we should not update the number of interaction length
+          //    left and the others done in G4HepEmGammaManager::Perform before
+          //    the interaction. HOWEVER, the interaction might not happen at
+          //    the end so we do what we need in that case anyway for clarity:
+          // - update the number of interaction length left for all process
+          // - reset the number of interaction lenght left for the winner
+          //   process, i.e. the one that will be invoked (the gamnma nuclear)
+          // - reset the energy deposit to zero
+          G4HepEmGammaManager::UpdateNumIALeft(thePrimaryTrack);
+          thePrimaryTrack->SetNumIALeft(-1.0, iDProc);
+          thePrimaryTrack->SetEnergyDeposit(0.0);
+          // Invoke the gamma-nuclear interaction using the Geant4 process
+          G4VParticleChange* particleChangeGNuc = nullptr;
+          if (fGNucProcess != nullptr) {
+            // call to set some fields of the process like material, energy etc...
+            G4ForceCondition forceCondition;
+            fGNucProcess->PostStepGetPhysicalInteractionLength(*aTrack, 0.0, &forceCondition);
+            //
+            postStepPoint.SetStepStatus(fPostStepDoItProc);
+            particleChangeGNuc = fGNucProcess->PostStepDoIt(*aTrack, step);
+            // update the track and stack according to the result of the interaction
+            particleChangeGNuc->UpdateStepForPostStep(&step);
+            step.UpdateTrack();
+            aTrack->SetTrackStatus(particleChangeGNuc->GetTrackStatus());
+            // need to add secondaries to the secondary vector of the current track
+            // NOTE: as we use Geant4, we should care only those changes that are
+            //   not included in the above update step and track, i.e. the energy
+            //   deposited due to applying the cut when stacking the secondaries
+            edep = StackG4Secondaries(particleChangeGNuc, aTrack, fGammaNoProcessVector[iDProc], g4IMC);
+            // done: clear the particle change
+            particleChangeGNuc->Clear();
+          }
+        }
+        // Set process defined setp and add edep to the step
+        proc = fGammaNoProcessVector[iDProc];
         step.AddTotalEnergyDeposit(edep);
       } // END if NOT onBoundary
 
@@ -934,8 +991,9 @@ void G4HepEmTrackingManager::HandOverOneTrack(G4Track *aTrack) {
   delete aTrack;
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+// Helper that can be used to stack secondary e-/e+ and gamma i.e. everything
+// that HepEm physics can produce
 double G4HepEmTrackingManager::StackSecondaries(G4HepEmTLData* aTLData, G4Track* aG4PrimaryTrack, const G4VProcess* aG4CreatorProcess, int aG4IMC) {
   const int numSecElectron = aTLData->GetNumSecondaryElectronTrack();
   const int numSecGamma    = aTLData->GetNumSecondaryGammaTrack();
@@ -1014,7 +1072,74 @@ double G4HepEmTrackingManager::StackSecondaries(G4HepEmTLData* aTLData, G4Track*
   return edep;
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Helper that can be used to stack secondary e-/e+ and gamma i.e. everything
+// that HepEm physics can produce
+double G4HepEmTrackingManager::StackG4Secondaries(G4VParticleChange* particleChange, G4Track* aG4PrimaryTrack, const G4VProcess* aG4CreatorProcess, int aG4IMC) {
+  const int numSecondaries = particleChange->GetNumberOfSecondaries();
+  // return early if there are no secondaries created by the physics interaction
+  double edep = 0.0;
+  if (numSecondaries == 0) {
+    return edep;
+  }
+
+  G4Step&        step           = *fStep;
+  G4TrackVector& secondaries    = *step.GetfSecondary();
+  G4StepPoint&   postStepPoint  = *step.GetPostStepPoint();
+
+  const G4ThreeVector&     theG4PostStepPointPosition = postStepPoint.GetPosition();
+  const G4double           theG4PostStepGlobalTime    = postStepPoint.GetGlobalTime();
+  const G4TouchableHandle& theG4TouchableHandle       = aG4PrimaryTrack->GetTouchableHandle();
+  const double             theG4ParentTrackWeight     = aG4PrimaryTrack->GetWeight();
+  const int                theG4ParentTrackID         = aG4PrimaryTrack->GetTrackID();
+
+  for (int isec=0; isec<particleChange->GetNumberOfSecondaries(); ++isec) {
+    G4Track *secTrack = particleChange->GetSecondary(isec);
+    double   secEKin  = secTrack->GetKineticEnergy();
+    const G4ParticleDefinition* secPartDef = secTrack->GetParticleDefinition();
+//    std::cout << "applycutS = " << applyCuts << " Ekin = " << secEKin  << " imc = " << aG4IMC << " par = " << secPartDef->GetParticleName()<<std::endl;
+
+    if (applyCuts) {
+      if (secPartDef == G4Gamma::Definition() && secEKin < (*theCutsGamma)[aG4IMC]) {
+        edep += secEKin;
+        continue;
+      } else if (secPartDef == G4Electron::Definition() && secEKin < (*theCutsElectron)[aG4IMC]) {
+        edep += secEKin;
+        continue;
+      } else if (secPartDef == G4Positron::Definition() && CLHEP::electron_mass_c2 < (*theCutsGamma)[aG4IMC]
+                 && secEKin < (*theCutsPositron)[aG4IMC]) {
+        edep += secEKin + 2 * CLHEP::electron_mass_c2;
+        continue;
+      }
+    }
+    secTrack->SetParentID(theG4ParentTrackID);
+    secTrack->SetCreatorProcess(aG4CreatorProcess);
+    secTrack->SetTouchableHandle(theG4TouchableHandle);
+    secTrack->SetWeight(theG4ParentTrackWeight);
+    secondaries.push_back(secTrack);
+  }
+
+  return edep;
+}
+
+
+// Try to get the nuclear process pointer from the process manager of the particle
+void G4HepEmTrackingManager::InitNuclearProcesses(int particleID) {
+  if (particleID == 2) {
+    G4ParticleDefinition* partDef = G4Gamma::Definition();
+    const G4ProcessVector* processVector = partDef->GetProcessManager()->GetProcessList();
+    for (std::size_t ip=0; ip<processVector->entries(); ip++) {
+      if( (*processVector)[ip]->GetProcessName()=="photonNuclear") {
+        fGNucProcess = (*processVector)[ip];
+        // make sure the process is initialised (element selectors needs to be built)
+        fGNucProcess->PreparePhysicsTable(*partDef);
+        fGNucProcess->BuildPhysicsTable(*partDef);
+        break;
+      }
+    }
+  }
+}
+
 
 // Try to get the fast sim mamanger process for the particle (e- 0, e+ 1, gamma 2)
 void G4HepEmTrackingManager::InitFastSimRelated(int particleID) {

G4HepEm/G4HepEmData/include/G4HepEmGammaData.hh

@@ -28,9 +28,17 @@ struct G4HepEmGammaData {
   double        fCompEILDelta = 0.0;       // =  3.040061373322763; // 1./[log(emax/emin)/84]
   double*       fCompEnergyGrid = nullptr;     // [fCompEnergyGridSize]
 
-  // the macroscopic cross sections for all materials and for [conversion,compton]
+//// === gamma nuclear related data. 146 bins form  2mc^2 - 100 TeV
+  const int     fGNucEnergyGridSize = 147;
+  double        fGNucLogMinEkin = 0.0;     // =  0.021759358706830;  // log(2mc^2)
+  double        fGNucEILDelta = 0.0;       // =  7.935247775833226;  // 1./[log(emax/emin)/146]
+  double*       fGNucEnergyGrid = nullptr;     // [fGNucEnergyGridSize]
+
+
+  // the macroscopic cross sections for all materials and for [conversion,compton,gamma-nuclear]
   // at each material
-  double*       fConvCompMacXsecData = nullptr;   // [#materials*2*(fConvEnergyGridSize+fCompEnergyGridSize)]
+  double*       fConvCompGNucMacXsecData = nullptr;   // [#materials*2*(fConvEnergyGridSize+fCompEnergyGridSize+fGNucEnergyGridSize)]
+
 
 //// === element selector for conversion (note: KN compton interaction do not know anything about Z)
   int           fElemSelectorConvEgridSize = 0;

G4HepEm/G4HepEmData/src/G4HepEmGammaData.cc

@@ -17,7 +17,8 @@ void FreeGammaData (struct G4HepEmGammaData** theGammaData)  {
   if (*theGammaData != nullptr) {
     delete[] (*theGammaData)->fConvEnergyGrid ;
     delete[] (*theGammaData)->fCompEnergyGrid;
-    delete[] (*theGammaData)->fConvCompMacXsecData;
+    delete[] (*theGammaData)->fGNucEnergyGrid;
+    delete[] (*theGammaData)->fConvCompGNucMacXsecData;
     delete[] (*theGammaData)->fElemSelectorConvStartIndexPerMat;
     delete[] (*theGammaData)->fElemSelectorConvEgrid;
     delete[] (*theGammaData)->fElemSelectorConvData;
@@ -56,10 +57,15 @@ void CopyGammaDataToDevice(struct G4HepEmGammaData* onHOST, struct G4HepEmGammaD
   // allocate memory on _d for the Compton energy grid and copy them form _h
   gpuErrchk ( cudaMalloc ( &(gmDataHTo_d->fCompEnergyGrid), sizeof( double ) * numCompData ) );
   gpuErrchk ( cudaMemcpy (   gmDataHTo_d->fCompEnergyGrid,  onHOST->fCompEnergyGrid, sizeof( double ) * numCompData, cudaMemcpyHostToDevice ) );
+  // -- go for the gamma-nuclear related data
+  int numGNucData = onHOST->fGNucEnergyGridSize;
+  // allocate memory on _d for the gamma-nuclear energy grid and copy them form _h
+  gpuErrchk ( cudaMalloc ( &(gmDataHTo_d->fGNucEnergyGrid), sizeof( double ) * numGNucData ) );
+  gpuErrchk ( cudaMemcpy (   gmDataHTo_d->fGNucEnergyGrid,  onHOST->fGNucEnergyGrid, sizeof( double ) * numGNucData, cudaMemcpyHostToDevice ) );
   // allocate memory on _d for the conversion and Compton macroscopic x-section data and copy them form _h
-  int numConvCompData = numHepEmMat*2*(numConvData+numCompData);
-  gpuErrchk ( cudaMalloc ( &(gmDataHTo_d->fConvCompMacXsecData), sizeof( double ) * numConvCompData ) );
-  gpuErrchk ( cudaMemcpy (   gmDataHTo_d->fConvCompMacXsecData,  onHOST->fConvCompMacXsecData, sizeof( double ) * numConvCompData, cudaMemcpyHostToDevice ) );
+  int numConvCompGNucData = numHepEmMat*2*(numConvData+numCompData+numGNucData);
+  gpuErrchk ( cudaMalloc ( &(gmDataHTo_d->fConvCompGNucMacXsecData), sizeof( double ) * numConvCompGNucData ) );
+  gpuErrchk ( cudaMemcpy (   gmDataHTo_d->fConvCompGNucMacXsecData,  onHOST->fConvCompGNucMacXsecData, sizeof( double ) * numConvCompGNucData, cudaMemcpyHostToDevice ) );
   //
   // -- go for the conversion element selector related data
   int numElSelE   = onHOST->fElemSelectorConvEgridSize;
@@ -91,10 +97,11 @@ void FreeGammaDataOnDevice(struct G4HepEmGammaData** onDEVICE) {
     // side dynamically allocated memories
     struct G4HepEmGammaData* onHostTo_d = new G4HepEmGammaData;
     gpuErrchk ( cudaMemcpy( onHostTo_d, *onDEVICE, sizeof( struct G4HepEmGammaData ), cudaMemcpyDeviceToHost ) );
-    // conversion and Compton macroscopic x-section related data
+    // conversion, Compton and gamma-nuclear macroscopic x-section related data
     cudaFree( onHostTo_d->fConvEnergyGrid );
     cudaFree( onHostTo_d->fCompEnergyGrid );
-    cudaFree( onHostTo_d->fConvCompMacXsecData );
+    cudaFree( onHostTo_d->fGNucEnergyGrid );
+    cudaFree( onHostTo_d->fConvCompGNucMacXsecData );
     // conversion element selector related data
     cudaFree( onHostTo_d->fElemSelectorConvStartIndexPerMat );
     cudaFree( onHostTo_d->fElemSelectorConvEgrid );

G4HepEm/G4HepEmDataJsonIO/src/G4HepEmDataJsonIOImpl.hh

@@ -737,11 +737,19 @@ namespace nlohmann
         j["fCompEnergyGrid"] =
           make_span(d->fCompEnergyGridSize, d->fCompEnergyGrid);
 
+        //// === compton related data. 84 bins (7 per decades) from 100 eV - 100
+        /// TeV
+        j["fGNucLogMinEkin"] = d->fGNucLogMinEkin;
+        j["fGNucEILDelta"]   = d->fGNucEILDelta;
+        j["fGNucEnergyGrid"] =
+          make_span(d->fGNucEnergyGridSize, d->fGNucEnergyGrid);
+
+
         const int macXsecDataSize =
           d->fNumMaterials * 2 *
-          (d->fConvEnergyGridSize + d->fCompEnergyGridSize);
-        j["fConvCompMacXsecData"] =
-          make_span(macXsecDataSize, d->fConvCompMacXsecData);
+          (d->fConvEnergyGridSize + d->fCompEnergyGridSize + d->fGNucEnergyGridSize);
+        j["fConvCompGNucMacXsecData"] =
+          make_span(macXsecDataSize, d->fConvCompGNucMacXsecData);
 
         //// === element selector for conversion (note: KN compton interaction
         /// do not know anything about Z)
@@ -787,12 +795,20 @@ namespace nlohmann
           j.at("fCompEnergyGrid").get<dynamic_array<double>>();
         d->fCompEnergyGrid = tmpCompEnergyGrid.data;
 
+        j.at("fGNucLogMinEkin").get_to(d->fGNucLogMinEkin);
+        j.at("fGNucEILDelta").get_to(d->fGNucEILDelta);
+        // Get the array but ignore the size (fGNucEnergyGridSize) as this is a
+        // const (at time of writing)
+        auto tmpGNucEnergyGrid =
+          j.at("fGNucEnergyGrid").get<dynamic_array<double>>();
+        d->fGNucEnergyGrid = tmpGNucEnergyGrid.data;
+
         // We don't store the size of the following array, rather should
         // validate that it is expected size: d->fNumMaterials * 2 *
-        // (d->fConvEnergyGridSize + d->fCompEnergyGridSize);
-        auto tmpConvCompXsecData =
-          j.at("fConvCompMacXsecData").get<dynamic_array<double>>();
-        d->fConvCompMacXsecData = tmpConvCompXsecData.data;
+        // (d->fConvEnergyGridSize + d->fCompEnergyGridSize + d->fGNucEnergyGridSize);
+        auto tmpConvCompGNucXsecData =
+          j.at("fConvCompGNucMacXsecData").get<dynamic_array<double>>();
+        d->fConvCompGNucMacXsecData = tmpConvCompGNucXsecData.data;
 
         j.at("fElemSelectorConvLogMinEkin")
           .get_to(d->fElemSelectorConvLogMinEkin);

G4HepEm/G4HepEmInit/include/G4HepEmGammaTableBuilder.hh

@@ -8,6 +8,7 @@
 //class G4ParticleDefinition;
 class G4PairProductionRelModel;
 class G4KleinNishinaCompton;
+class G4CrossSectionDataStore;
 
 struct G4HepEmData;
 struct G4HepEmParameters;
@@ -16,7 +17,7 @@ struct G4HepEmParameters;
 
 // Should receive pointers to G4 models that are already initialised
 void BuildLambdaTables(G4PairProductionRelModel* ppModel, G4KleinNishinaCompton* knModel,
-                     struct G4HepEmData* hepEmData);
+                       G4CrossSectionDataStore* hadGNucXSDataStore, struct G4HepEmData* hepEmData);
 
 void BuildElementSelectorTables(G4PairProductionRelModel* ppModel, struct G4HepEmData* hepEmData);