source: trunk/examples/advanced/radioprotection/src/RemSimSteppingAction.cc@ 1302

Last change on this file since 1302 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 14.1 KB
RevLine 
[807]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// $Id: RemSimSteppingAction.cc,v 1.10 2006/06/29 16:24:25 gunter Exp $
[1230]28// GEANT4 tag $Name: geant4-09-03-cand-01 $
[807]29//
30//
31// Author: Susanna Guatelli (guatelli@ge.infn.it)
32//
33
34#include "G4ios.hh"
35#include "G4SteppingManager.hh"
36#include "G4Step.hh"
37#include "G4Track.hh"
38#include "G4StepPoint.hh"
39#include "G4ParticleDefinition.hh"
40#include "G4VPhysicalVolume.hh"
41#include "RemSimSteppingAction.hh"
42#include "RemSimPrimaryGeneratorAction.hh"
43#include "G4TrackStatus.hh"
44#include "RemSimSteppingActionMessenger.hh"
45#ifdef G4ANALYSIS_USE
46#include "RemSimAnalysisManager.hh"
47#endif
48
49RemSimSteppingAction::RemSimSteppingAction(RemSimPrimaryGeneratorAction* primary):
50 primaryAction(primary)
51{
52 hadronic = "Off";
53 messenger = new RemSimSteppingActionMessenger(this);
54}
55
56RemSimSteppingAction::~RemSimSteppingAction()
57{
58 delete messenger;
59}
60
61void RemSimSteppingAction::UserSteppingAction(const G4Step* aStep)
62{
63 G4String oldVolumeName = aStep -> GetPreStepPoint()->
64 GetPhysicalVolume()-> GetName();
65
66 // Store in histograms useful information concerning primary particles
67
68#ifdef G4ANALYSIS_USE
69 RemSimAnalysisManager* analysis = RemSimAnalysisManager::getInstance();
70
71 G4VPhysicalVolume* volume = aStep -> GetPostStepPoint() -> GetPhysicalVolume();
72
73
74 // Particle reaching the astronaut
75 if(oldVolumeName != "phantom")
76 {
77 if (volume)
78 {
79 G4String volumeName = volume -> GetName();
80 G4String particleName = (aStep -> GetTrack() -> GetDynamicParticle()
81 -> GetDefinition() -> GetParticleName());
82
83 if (volumeName == "phantom")
84 {
85 G4double particleEnergy = aStep -> GetTrack() -> GetKineticEnergy();
86
87 if(aStep -> GetTrack() -> GetTrackID()== 1) // primary particles
88 {
89 G4double initialEnergy = primaryAction -> GetInitialEnergy();
90
91 if ( particleName == "proton" ||
92 particleName == "alpha" ||
93 particleName == "IonO16" ||
94 particleName == "IonC12" ||
95 particleName == "IonFe52"||
96 particleName == "IonSi28")
97 {
98 G4int baryon = aStep -> GetTrack() -> GetDefinition() -> GetBaryonNumber();
99
100 // Initial energy of primary particles impinging on the phantom
101 analysis -> PrimaryInitialEnergyIn((initialEnergy/baryon)/MeV);
102
103 // Energy of primary particles impinging on the phantom
104 analysis -> PrimaryEnergyIn((particleEnergy/baryon)/MeV);
105 }
106 }
107 // secondary particle reaching the astronaut
108 else {
109
110 // if i =0 secondary proton, i =1 neutron, i=2 pion, i=3 alpha,
111 // i =4 other, i=5 electron, i = 6 gamma, i=7 positrons,
112 // i=8 muons, i=9 neutrinos
113
114 if (particleName == "proton")
115 {
116 analysis -> SecondaryReachingThePhantom(0);
117 analysis -> SecondaryProtonReachingThePhantom(particleEnergy/MeV);
118 }
119 else
120 {
121 if (particleName == "neutron")
122 {
123 analysis -> SecondaryReachingThePhantom(1);
124 analysis -> SecondaryNeutronReachingThePhantom(particleEnergy/MeV);
125 }
126 else{
127 if (particleName == "pi+" ||
128 particleName == "pi-" ||particleName == "pi0" )
129 {
130 analysis -> SecondaryReachingThePhantom(2);
131 analysis -> SecondaryPionReachingThePhantom(particleEnergy/MeV);
132 }
133 else{
134 if (particleName == "alpha")
135 {
136 analysis -> SecondaryReachingThePhantom(3);
137 analysis -> SecondaryAlphaReachingThePhantom(particleEnergy/MeV);
138 }
139 else{
140 if (particleName == "e+")
141 {analysis -> SecondaryReachingThePhantom(7);
142 analysis -> SecondaryPositronReachingThePhantom(particleEnergy/MeV);
143 }
144 else{
145 if (particleName == "e-")
146 {analysis -> SecondaryReachingThePhantom(5);
147 analysis -> SecondaryElectronReachingThePhantom(particleEnergy/MeV);
148 }
149 else{
150 if (particleName == "gamma")
151 {analysis -> SecondaryReachingThePhantom(6);
152 analysis -> SecondaryGammaReachingThePhantom(particleEnergy/MeV);
153 }
154 else{
155 if (particleName == "mu+"
156 ||particleName == "mu-" )
157 {analysis -> SecondaryReachingThePhantom(8);
158 analysis -> SecondaryMuonReachingThePhantom(particleEnergy/MeV);
159 }
160 else {
161 if (particleName == "nu_e" || particleName == "nu_mu" ||
162 particleName == "anti_nu_e" || particleName == "anti_nu_mu")
163 { analysis -> SecondaryReachingThePhantom(9);
164 }
165 else{
166 analysis -> SecondaryReachingThePhantom(4);
167 analysis -> SecondaryOtherReachingThePhantom(particleEnergy/MeV);
168 }
169 }
170 }
171 }
172 }
173 }
174 }
175 }
176 }
177 }
178 }
179 }
180 }
181
182
183 // Primar particles outgoing the phantom
184 if(oldVolumeName == "phantom")
185 {
186 if (volume)
187 {
188 G4String volumeName = volume -> GetName();
189 G4String particleName = (aStep -> GetTrack() -> GetDynamicParticle()
190 -> GetDefinition() -> GetParticleName());
191 if (volumeName != "phantom")
192 {
193 if(aStep -> GetTrack() -> GetTrackID()== 1) // primary particles
194 {
195 G4double initialEnergy = primaryAction -> GetInitialEnergy();
196 G4double particleEnergy = aStep->GetTrack() -> GetKineticEnergy();
197 if ( particleName == "proton" ||
198 particleName == "alpha" ||
199 particleName == "IonO16" ||
200 particleName == "IonC12" ||
201 particleName == "IonFe52"||
202 particleName == "IonSi28" )
203 {
204 G4int baryon = aStep -> GetTrack() -> GetDefinition() -> GetBaryonNumber();
205 // plot of MeV/nucl
206 analysis -> PrimaryInitialEnergyOut((initialEnergy/baryon)/MeV);
207 analysis -> PrimaryEnergyOut((particleEnergy/baryon)/MeV);
208 }
209 }
210 }
211 }
212 }
213
214 // Analysis of the secondary particles generated in the phantom
215
216 G4SteppingManager* steppingManager = fpSteppingManager;
217 G4Track* theTrack = aStep -> GetTrack();
218
219 // check if it is alive
220 if(theTrack-> GetTrackStatus() == fAlive) { return; }
221
222 // Retrieve the secondary particles
223 G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
224
225 for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
226 {
227 // Retrieve the info about the generation of secondary particles in the phantom and
228 // in the vehicle
229 G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName();
230 G4String secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();
231 G4double secondaryParticleKineticEnergy = (*fSecondary)[lp1] -> GetKineticEnergy();
232 G4String process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName();
233
234 // If the secondaries are originated in the phantom....
235 if (volumeName == "phantom")
236 {
237
238 G4double slice = (*fSecondary)[lp1] -> GetPosition().z() ;
239 //G4cout << "Secondary particle in phantom: " << secondaryParticleName
240 // << "energy (MeV): " << secondaryParticleKineticEnergy << G4endl;
241 // << "creation slice: " << slice << G4endl;
242 // if i =0 secondary proton, i =1 neutron, i=2 pion, i=3 alpha,
243 // i =4 other, i=5 electron, i = 6 gamma, i=7 positrons,
244 // i=8 muons, i=9 neutrinos
245 G4double translation = 15. * cm;
246 if (secondaryParticleName == "proton")
247 {
248 analysis -> SecondaryInPhantom(0);
249 analysis -> SecondaryProtonInPhantom(secondaryParticleKineticEnergy/MeV);
250 slice = slice + translation;
251 analysis -> SecondaryProtonInPhantomSlice(slice/cm); }
252 else
253 {
254 if (secondaryParticleName == "neutron")
255 {
256 analysis -> SecondaryInPhantom(1);
257 analysis -> SecondaryNeutronInPhantom(secondaryParticleKineticEnergy/MeV);
258 slice = slice + translation;
259 analysis -> SecondaryNeutronInPhantomSlice(slice/cm);
260 }
261 else{
262 if (secondaryParticleName == "pi+" ||
263 secondaryParticleName == "pi-"||secondaryParticleName == "pi0")
264 {
265 analysis -> SecondaryInPhantom(2);
266 analysis -> SecondaryPionInPhantom(secondaryParticleKineticEnergy/MeV);
267 slice = slice + translation;
268 analysis -> SecondaryPionInPhantomSlice(slice/cm);
269 }
270 else{
271 if (secondaryParticleName == "alpha")
272 {
273 analysis -> SecondaryInPhantom(3);
274 analysis -> SecondaryAlphaInPhantom(secondaryParticleKineticEnergy/MeV);
275 slice = slice + translation;
276 analysis -> SecondaryAlphaInPhantomSlice(slice/cm);
277 }
278 else{
279 if (secondaryParticleName == "e+")
280 {analysis -> SecondaryInPhantom(7);
281 analysis -> SecondaryPositronInPhantom(secondaryParticleKineticEnergy/MeV);
282 slice = slice + translation;
283 analysis -> SecondaryPositronInPhantomSlice(slice/cm);
284 }
285 else{
286 if (secondaryParticleName == "e-")
287 {analysis -> SecondaryInPhantom(5);
288 analysis -> SecondaryElectronInPhantom(secondaryParticleKineticEnergy/MeV);
289 slice = slice + translation;
290 analysis -> SecondaryElectronInPhantomSlice(slice/cm);
291 }
292 else{
293 if (secondaryParticleName == "gamma")
294 {analysis -> SecondaryInPhantom(6);
295 analysis -> SecondaryGammaInPhantom(secondaryParticleKineticEnergy/MeV);
296 slice = slice + translation;
297 analysis -> SecondaryGammaInPhantomSlice(slice/cm);
298 }
299 else{
300 if (secondaryParticleName == "mu+"
301 ||secondaryParticleName == "mu-" )
302 {analysis -> SecondaryInPhantom(8);
303 analysis -> SecondaryMuonInPhantom(secondaryParticleKineticEnergy/MeV);
304 slice = slice + translation;
305 analysis -> SecondaryMuonInPhantomSlice(slice/cm);
306 }
307 else{
308 if (secondaryParticleName == "nu_e" || secondaryParticleName == "nu_mu"
309 || secondaryParticleName == "anti_nu_e" || secondaryParticleName == "anti_nu_mu")
310 {
311 analysis -> SecondaryInPhantom(9);
312 analysis -> SecondaryNeutrinoInPhantom(secondaryParticleKineticEnergy/MeV);}
313 else{
314 analysis -> SecondaryInPhantom(4);
315 analysis -> SecondaryOtherInPhantom(secondaryParticleKineticEnergy/MeV);
316 slice = slice + translation;
317 analysis -> SecondaryOtherInPhantomSlice(slice/cm);
318 }
319 }
320 }
321 }
322 }
323 }
324
325 }
326 }
327 }
328 }
329
330 // secondary particles produced in the multilayer + shielding
331 else{
332 if (volumeName != "world")
333 {
334 if (secondaryParticleName == "proton") analysis -> SecondaryInVehicle(0);
335 else
336 {
337 if (secondaryParticleName == "neutron") analysis -> SecondaryInVehicle(1);
338 else{
339 if (secondaryParticleName == "pi+" ||
340 secondaryParticleName == "pi-"||secondaryParticleName == "pi0") analysis -> SecondaryInVehicle(2);
341 else{
342 if (secondaryParticleName == "alpha") analysis -> SecondaryInVehicle(3);
343
344 else{
345 if (secondaryParticleName == "e+") analysis -> SecondaryInVehicle(7);
346
347 else{
348 if (secondaryParticleName == "e-") analysis -> SecondaryInVehicle(5);
349 else{
350 if (secondaryParticleName == "gamma") analysis -> SecondaryInVehicle(6);
351 else{
352 if (secondaryParticleName == "mu+"
353 ||secondaryParticleName == "mu-" ) analysis -> SecondaryInVehicle(8);
354 else{
355 if (secondaryParticleName == "nu_e"|| secondaryParticleName == "nu_mu" ||
356 secondaryParticleName == "anti_nu_e" || secondaryParticleName == "anti_nu_mu")
357 {analysis -> SecondaryInVehicle(9);}
358
359 else{
360 analysis -> SecondaryInVehicle(4);
361 }
362 }
363 }
364 }
365 }
366 }
367 }
368 }
369 }
370 }
371 }
372 }
373
374#endif
375
376 //analysis of hadronic physics
377 if (hadronic == "On")
378 {
379 if(aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL)
380 {
381 G4String process = aStep -> GetPostStepPoint() ->
382 GetProcessDefinedStep() ->GetProcessName();
383
384 if ((process != "Transportation") &&
385 (process != "msc") &&
386 (process != "LowEnergyIoni") &&
387 (process != "LowEnergyBrem") &&
388 (process != "eIoni") &&
389 (process != "hIoni") &&
390 (process != "eBrem") &&
391 (process != "compt") &&
392 (process != "phot") &&
393 (process != "conv") &&
394 (process != "annihil") &&
395 (process != "hLowEIoni") &&
396 (process != "LowEnBrem") &&
397 (process != "LowEnCompton") &&
398 (process != "LowEnPhotoElec") &&
399 (process != "LowEnRayleigh") &&
400 (process != "LowEnConversion"))
401 G4cout << "Hadronic Process:" << process << G4endl;
402 }
403 }
404}
405void RemSimSteppingAction::SetHadronicAnalysis(G4String value)
406{
407 hadronic = value;
408}
Note: See TracBrowser for help on using the repository browser.