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