Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-18 07:54:43

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file SteppingAction.cc
0027 /// \brief Implementation of the SteppingAction class
0028 
0029 // This example is provided by the Geant4-DNA collaboration
0030 // Any report or published results obtained using the Geant4-DNA software
0031 // shall cite the following Geant4-DNA collaboration publications:
0032 // Med. Phys. 45 (2018) e722-e739
0033 // Phys. Med. 31 (2015) 861-874
0034 // Med. Phys. 37 (2010) 4692-4708
0035 // Int. J. Model. Simul. Sci. Comput. 1 (2010) 157–178
0036 //
0037 // The Geant4-DNA web site is available at http://geant4-dna.org
0038 //
0039 
0040 #include "SteppingAction.hh"
0041 #include "SteppingMessenger.hh"
0042 
0043 #include "DetectorConstruction.hh"
0044 #include "PrimaryGeneratorAction.hh"
0045 #include "RunAction.hh"
0046 
0047 #include "G4Alpha.hh"
0048 #include "G4AnalysisManager.hh"
0049 #include "G4DNAGenericIonsManager.hh"
0050 #include "G4Electron.hh"
0051 #include "G4Event.hh"
0052 #include "G4EventManager.hh"
0053 #include "G4Gamma.hh"
0054 #include "G4Proton.hh"
0055 #include "G4SteppingManager.hh"
0056 #include "G4SystemOfUnits.hh"
0057 
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059 
0060 SteppingAction::SteppingAction() : G4UserSteppingAction() 
0061 {
0062   fSteppingMessenger = new SteppingMessenger(this);
0063 }
0064 
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066 
0067 SteppingAction::~SteppingAction() 
0068 {
0069   delete fSteppingMessenger;
0070 }
0071 
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073 
0074 void SteppingAction::UserSteppingAction(const G4Step* step)
0075 {
0076   // Protection
0077   if (!step->GetPostStepPoint()) return;
0078   if (!step->GetPostStepPoint()->GetProcessDefinedStep()) return;
0079 
0080   // ***** FOR FIRST STEP RECORD ONLY   
0081   if (step->GetTrack()->GetCurrentStepNumber()>=1 && fKill)
0082     step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
0083   //
0084   
0085   //
0086   G4double flagParticle = -1.;
0087   G4double flagProcess = -1.;
0088   G4double x, y, z, xp, yp, zp;
0089 
0090   // 1) Particle identification
0091 
0092   // The following method avoids the usage of string comparison
0093 
0094   G4ParticleDefinition* partDef =
0095     step->GetTrack()->GetDynamicParticle()->GetDefinition();
0096 
0097   if (partDef == G4Gamma::GammaDefinition()) flagParticle = 0;
0098 
0099   if (partDef == G4Electron::ElectronDefinition()) flagParticle = 1;
0100 
0101   if (partDef == G4Proton::ProtonDefinition()) flagParticle = 2;
0102 
0103   if (partDef == G4Alpha::AlphaDefinition()) flagParticle = 4;
0104 
0105   G4DNAGenericIonsManager* instance;
0106   instance = G4DNAGenericIonsManager::Instance();
0107 
0108   // G4DNAGenericIonsManager usage example
0109   /*
0110   G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
0111   G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
0112   G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
0113   G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
0114   G4ParticleDefinition* heliumDef = instance->GetIon("helium");
0115   */
0116 
0117   if (partDef == instance->GetIon("hydrogen")) flagParticle = 3;
0118 
0119   if (partDef == instance->GetIon("alpha+")) flagParticle = 5;
0120 
0121   if (partDef == instance->GetIon("helium")) flagParticle = 6;
0122 
0123   // Heavier ions
0124   if (partDef->GetPDGCharge()>4) flagParticle = 7;
0125 
0126   // Alternative method (based on string comparison ) - not recommended -
0127   /*
0128   const G4String& particleName = step->GetTrack()->GetDynamicParticle()->
0129       GetDefinition()->GetParticleName();
0130 
0131   if (particleName == "gamma")         flagParticle = 0;
0132   else if (particleName == "e-")       flagParticle = 1;
0133   else if (particleName == "proton")   flagParticle = 2;
0134   else if (particleName == "hydrogen") flagParticle = 3;
0135   else if (particleName == "alpha")    flagParticle = 4;
0136   else if (particleName == "alpha+")   flagParticle = 5;
0137   else if (particleName == "helium")   flagParticle = 6;
0138   */
0139 
0140   // 2) Process identification
0141   //
0142   // Reminder
0143   // Process sub-types are listed in G4PhysicsListHelper.cc
0144   // or in Geant4-DNA process class implementation files (*.cc)
0145 
0146   G4StepPoint* preStep = step->GetPreStepPoint();
0147   G4StepPoint* postStep = step->GetPostStepPoint();
0148   G4int procID = postStep->GetProcessDefinedStep()->GetProcessSubType();
0149 
0150   const G4String& processName =
0151     postStep->GetProcessDefinedStep()->GetProcessName();
0152 
0153   // For gammas
0154   if (flagParticle == 0) {
0155     if (procID == 12)
0156       flagProcess = 81;
0157     else if (procID == 13)
0158       flagProcess = 82;
0159     else if (procID == 14)
0160       flagProcess = 83;
0161     else if (procID == 11)
0162       flagProcess = 84;
0163   }
0164 
0165   // For electrons
0166   else if (flagParticle == 1) {
0167     if (procID == 58)
0168       flagProcess = 10;
0169     else if (procID == 51)
0170       flagProcess = 11;
0171     else if (procID == 52)
0172       flagProcess = 12;
0173     else if (procID == 53)
0174       flagProcess = 13;
0175     else if (procID == 55)
0176       flagProcess = 14;
0177     else if (procID == 54)
0178       flagProcess = 15;
0179     else if (procID == 10)
0180       flagProcess = 110;
0181     else if (procID == 1)
0182       flagProcess = 120;
0183     else if (procID == 2)
0184       flagProcess = 130;
0185   }
0186 
0187   // For protons
0188   else if (flagParticle == 2) {
0189     if (procID == 51)
0190       flagProcess = 21;
0191     else if (procID == 52)
0192       flagProcess = 22;
0193     else if (procID == 53)
0194       flagProcess = 23;
0195     else if (procID == 56)
0196       flagProcess = 24;
0197     else if (procID == 10)
0198       flagProcess = 210;
0199     else if (procID == 1)
0200       flagProcess = 220;
0201     else if (procID == 2)
0202       flagProcess = 230;
0203     else if (procID == 8)
0204       flagProcess = 240;
0205   }
0206 
0207   // For hydrogen
0208   else if (flagParticle == 3) {
0209     if (procID == 51)
0210       flagProcess = 31;
0211     else if (procID == 52)
0212       flagProcess = 32;
0213     else if (procID == 53)
0214       flagProcess = 33;
0215     else if (procID == 57)
0216       flagProcess = 35;
0217   }
0218 
0219   // For alpha
0220   else if (flagParticle == 4) {
0221     if (procID == 51)
0222       flagProcess = 41;
0223     else if (procID == 52)
0224       flagProcess = 42;
0225     else if (procID == 53)
0226       flagProcess = 43;
0227     else if (procID == 56)
0228       flagProcess = 44;
0229     else if (procID == 10)
0230       flagProcess = 410;
0231     else if (procID == 1)
0232       flagProcess = 420;
0233     else if (procID == 2)
0234       flagProcess = 430;
0235     else if (procID == 8)
0236       flagProcess = 440;
0237   }
0238 
0239   // For alpha+
0240   else if (flagParticle == 5) {
0241     if (procID == 51)
0242       flagProcess = 51;
0243     else if (procID == 52)
0244       flagProcess = 52;
0245     else if (procID == 53)
0246       flagProcess = 53;
0247     else if (procID == 56)
0248       flagProcess = 54;
0249     else if (procID == 57)
0250       flagProcess = 55;
0251     else if (procID == 10)
0252       flagProcess = 510;
0253     else if (procID == 1)
0254       flagProcess = 520;
0255     else if (procID == 2)
0256       flagProcess = 530;
0257     else if (procID == 8)
0258       flagProcess = 540;
0259   }
0260 
0261   // For helium
0262   else if (flagParticle == 6) {
0263     if (procID == 51)
0264       flagProcess = 61;
0265     else if (procID == 52)
0266       flagProcess = 62;
0267     else if (procID == 53)
0268       flagProcess = 63;
0269     else if (procID == 57)
0270       flagProcess = 65;
0271   }
0272 
0273   else if (flagParticle == 7) {
0274     // For generic ions
0275     if (processName == "Capture")
0276       flagProcess = 1;
0277       // This process is used to kill ions below tracking cut
0278       // No subType and procID exists at the moment for this process,
0279       // so testing is performed using string comparison
0280     else if (procID == 210)
0281       flagProcess = 2;
0282       // This is the RadioactiveDecay process
0283     else if (procID == 53) // GenericIon_G4DNAIonisation
0284       flagProcess = 73;
0285     else if (procID == 10) // msc
0286       flagProcess = 710;
0287     else if (procID == 1) // CoulombScat
0288       flagProcess = 720;
0289     else if (procID == 2) // ionIoni
0290       flagProcess = 730;
0291     else if (procID == 8) // nuclearStopping
0292       flagProcess = 740;
0293   }
0294 
0295   // Alternative method (based on string comparison ) - not recommended -
0296   /*
0297   else if (processName=="e-_G4DNAElectronSolvation")    flagProcess =10;
0298   else if (processName=="e-_G4DNAElastic")              flagProcess =11;
0299   else if (processName=="e-_G4DNAExcitation")           flagProcess =12;
0300   else if (processName=="e-_G4DNAIonisation")           flagProcess =13;
0301   else if (processName=="e-_G4DNAAttachment")           flagProcess =14;
0302   else if (processName=="e-_G4DNAVibExcitation")        flagProcess =15;
0303 
0304   else if (processName=="proton_G4DNAElastic")          flagProcess =21;
0305   else if (processName=="proton_G4DNAExcitation")       flagProcess =22;
0306   else if (processName=="proton_G4DNAIonisation")       flagProcess =23;
0307   else if (processName=="proton_G4DNAChargeDecrease")   flagProcess =24;
0308 
0309   else if (processName=="hydrogen_G4DNAElastic")        flagProcess =31;
0310   else if (processName=="hydrogen_G4DNAExcitation")     flagProcess =32;
0311   else if (processName=="hydrogen_G4DNAIonisation")     flagProcess =33;
0312   else if (processName=="hydrogen_G4DNAChargeIncrease") flagProcess =35;
0313 
0314   else if (processName=="alpha_G4DNAElastic")           flagProcess =41;
0315   else if (processName=="alpha_G4DNAExcitation")        flagProcess =42;
0316   else if (processName=="alpha_G4DNAIonisation")        flagProcess =43;
0317   else if (processName=="alpha_G4DNAChargeDecrease")    flagProcess =44;
0318 
0319   else if (processName=="alpha+_G4DNAElastic")          flagProcess =51;
0320   else if (processName=="alpha+_G4DNAExcitation")       flagProcess =52;
0321   else if (processName=="alpha+_G4DNAIonisation")       flagProcess =53;
0322   else if (processName=="alpha+_G4DNAChargeDecrease")   flagProcess =54;
0323   else if (processName=="alpha+_G4DNAChargeIncrease")   flagProcess =55;
0324 
0325   else if (processName=="helium_G4DNAElastic")          flagProcess =61;
0326   else if (processName=="helium_G4DNAExcitation")       flagProcess =62;
0327   else if (processName=="helium_G4DNAIonisation")       flagProcess =63;
0328   else if (processName=="helium_G4DNAChargeIncrease")   flagProcess =65;
0329 
0330   else if (processName=="GenericIon_G4DNAIonisation")   flagProcess =73;
0331   */
0332 
0333   // 3) Fill ntuples
0334 
0335   if (processName != "Transportation") {
0336     x = preStep->GetPosition().x() / nanometer;
0337     y = preStep->GetPosition().y() / nanometer;
0338     z = preStep->GetPosition().z() / nanometer;
0339 
0340     xp = postStep->GetPosition().x() / nanometer;
0341     yp = postStep->GetPosition().y() / nanometer;
0342     zp = postStep->GetPosition().z() / nanometer;
0343 
0344     G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0345 
0346     analysisManager->FillNtupleDColumn(0, flagParticle);
0347     analysisManager->FillNtupleDColumn(1, flagProcess);
0348     analysisManager->FillNtupleDColumn(2, xp);
0349     analysisManager->FillNtupleDColumn(3, yp);
0350     analysisManager->FillNtupleDColumn(4, zp);
0351     analysisManager->FillNtupleDColumn(5, step->GetTotalEnergyDeposit() / eV);
0352 
0353     analysisManager->FillNtupleDColumn(
0354       6, std::sqrt((x - xp)*(x - xp) + (y - yp)*(y - yp) + (z - zp)*(z - zp)));
0355 
0356     analysisManager->FillNtupleDColumn(
0357       7, (preStep->GetKineticEnergy() - postStep->GetKineticEnergy()) / eV);
0358 
0359     analysisManager->FillNtupleDColumn(8, preStep->GetKineticEnergy() / eV);
0360 
0361     analysisManager->FillNtupleDColumn(9, preStep->GetMomentumDirection()
0362                                             * postStep->GetMomentumDirection());
0363 
0364     analysisManager->FillNtupleIColumn(
0365       10, G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID());
0366 
0367     analysisManager->FillNtupleIColumn(11, step->GetTrack()->GetTrackID());
0368 
0369     analysisManager->FillNtupleIColumn(12, step->GetTrack()->GetParentID());
0370 
0371     analysisManager->FillNtupleIColumn(13, step->GetTrack()->GetCurrentStepNumber());
0372 
0373     analysisManager->AddNtupleRow();
0374   }
0375 }