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