Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-11-04 09:29:51

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