dsbandrepair/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 //
0027 /// \file SteppingAction.cc
0028 /// \brief Implementation of the SteppingAction class
0029 /// \file SteppingAction.cc
0030 /// \brief Implementation of the SteppingAction class
0031
0032 #include "SteppingAction.hh"
0033 #include "Analysis.hh"
0034 #include "G4SteppingManager.hh"
0035 #include "G4VTouchable.hh"
0036 #include "G4VPhysicalVolume.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4LogicalVolumeStore.hh"
0039 #include "G4SteppingManager.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4Track.hh"
0042 #include "G4RunManager.hh"
0043 #include "G4Proton.hh"
0044 #include "G4Electron.hh"
0045 #include "G4Alpha.hh"
0046 #include "G4DNAGenericIonsManager.hh"
0047
0048 #ifdef USE_MPI
0049 #include "G4MPImanager.hh"
0050 #endif
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0053
0054 SteppingAction::SteppingAction(EventAction* pEvent)
0055 {
0056     fEventAction = pEvent;
0057 }
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0060
0061 void SteppingAction::UserSteppingAction(const G4Step* step)
0062 {
0063     SetupFlags(step);
0064
0065     SetupVoxelCopyNumber(step);
0066
0067     if(fFlagVolume>0)
0068     {
0069         fEventAction->AddEdep(step->GetTotalEnergyDeposit()/eV);
0070     }
0071     //G4cout<<step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName()<<" ; "<<step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName()<<G4endl;
0072     // Check we are in a DNA volume
0073     if(fFlagVolume == 1 // d1
0074             || fFlagVolume == 11 // p1
0075             || fFlagVolume == 2 // d2
0076             || fFlagVolume == 22 // p2
0077             || fFlagVolume == 3 // cyto
0078             || fFlagVolume == 4 // gua
0079             || fFlagVolume == 5 // thy
0080             || fFlagVolume == 6 // ade
0081             || fFlagVolume == 7 // d1_w
0082             || fFlagVolume == 71 // p1_w
0083             || fFlagVolume == 8 // d2_w
0084             || fFlagVolume == 81 // p2_w
0085             || fFlagVolume == 9 // ade_w
0086             || fFlagVolume == 10 // gua_w
0087             || fFlagVolume == 13 // cyto_w
0088             || fFlagVolume == 12) // thy_w
0089     {
0090         // *****************************************************
0091         // Saving physical stage informations
0092         // *****************************************************
0093         if(step->GetPostStepPoint()->GetTouchable()->GetVolume() ) // avoid asking non existing information in postStep
0094         {
0095             G4double x=step->GetPreStepPoint()->GetPosition().x()/nanometer;
0096             G4double y=step->GetPreStepPoint()->GetPosition().y()/nanometer;
0097             G4double z=step->GetPreStepPoint()->GetPosition().z()/nanometer;
0098
0099             G4double dE = step->GetTotalEnergyDeposit()/eV;
0100             G4double copyNo = G4double (step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetCopyNo() );
0101
0102             G4int eventId = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
0103 #ifdef USE_MPI
0104             auto g4MPI = G4MPImanager::GetManager();
0105             if (g4MPI->IsSlave()) { // update eventID only for slave, cause rank_master=0
0106                 G4int rank = g4MPI->GetRank();
0107                 eventId += g4MPI->GetEventsInMaster() + (rank-1)*g4MPI->GetEventsInSlave();
0108             }
0109 #endif
0110             // ***************************
0111             // put to vector
0112             // ***************************
0113             InfoInPhysStage aInfo;
0114             aInfo.fFlagParticle = fFlagParticle;
0115             aInfo.fFlagParentID = fFlagParentID;
0116             aInfo.fFlagProcess = fFlagProcess;
0117             aInfo.fX = x;
0118             aInfo.fY = y;
0119             aInfo.fZ = z;
0120             aInfo.fEdep = dE;
0121             aInfo.fEventNumber = eventId;
0122             aInfo.fVolumeName = fFlagVolume;
0123             aInfo.fCopyNumber = copyNo;
0124             aInfo.fLastMetVoxelCopyNum = fLastMetVoxelCopyNumber;
0125             Analysis::GetAnalysis()->AddInfoInPhysStage(aInfo);
0126         }
0127     }
0128 }
0129
0130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0131
0132 void SteppingAction::SetupFlags(const G4Step* step)
0133 {
0134     fFlagParticle=0;
0135     fFlagProcess=0;
0136     fFlagParentID=0;
0137     fFlagVolume=0;
0138
0139     fFlagParentID = step->GetTrack()->GetParentID();
0140
0141     SetupParticleAndProcessFlags(step);
0142
0143     fFlagVolume = SetupVolumeFlag(step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName() );
0144 }
0145
0146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0147
0148 void SteppingAction::SetupParticleAndProcessFlags(const G4Step *step)
0149 {
0150     fFlagParticle = 0;
0151     fFlagProcess = 0;
0152     auto partDef = step->GetTrack()->GetDynamicParticle()->GetDefinition();
0153     auto* instance = G4DNAGenericIonsManager::Instance();
0154
0155     if (partDef == G4Electron::ElectronDefinition())       fFlagParticle = 1;
0156     if (partDef == G4Proton::ProtonDefinition())   fFlagParticle = 2;
0157     if (partDef == instance->GetIon("hydrogen")) fFlagParticle = 3;
0158     if (partDef == G4Alpha::AlphaDefinition())    fFlagParticle = 4;
0159     if (partDef == instance->GetIon("alpha+"))   fFlagParticle = 5;
0160     if (partDef == instance->GetIon("helium"))   fFlagParticle = 6;
0161
0162     G4int procSubtype = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessSubType();
0163
0164
0165     if (fFlagParticle == 1) { // e-
0166         if (procSubtype == 58) fFlagProcess = 10; //"e-_G4DNAElectronSolvation"
0167         if (procSubtype == 51) fFlagProcess = 11; //"e-_G4DNAElastic"
0168         if (procSubtype == 52) fFlagProcess = 12; //"e-_G4DNAExcitation"
0169         if (procSubtype == 53) fFlagProcess = 13; //"e-_G4DNAIonisation"
0170         if (procSubtype == 55) fFlagProcess = 14; //"e-_G4DNAAttachment"
0171         if (procSubtype == 54) fFlagProcess = 15; //"e-_G4DNAVibExcitation"
0172     }
0173
0174     if (fFlagParticle == 2) { // Proton
0175         if (procSubtype == 52) fFlagProcess = 17; //"proton_G4DNAExcitation"
0176         if (procSubtype == 53) fFlagProcess = 18; //"proton_G4DNAIonisation"
0177         if (procSubtype == 56) fFlagProcess = 19; //"proton_G4DNAChargeDecrease"
0178     }
0179
0180     if (fFlagParticle == 3) { // hydrogen
0181         if (procSubtype == 52) fFlagProcess = 20; //"hydrogen_G4DNAExcitation"
0182         if (procSubtype == 53) fFlagProcess = 21; //"hydrogen_G4DNAIonisation"
0183         if (procSubtype == 57) fFlagProcess = 22; //"hydrogen_G4DNAChargeIncrease"
0184     }
0185
0186     if (fFlagParticle == 4) { // alpha
0187         if (procSubtype == 52) fFlagProcess = 23; //"alpha_G4DNAExcitation"
0188         if (procSubtype == 53) fFlagProcess = 24; //"alpha_G4DNAIonisation"
0189         if (procSubtype == 56) fFlagProcess = 25; //"alpha_G4DNAChargeDecrease"
0190     }
0191
0192     if (fFlagParticle == 5) { // alpha+
0193         if (procSubtype == 52) fFlagProcess = 26; //"alpha+_G4DNAExcitation"
0194         if (procSubtype == 53) fFlagProcess = 27; //"alpha+_G4DNAIonisation"
0195         if (procSubtype == 56) fFlagProcess = 28; //"alpha+_G4DNAChargeDecrease"
0196         if (procSubtype == 57) fFlagProcess = 29; //"alpha+_G4DNAChargeIncrease"
0197     }
0198
0199     if (fFlagParticle == 6) { // helium
0200         if (procSubtype == 52) fFlagProcess = 30; //"helium__G4DNAExcitation"
0201         if (procSubtype == 53) fFlagProcess = 31; //"helium__G4DNAIonisation"
0202         if (procSubtype == 57) fFlagProcess = 32; //"helium__G4DNAChargeIncrease"
0203     }
0204
0205 }
0206
0207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0208
0209 G4int SteppingAction::SetupVolumeFlag(const G4String& volumeName)
0210 {
0211     G4int flagVolume(-1);
0212
0213     if(volumeName=="deoxyribose1_phys") flagVolume = 1;
0214     else if(volumeName=="phosphate1_phys") flagVolume = 11;
0215     else if(volumeName=="deoxyribose2_phys") flagVolume = 2;
0216     else if(volumeName=="phosphate2_phys") flagVolume = 22;
0217     else if(volumeName=="base_cytosine_phys") flagVolume = 3;
0218     else if(volumeName=="base_guanine_phys") flagVolume = 4;
0219     else if(volumeName=="base_thymine_phys") flagVolume = 5;
0220     else if(volumeName=="base_adenine_phys") flagVolume = 6;
0221     else if(volumeName=="deoxyribose1_water_phys") flagVolume = 7;
0222     else if(volumeName=="phosphate1_water_phys") flagVolume = 71;
0223     else if(volumeName=="deoxyribose2_water_phys") flagVolume = 8;
0224     else if(volumeName=="phosphate2_water_phys") flagVolume = 81;
0225     else if(volumeName=="base_adenine_water_phys") flagVolume = 9;
0226     else if(volumeName=="base_guanine_water_phys") flagVolume = 10;
0227     else if(volumeName=="base_cytosine_water_phys") flagVolume = 13;
0228     else if(volumeName=="base_thymine_water_phys") flagVolume = 12;
0229     else if(volumeName=="fiber") flagVolume = 110;
0230     else if(volumeName=="voxelStraight" || volumeName=="VoxelStraight") flagVolume = 161;
0231     else if(volumeName=="voxelRight" || volumeName=="VoxelRight") flagVolume = 162;
0232     else if(volumeName=="voxelLeft" || volumeName=="VoxelLeft") flagVolume = 163;
0233     else if(volumeName=="voxelUp" || volumeName=="VoxelUp") flagVolume = 164;
0234     else if(volumeName=="voxelDown" || volumeName=="VoxelDown") flagVolume = 165;
0235     else if(volumeName=="voxelStraight2" || volumeName=="VoxelStraight2") flagVolume = 261;
0236     else if(volumeName=="voxelRight2" || volumeName=="VoxelRight2") flagVolume = 262;
0237     else if(volumeName=="voxelLeft2" || volumeName=="VoxelLeft2") flagVolume = 263;
0238     else if(volumeName=="voxelUp2" || volumeName=="VoxelUp2") flagVolume = 264;
0239     else if(volumeName=="voxelDown2" || volumeName=="VoxelDown2") flagVolume = 265;
0240     else if(volumeName=="physWorld") flagVolume = -1.;
0241     else if(volumeName=="wrapper") flagVolume = 130;
0242     else if(volumeName=="histone_phys") flagVolume = 140;
0243     else if(volumeName=="nucleus_pl") flagVolume = 200;
0244     return flagVolume;
0245 }
0246
0247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0248
0249 void SteppingAction::SetupVoxelCopyNumber(const G4Step* step)
0250 {
0251     // Each time we will do a step in a voxel, we will setup the lastMetVoxelCpNum flag.
0252     // This way, this number will always be the one of the last voxel met by the particle.
0253     // Therefore, in a DNA volume, the number will be the one of the container voxel.
0254
0255     if(step->GetPostStepPoint()->GetTouchable()->GetVolume() ) // avoid asking non existing information in postStep
0256     {
0257         const G4String& volPost = step->GetPostStepPoint()->GetTouchable()->GetVolume()->GetName();
0258
0259         // If we enter a voxel or stay in it
0260         if(    volPost == "VoxelStraight" || volPost=="voxelStraight"
0261             || volPost == "VoxelRight" || volPost=="voxelRight"
0262             || volPost == "VoxelLeft" || volPost=="voxelLeft"
0263             || volPost == "VoxelUp" || volPost=="voxelUp"
0264             || volPost == "VoxelDown" || volPost=="voxelDown"
0265             || volPost == "VoxelStraight2" || volPost=="voxelStraight2"
0266             || volPost == "VoxelRight2" || volPost=="voxelRight2"
0267             || volPost == "VoxelLeft2" || volPost=="voxelLeft2"
0268             || volPost == "VoxelUp2" || volPost=="voxelUp2"
0269             || volPost == "VoxelDown2" || volPost=="voxelDown2")
0270         {
0271             fLastMetVoxelCopyNumber = step->GetPostStepPoint()->GetTouchable()->GetCopyNumber();
0272         }
0273     }
0274 }
0275
0276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0277