| 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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publications:
29 // Med. Phys. 45, (2018) e722-e739
30 // Phys. Med. 31 (2015) 861-874
31 // Med. Phys. 37 (2010) 4692-4708
32 // Int. J. Model. Simul. Sci. Comput. 1 (2010) 157–178
33 // The Geant4-DNA web site is available at http://geant4-dna.org
34 //
35 //
36 /// \file SteppingAction.cc
37 /// \brief Implementation of the SteppingAction class
38
39 #include "SteppingAction.hh"
40
41 #include "G4Alpha.hh"
42 #include "G4AnalysisManager.hh"
43 #include "G4DNAGenericIonsManager.hh"
44 #include "G4Electron.hh"
45 #include "G4EventManager.hh"
46 #include "G4Gamma.hh"
47 #include "G4Proton.hh"
48 #include "G4SystemOfUnits.hh"
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 SteppingAction::SteppingAction() : G4UserSteppingAction() {}
53
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55
56 SteppingAction::~SteppingAction() {}
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59
60 void SteppingAction::UserSteppingAction(const G4Step* step)
61 {
62 // Protection
63 if (!step->GetPostStepPoint()) return;
64 if (!step->GetPostStepPoint()->GetProcessDefinedStep()) return;
65
66 //
67 G4double flagParticle = -1.;
68 G4double flagProcess = -1.;
69 G4double x, y, z, xp, yp, zp;
70
71 // Particle identification
72
73 // The following method avoids the usage of string comparison
74
75 G4ParticleDefinition* partDef = step->GetTrack()->GetDynamicParticle()->GetDefinition();
76
77 if (partDef == G4Gamma::GammaDefinition()) flagParticle = 0;
78
79 if (partDef == G4Electron::ElectronDefinition()) flagParticle = 1;
80
81 if (partDef == G4Proton::ProtonDefinition()) flagParticle = 2;
82
83 if (partDef == G4Alpha::AlphaDefinition()) flagParticle = 4;
84
85 G4DNAGenericIonsManager* instance;
86 instance = G4DNAGenericIonsManager::Instance();
87
88 // Usage example
89 /*
90 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
91 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
92 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
93 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
94 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
95 */
96
97 if (partDef == instance->GetIon("hydrogen")) flagParticle = 3;
98
99 if (partDef == instance->GetIon("alpha+")) flagParticle = 5;
100
101 if (partDef == instance->GetIon("helium")) flagParticle = 6;
102
103 // Alternative method (based on string comparison)
104
105 /*
106 const G4String& particleName = step->GetTrack()->GetDynamicParticle()->
107 GetDefinition()->GetParticleName();
108
109 if (particleName == "gamma") flagParticle = 0;
110 else if (particleName == "e-") flagParticle = 1;
111 else if (particleName == "proton") flagParticle = 2;
112 else if (particleName == "hydrogen") flagParticle = 3;
113 else if (particleName == "alpha") flagParticle = 4;
114 else if (particleName == "alpha+") flagParticle = 5;
115 else if (particleName == "helium") flagParticle = 6;
116 */
117
118 // Process identification
119
120 // Process sub-types are listed in G4PhysicsListHelper.cc
121 // or in Geant4-DNA process class implementation files (*.cc)
122
123 G4StepPoint* preStep = step->GetPreStepPoint();
124 G4StepPoint* postStep = step->GetPostStepPoint();
125 G4int procID = postStep->GetProcessDefinedStep()->GetProcessSubType();
126
127 const G4String& processName = postStep->GetProcessDefinedStep()->GetProcessName();
128
129 if (processName == "Capture") flagProcess = 1;
130 // (no subType and procID exists at the moment for this process)
131
132 else if (flagParticle == 0) {
133 if (procID == 12)
134 flagProcess = 81;
135 else if (procID == 13)
136 flagProcess = 82;
137 else if (procID == 14)
138 flagProcess = 83;
139 else if (procID == 11)
140 flagProcess = 84;
141 }
142
143 else if (flagParticle == 1) {
144 if (procID == 58)
145 flagProcess = 10;
146 else if (procID == 51)
147 flagProcess = 11;
148 else if (procID == 52)
149 flagProcess = 12;
150 else if (procID == 53)
151 flagProcess = 13;
152 else if (procID == 55)
153 flagProcess = 14;
154 else if (procID == 54)
155 flagProcess = 15;
156 else if (procID == 10)
157 flagProcess = 110;
158 else if (procID == 1)
159 flagProcess = 120;
160 else if (procID == 2)
161 flagProcess = 130;
162 }
163
164 else if (flagParticle == 2) {
165 if (procID == 51)
166 flagProcess = 21;
167 else if (procID == 52)
168 flagProcess = 22;
169 else if (procID == 53)
170 flagProcess = 23;
171 else if (procID == 56)
172 flagProcess = 24;
173 else if (procID == 10)
174 flagProcess = 210;
175 else if (procID == 1)
176 flagProcess = 220;
177 else if (procID == 2)
178 flagProcess = 230;
179 else if (procID == 8)
180 flagProcess = 240;
181 }
182
183 else if (flagParticle == 3) {
184 if (procID == 51)
185 flagProcess = 31;
186 else if (procID == 52)
187 flagProcess = 32;
188 else if (procID == 53)
189 flagProcess = 33;
190 else if (procID == 57)
191 flagProcess = 35;
192 }
193
194 else if (flagParticle == 4) {
195 if (procID == 51)
196 flagProcess = 41;
197 else if (procID == 52)
198 flagProcess = 42;
199 else if (procID == 53)
200 flagProcess = 43;
201 else if (procID == 56)
202 flagProcess = 44;
203 else if (procID == 10)
204 flagProcess = 410;
205 else if (procID == 1)
206 flagProcess = 420;
207 else if (procID == 2)
208 flagProcess = 430;
209 else if (procID == 8)
210 flagProcess = 440;
211 }
212
213 else if (flagParticle == 5) {
214 if (procID == 51)
215 flagProcess = 51;
216 else if (procID == 52)
217 flagProcess = 52;
218 else if (procID == 53)
219 flagProcess = 53;
220 else if (procID == 56)
221 flagProcess = 54;
222 else if (procID == 57)
223 flagProcess = 55;
224 else if (procID == 10)
225 flagProcess = 510;
226 else if (procID == 1)
227 flagProcess = 520;
228 else if (procID == 2)
229 flagProcess = 530;
230 else if (procID == 8)
231 flagProcess = 540;
232 }
233
234 else if (flagParticle == 6) {
235 if (procID == 51)
236 flagProcess = 61;
237 else if (procID == 52)
238 flagProcess = 62;
239 else if (procID == 53)
240 flagProcess = 63;
241 else if (procID == 57)
242 flagProcess = 65;
243 }
244
245 else if (processName == "GenericIon_G4DNAIonisation")
246 flagProcess = 73;
247 else if (processName == "msc")
248 flagProcess = 710;
249 else if (processName == "CoulombScat")
250 flagProcess = 720;
251 else if (processName == "ionIoni")
252 flagProcess = 730;
253 else if (processName == "nuclearStopping")
254 flagProcess = 740;
255 // (for all GenericIons)
256
257 // Alternatively, using process names
258
259 /*
260 else if (processName=="e-_G4DNAElectronSolvation") flagProcess =10;
261 else if (processName=="e-_G4DNAElastic") flagProcess =11;
262 else if (processName=="e-_G4DNAExcitation") flagProcess =12;
263 else if (processName=="e-_G4DNAIonisation") flagProcess =13;
264 else if (processName=="e-_G4DNAAttachment") flagProcess =14;
265 else if (processName=="e-_G4DNAVibExcitation") flagProcess =15;
266
267 else if (processName=="proton_G4DNAElastic") flagProcess =21;
268 else if (processName=="proton_G4DNAExcitation") flagProcess =22;
269 else if (processName=="proton_G4DNAIonisation") flagProcess =23;
270 else if (processName=="proton_G4DNAChargeDecrease") flagProcess =24;
271
272 else if (processName=="hydrogen_G4DNAElastic") flagProcess =31;
273 else if (processName=="hydrogen_G4DNAExcitation") flagProcess =32;
274 else if (processName=="hydrogen_G4DNAIonisation") flagProcess =33;
275 else if (processName=="hydrogen_G4DNAChargeIncrease") flagProcess =35;
276
277 else if (processName=="alpha_G4DNAElastic") flagProcess =41;
278 else if (processName=="alpha_G4DNAExcitation") flagProcess =42;
279 else if (processName=="alpha_G4DNAIonisation") flagProcess =43;
280 else if (processName=="alpha_G4DNAChargeDecrease") flagProcess =44;
281
282 else if (processName=="alpha+_G4DNAElastic") flagProcess =51;
283 else if (processName=="alpha+_G4DNAExcitation") flagProcess =52;
284 else if (processName=="alpha+_G4DNAIonisation") flagProcess =53;
285 else if (processName=="alpha+_G4DNAChargeDecrease") flagProcess =54;
286 else if (processName=="alpha+_G4DNAChargeIncrease") flagProcess =55;
287
288 else if (processName=="helium_G4DNAElastic") flagProcess =61;
289 else if (processName=="helium_G4DNAExcitation") flagProcess =62;
290 else if (processName=="helium_G4DNAIonisation") flagProcess =63;
291 else if (processName=="helium_G4DNAChargeIncrease") flagProcess =65;
292
293 else if (processName=="GenericIon_G4DNAIonisation") flagProcess =73;
294
295 */
296
297 if (processName != "Transportation") {
298 x = preStep->GetPosition().x() / nanometer;
299 y = preStep->GetPosition().y() / nanometer;
300 z = preStep->GetPosition().z() / nanometer;
301
302 xp = postStep->GetPosition().x() / nanometer;
303 yp = postStep->GetPosition().y() / nanometer;
304 zp = postStep->GetPosition().z() / nanometer;
305
306 // get analysis manager
307
308 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
309
310 // fill ntuple
311 analysisManager->FillNtupleDColumn(0, flagParticle);
312 analysisManager->FillNtupleDColumn(1, flagProcess);
313 analysisManager->FillNtupleDColumn(2, xp);
314 analysisManager->FillNtupleDColumn(3, yp);
315 analysisManager->FillNtupleDColumn(4, zp);
316 analysisManager->FillNtupleDColumn(5, step->GetTotalEnergyDeposit() / eV);
317
318 analysisManager->FillNtupleDColumn(
319 6, std::sqrt((x - xp) * (x - xp) + (y - yp) * (y - yp) + (z - zp) * (z - zp)));
320
321 analysisManager->FillNtupleDColumn(
322 7, (preStep->GetKineticEnergy() - postStep->GetKineticEnergy()) / eV);
323
324 analysisManager->FillNtupleDColumn(8, preStep->GetKineticEnergy() / eV);
325
326 analysisManager->FillNtupleDColumn(9, preStep->GetMomentumDirection()
327 * postStep->GetMomentumDirection());
328
329 analysisManager->FillNtupleIColumn(
330 10, G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID());
331
332 analysisManager->FillNtupleIColumn(11, step->GetTrack()->GetTrackID());
333
334 analysisManager->FillNtupleIColumn(12, step->GetTrack()->GetParentID());
335
336 analysisManager->FillNtupleIColumn(13, step->GetTrack()->GetCurrentStepNumber());
337
338 analysisManager->AddNtupleRow();
339 }
340 }
341