microdosimetry/src/SteppingAction.cc

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