| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file SteppingAction.cc
28 /// \brief Implementation of the SteppingAction class
29 /// \file SteppingAction.cc
30 /// \brief Implementation of the SteppingAction class
31
32 #include "SteppingAction.hh"
33 #include "Analysis.hh"
34 #include "G4SteppingManager.hh"
35 #include "G4VTouchable.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4RunManager.hh"
38 #include "G4LogicalVolumeStore.hh"
39 #include "G4SteppingManager.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Track.hh"
42 #include "G4RunManager.hh"
43 #include "G4Proton.hh"
44 #include "G4Electron.hh"
45 #include "G4Alpha.hh"
46 #include "G4DNAGenericIonsManager.hh"
47
48 #ifdef USE_MPI
49 #include "G4MPImanager.hh"
50 #endif
51
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54 SteppingAction::SteppingAction(EventAction* pEvent)
55 {
56 fEventAction = pEvent;
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61 void SteppingAction::UserSteppingAction(const G4Step* step)
62 {
63 SetupFlags(step);
64
65 SetupVoxelCopyNumber(step);
66
67 if(fFlagVolume>0)
68 {
69 fEventAction->AddEdep(step->GetTotalEnergyDeposit()/eV);
70 }
71 //G4cout<<step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName()<<" ; "<<step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName()<<G4endl;
72 // Check we are in a DNA volume
73 if(fFlagVolume == 1 // d1
74 || fFlagVolume == 11 // p1
75 || fFlagVolume == 2 // d2
76 || fFlagVolume == 22 // p2
77 || fFlagVolume == 3 // cyto
78 || fFlagVolume == 4 // gua
79 || fFlagVolume == 5 // thy
80 || fFlagVolume == 6 // ade
81 || fFlagVolume == 7 // d1_w
82 || fFlagVolume == 71 // p1_w
83 || fFlagVolume == 8 // d2_w
84 || fFlagVolume == 81 // p2_w
85 || fFlagVolume == 9 // ade_w
86 || fFlagVolume == 10 // gua_w
87 || fFlagVolume == 13 // cyto_w
88 || fFlagVolume == 12) // thy_w
89 {
90 // *****************************************************
91 // Saving physical stage informations
92 // *****************************************************
93 if(step->GetPostStepPoint()->GetTouchable()->GetVolume() ) // avoid asking non existing information in postStep
94 {
95 G4double x=step->GetPreStepPoint()->GetPosition().x()/nanometer;
96 G4double y=step->GetPreStepPoint()->GetPosition().y()/nanometer;
97 G4double z=step->GetPreStepPoint()->GetPosition().z()/nanometer;
98
99 G4double dE = step->GetTotalEnergyDeposit()/eV;
100 G4double copyNo = G4double (step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetCopyNo() );
101
102 G4int eventId = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
103 #ifdef USE_MPI
104 auto g4MPI = G4MPImanager::GetManager();
105 if (g4MPI->IsSlave()) { // update eventID only for slave, cause rank_master=0
106 G4int rank = g4MPI->GetRank();
107 eventId += g4MPI->GetEventsInMaster() + (rank-1)*g4MPI->GetEventsInSlave();
108 }
109 #endif
110 // ***************************
111 // put to vector
112 // ***************************
113 InfoInPhysStage aInfo;
114 aInfo.fFlagParticle = fFlagParticle;
115 aInfo.fFlagParentID = fFlagParentID;
116 aInfo.fFlagProcess = fFlagProcess;
117 aInfo.fX = x;
118 aInfo.fY = y;
119 aInfo.fZ = z;
120 aInfo.fEdep = dE;
121 aInfo.fEventNumber = eventId;
122 aInfo.fVolumeName = fFlagVolume;
123 aInfo.fCopyNumber = copyNo;
124 aInfo.fLastMetVoxelCopyNum = fLastMetVoxelCopyNumber;
125 Analysis::GetAnalysis()->AddInfoInPhysStage(aInfo);
126 }
127 }
128 }
129
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132 void SteppingAction::SetupFlags(const G4Step* step)
133 {
134 fFlagParticle=0;
135 fFlagProcess=0;
136 fFlagParentID=0;
137 fFlagVolume=0;
138
139 fFlagParentID = step->GetTrack()->GetParentID();
140
141 SetupParticleAndProcessFlags(step);
142
143 fFlagVolume = SetupVolumeFlag(step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName() );
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148 void SteppingAction::SetupParticleAndProcessFlags(const G4Step *step)
149 {
150 fFlagParticle = 0;
151 fFlagProcess = 0;
152 auto partDef = step->GetTrack()->GetDynamicParticle()->GetDefinition();
153 auto* instance = G4DNAGenericIonsManager::Instance();
154
155 if (partDef == G4Electron::ElectronDefinition()) fFlagParticle = 1;
156 if (partDef == G4Proton::ProtonDefinition()) fFlagParticle = 2;
157 if (partDef == instance->GetIon("hydrogen")) fFlagParticle = 3;
158 if (partDef == G4Alpha::AlphaDefinition()) fFlagParticle = 4;
159 if (partDef == instance->GetIon("alpha+")) fFlagParticle = 5;
160 if (partDef == instance->GetIon("helium")) fFlagParticle = 6;
161
162 G4int procSubtype = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessSubType();
163
164
165 if (fFlagParticle == 1) { // e-
166 if (procSubtype == 58) fFlagProcess = 10; //"e-_G4DNAElectronSolvation"
167 if (procSubtype == 51) fFlagProcess = 11; //"e-_G4DNAElastic"
168 if (procSubtype == 52) fFlagProcess = 12; //"e-_G4DNAExcitation"
169 if (procSubtype == 53) fFlagProcess = 13; //"e-_G4DNAIonisation"
170 if (procSubtype == 55) fFlagProcess = 14; //"e-_G4DNAAttachment"
171 if (procSubtype == 54) fFlagProcess = 15; //"e-_G4DNAVibExcitation"
172 }
173
174 if (fFlagParticle == 2) { // Proton
175 if (procSubtype == 52) fFlagProcess = 17; //"proton_G4DNAExcitation"
176 if (procSubtype == 53) fFlagProcess = 18; //"proton_G4DNAIonisation"
177 if (procSubtype == 56) fFlagProcess = 19; //"proton_G4DNAChargeDecrease"
178 }
179
180 if (fFlagParticle == 3) { // hydrogen
181 if (procSubtype == 52) fFlagProcess = 20; //"hydrogen_G4DNAExcitation"
182 if (procSubtype == 53) fFlagProcess = 21; //"hydrogen_G4DNAIonisation"
183 if (procSubtype == 57) fFlagProcess = 22; //"hydrogen_G4DNAChargeIncrease"
184 }
185
186 if (fFlagParticle == 4) { // alpha
187 if (procSubtype == 52) fFlagProcess = 23; //"alpha_G4DNAExcitation"
188 if (procSubtype == 53) fFlagProcess = 24; //"alpha_G4DNAIonisation"
189 if (procSubtype == 56) fFlagProcess = 25; //"alpha_G4DNAChargeDecrease"
190 }
191
192 if (fFlagParticle == 5) { // alpha+
193 if (procSubtype == 52) fFlagProcess = 26; //"alpha+_G4DNAExcitation"
194 if (procSubtype == 53) fFlagProcess = 27; //"alpha+_G4DNAIonisation"
195 if (procSubtype == 56) fFlagProcess = 28; //"alpha+_G4DNAChargeDecrease"
196 if (procSubtype == 57) fFlagProcess = 29; //"alpha+_G4DNAChargeIncrease"
197 }
198
199 if (fFlagParticle == 6) { // helium
200 if (procSubtype == 52) fFlagProcess = 30; //"helium__G4DNAExcitation"
201 if (procSubtype == 53) fFlagProcess = 31; //"helium__G4DNAIonisation"
202 if (procSubtype == 57) fFlagProcess = 32; //"helium__G4DNAChargeIncrease"
203 }
204
205 }
206
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208
209 G4int SteppingAction::SetupVolumeFlag(const G4String& volumeName)
210 {
211 G4int flagVolume(-1);
212
213 if(volumeName=="deoxyribose1_phys") flagVolume = 1;
214 else if(volumeName=="phosphate1_phys") flagVolume = 11;
215 else if(volumeName=="deoxyribose2_phys") flagVolume = 2;
216 else if(volumeName=="phosphate2_phys") flagVolume = 22;
217 else if(volumeName=="base_cytosine_phys") flagVolume = 3;
218 else if(volumeName=="base_guanine_phys") flagVolume = 4;
219 else if(volumeName=="base_thymine_phys") flagVolume = 5;
220 else if(volumeName=="base_adenine_phys") flagVolume = 6;
221 else if(volumeName=="deoxyribose1_water_phys") flagVolume = 7;
222 else if(volumeName=="phosphate1_water_phys") flagVolume = 71;
223 else if(volumeName=="deoxyribose2_water_phys") flagVolume = 8;
224 else if(volumeName=="phosphate2_water_phys") flagVolume = 81;
225 else if(volumeName=="base_adenine_water_phys") flagVolume = 9;
226 else if(volumeName=="base_guanine_water_phys") flagVolume = 10;
227 else if(volumeName=="base_cytosine_water_phys") flagVolume = 13;
228 else if(volumeName=="base_thymine_water_phys") flagVolume = 12;
229 else if(volumeName=="fiber") flagVolume = 110;
230 else if(volumeName=="voxelStraight" || volumeName=="VoxelStraight") flagVolume = 161;
231 else if(volumeName=="voxelRight" || volumeName=="VoxelRight") flagVolume = 162;
232 else if(volumeName=="voxelLeft" || volumeName=="VoxelLeft") flagVolume = 163;
233 else if(volumeName=="voxelUp" || volumeName=="VoxelUp") flagVolume = 164;
234 else if(volumeName=="voxelDown" || volumeName=="VoxelDown") flagVolume = 165;
235 else if(volumeName=="voxelStraight2" || volumeName=="VoxelStraight2") flagVolume = 261;
236 else if(volumeName=="voxelRight2" || volumeName=="VoxelRight2") flagVolume = 262;
237 else if(volumeName=="voxelLeft2" || volumeName=="VoxelLeft2") flagVolume = 263;
238 else if(volumeName=="voxelUp2" || volumeName=="VoxelUp2") flagVolume = 264;
239 else if(volumeName=="voxelDown2" || volumeName=="VoxelDown2") flagVolume = 265;
240 else if(volumeName=="physWorld") flagVolume = -1.;
241 else if(volumeName=="wrapper") flagVolume = 130;
242 else if(volumeName=="histone_phys") flagVolume = 140;
243 else if(volumeName=="nucleus_pl") flagVolume = 200;
244 return flagVolume;
245 }
246
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
248
249 void SteppingAction::SetupVoxelCopyNumber(const G4Step* step)
250 {
251 // Each time we will do a step in a voxel, we will setup the lastMetVoxelCpNum flag.
252 // This way, this number will always be the one of the last voxel met by the particle.
253 // Therefore, in a DNA volume, the number will be the one of the container voxel.
254
255 if(step->GetPostStepPoint()->GetTouchable()->GetVolume() ) // avoid asking non existing information in postStep
256 {
257 const G4String& volPost = step->GetPostStepPoint()->GetTouchable()->GetVolume()->GetName();
258
259 // If we enter a voxel or stay in it
260 if( volPost == "VoxelStraight" || volPost=="voxelStraight"
261 || volPost == "VoxelRight" || volPost=="voxelRight"
262 || volPost == "VoxelLeft" || volPost=="voxelLeft"
263 || volPost == "VoxelUp" || volPost=="voxelUp"
264 || volPost == "VoxelDown" || volPost=="voxelDown"
265 || volPost == "VoxelStraight2" || volPost=="voxelStraight2"
266 || volPost == "VoxelRight2" || volPost=="voxelRight2"
267 || volPost == "VoxelLeft2" || volPost=="voxelLeft2"
268 || volPost == "VoxelUp2" || volPost=="voxelUp2"
269 || volPost == "VoxelDown2" || volPost=="voxelDown2")
270 {
271 fLastMetVoxelCopyNumber = step->GetPostStepPoint()->GetTouchable()->GetCopyNumber();
272 }
273 }
274 }
275
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277
278