Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-31 09:22:03

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