ATLAS Offline Software
Loading...
Searching...
No Matches
RadiationMapsMaker.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "RadiationMapsMaker.h"
6#include <iostream>
7#include <string>
8#include <cmath>
9#include "G4Electron.hh"
10#include "G4Positron.hh"
11#include "G4PionPlus.hh"
12#include "G4PionMinus.hh"
13#include "G4PionZero.hh"
14#include "G4KaonPlus.hh"
15#include "G4KaonMinus.hh"
16#include "G4KaonZeroShort.hh"
17#include "G4KaonZeroLong.hh"
18#include "G4KaonZero.hh"
19#include "G4AntiKaonZero.hh"
20#include "G4MuonMinus.hh"
21#include "G4MuonPlus.hh"
22#include "G4Neutron.hh"
23#include "G4AntiNeutron.hh"
24#include "G4Proton.hh"
25#include "G4AntiProton.hh"
26#include "G4Lambda.hh"
27#include "G4AntiLambda.hh"
28#include "G4SigmaPlus.hh"
29#include "G4AntiSigmaPlus.hh"
30#include "G4SigmaMinus.hh"
31#include "G4AntiSigmaMinus.hh"
32#include "G4SigmaZero.hh"
33#include "G4AntiSigmaZero.hh"
34#include "G4XiZero.hh"
35#include "G4AntiXiZero.hh"
36#include "G4XiMinus.hh"
37#include "G4AntiXiMinus.hh"
38#include "G4OmegaMinus.hh"
39#include "G4AntiOmegaMinus.hh"
40#include "G4Gamma.hh"
41#include "G4Alpha.hh"
42#include "G4Deuteron.hh"
43#include "G4Triton.hh"
44#include "G4He3.hh"
45#include "G4Geantino.hh"
46#include "TGraph.h"
47#include "G4ParticleTable.hh"
48#include "G4IonTable.hh"
49#include "G4VProcess.hh"
50#include "G4Step.hh"
51#include "G4StepPoint.hh"
52#include "G4VSensitiveDetector.hh"
53#include "Randomize.hh"
54#include "PathResolver/PathResolver.h"
55
56namespace G4UA{
57
58
59 //----------------------------------------------------------------------------
60 // Constructor
61 //----------------------------------------------------------------------------
66
67 //----------------------------------------------------------------------------
68 // Destructor
69 //----------------------------------------------------------------------------
70 RadiationMapsMaker::~RadiationMapsMaker()
71 {
72 // close activation file if open
73 if ( m_activationFile.is_open() ) {
74 m_activationFile.close();
75 }
76
77 // delete TGraph objects
78
79 // NIEL weights
80 if ( m_tgpSiA ) {
81 m_tgpSiA->Delete();
82 m_tgpSiA = 0;
83 }
84 if ( m_tgpSiB ) {
85 m_tgpSiB->Delete();
86 m_tgpSiB = 0;
87 }
88 if ( m_tgnSiA ) {
89 m_tgnSiA->Delete();
90 m_tgnSiA = 0;
91 }
92 if ( m_tgnSiB ) {
93 m_tgnSiB->Delete();
94 m_tgnSiB = 0;
95 }
96 if ( m_tgnSiC ) {
97 m_tgnSiC->Delete();
98 m_tgnSiC = 0;
99 }
100 if ( m_tgpiSi ) {
101 m_tgpiSi->Delete();
102 m_tgpiSi = 0;
103 }
104 if ( m_tgeSi ) {
105 m_tgeSi->Delete();
106 m_tgeSi = 0;
107 }
108
109 // H_T weights
110 if ( m_tgHn ) {
111 m_tgHn->Delete();
112 m_tgHn = 0;
113 }
114 if ( m_tgHg ) {
115 m_tgHg->Delete();
116 m_tgHg = 0;
117 }
118 if ( m_tgHp ) {
119 m_tgHp->Delete();
120 m_tgHp = 0;
121 }
122 if ( m_tgHem ) {
123 m_tgHem->Delete();
124 m_tgHem = 0;
125 }
126 if ( m_tgHep ) {
127 m_tgHep->Delete();
128 m_tgHep = 0;
129 }
130 if ( m_tgHmm ) {
131 m_tgHmm->Delete();
132 m_tgHmm = 0;
133 }
134 if ( m_tgHmp ) {
135 m_tgHmp->Delete();
136 m_tgHmp = 0;
137 }
138 if ( m_tgHpm ) {
139 m_tgHpm->Delete();
140 m_tgHpm = 0;
141 }
142 if ( m_tgHpp ) {
143 m_tgHpp->Delete();
144 m_tgHpp = 0;
145 }
146 if ( m_tgHa ) {
147 m_tgHa->Delete();
148 m_tgHa = 0;
149 }
150 }
151
152 //---------------------------------------------------------------------------
153 // Merge results
154 //---------------------------------------------------------------------------
155 void RadiationMapsMaker::Report::merge(const RadiationMapsMaker::Report& maps)
156 {
157
158 // all 2d vectors have same size
159 for(unsigned int i=0;i<maps.m_rz_tid.size();i++) {
160 // zoom 2d
161 m_rz_tid [i] += maps.m_rz_tid [i];
162 m_rz_eion[i] += maps.m_rz_eion[i];
163 m_rz_niel[i] += maps.m_rz_niel[i];
164 m_rz_h20 [i] += maps.m_rz_h20 [i];
165 m_rz_neut[i] += maps.m_rz_neut[i];
166 m_rz_chad[i] += maps.m_rz_chad[i];
167
168 // full 2d
169 m_full_rz_tid [i] += maps.m_full_rz_tid [i];
170 m_full_rz_eion[i] += maps.m_full_rz_eion[i];
171 m_full_rz_niel[i] += maps.m_full_rz_niel[i];
172 m_full_rz_h20 [i] += maps.m_full_rz_h20 [i];
173 m_full_rz_neut[i] += maps.m_full_rz_neut[i];
174 m_full_rz_chad[i] += maps.m_full_rz_chad[i];
175 }
176
177 for(unsigned int i=0;i<maps.m_rz_vol.size();i++) {
178 // need volume fraction only if particular materials are selected
179 m_rz_vol [i] += maps.m_rz_vol [i];
180 m_rz_norm[i] += maps.m_rz_norm[i];
181 m_full_rz_vol [i] += maps.m_full_rz_vol [i];
182 m_full_rz_norm[i] += maps.m_full_rz_norm[i];
183 }
184
185 // all zoom 3d vectors have same size
186 for(unsigned int i=0;i<maps.m_3d_tid.size();i++) {
187 // zoom 3d
188 m_3d_tid [i] += maps.m_3d_tid [i];
189 m_3d_eion[i] += maps.m_3d_eion[i];
190 m_3d_niel[i] += maps.m_3d_niel[i];
191 m_3d_h20 [i] += maps.m_3d_h20 [i];
192 m_3d_neut[i] += maps.m_3d_neut[i];
193 m_3d_chad[i] += maps.m_3d_chad[i];
194 }
195
196 // all time 2d vectors have same size
197 for(unsigned int i=0;i<maps.m_rz_tid_time.size();i++) {
198 // zoom 2d
199 m_rz_tid_time[i] += maps.m_rz_tid_time[i];
200 m_rz_ht_time [i] += maps.m_rz_ht_time [i];
201
202 // full 2d
203 m_full_rz_tid_time[i] += maps.m_full_rz_tid_time[i];
204 m_full_rz_ht_time [i] += maps.m_full_rz_ht_time [i];
205 }
206
207 // all element fraction 2d vectors have same size
208 for(unsigned int i=0;i<maps.m_rz_element.size();i++) {
209 m_rz_element [i] += maps.m_rz_element [i];
210 m_full_rz_element[i] += maps.m_full_rz_element[i];
211 }
212
213 for(unsigned int i=0;i<maps.m_3d_vol.size();i++) {
214 // need volume fraction only if particular materials are selected
215 m_3d_vol [i] += maps.m_3d_vol [i];
216 m_3d_norm[i] += maps.m_3d_norm[i];
217 }
218
219 // neutron spectra have different size from all other particle's spectra
220 for(unsigned int i=0;i<maps.m_rz_neut_spec.size();i++) {
221 m_rz_neut_spec [i] += maps.m_rz_neut_spec [i];
222 m_full_rz_neut_spec[i] += maps.m_full_rz_neut_spec[i];
223 }
224 // x theta bins
225 for(unsigned int i=0;i<maps.m_theta_full_rz_neut_spec.size();i++) {
226 m_theta_full_rz_neut_spec[i] += maps.m_theta_full_rz_neut_spec[i];
227 }
228 // all other particle's spectra have same size
229 for(unsigned int i=0;i<maps.m_rz_gamm_spec.size();i++) {
230 m_rz_gamm_spec [i] += maps.m_rz_gamm_spec [i];
231 m_full_rz_gamm_spec[i] += maps.m_full_rz_gamm_spec[i];
232 m_rz_elec_spec [i] += maps.m_rz_elec_spec [i];
233 m_full_rz_elec_spec[i] += maps.m_full_rz_elec_spec[i];
234 m_rz_muon_spec [i] += maps.m_rz_muon_spec [i];
235 m_full_rz_muon_spec[i] += maps.m_full_rz_muon_spec[i];
236 m_rz_pion_spec [i] += maps.m_rz_pion_spec [i];
237 m_full_rz_pion_spec[i] += maps.m_full_rz_pion_spec[i];
238 m_rz_prot_spec [i] += maps.m_rz_prot_spec [i];
239 m_full_rz_prot_spec[i] += maps.m_full_rz_prot_spec[i];
240 m_rz_rest_spec [i] += maps.m_rz_rest_spec [i];
241 m_full_rz_rest_spec[i] += maps.m_full_rz_rest_spec[i];
242 }
243 // x theta bins
244 for(unsigned int i=0;i<maps.m_theta_full_rz_gamm_spec.size();i++) {
245 m_theta_full_rz_gamm_spec[i] += maps.m_theta_full_rz_gamm_spec[i];
246 m_theta_full_rz_elec_spec[i] += maps.m_theta_full_rz_elec_spec[i];
247 m_theta_full_rz_muon_spec[i] += maps.m_theta_full_rz_muon_spec[i];
248 m_theta_full_rz_pion_spec[i] += maps.m_theta_full_rz_pion_spec[i];
249 m_theta_full_rz_prot_spec[i] += maps.m_theta_full_rz_prot_spec[i];
250 m_theta_full_rz_rchgd_spec[i] += maps.m_theta_full_rz_rchgd_spec[i];
251 m_theta_full_rz_rneut_spec[i] += maps.m_theta_full_rz_rneut_spec[i];
252 }
253 }
254
255 void RadiationMapsMaker::BeginOfRunAction(const G4Run*){
256
257 // open activation file if wanted
258 if ( !m_config.activationFileName.empty() ) {
259 m_activationFile.open(m_config.activationFileName,std::ios_base::app);
260 }
261 // first make sure the vectors are empty
262
263 m_maps.m_rz_tid .resize(0);
264 m_maps.m_rz_eion.resize(0);
265 m_maps.m_rz_niel.resize(0);
266 m_maps.m_rz_h20 .resize(0);
267 m_maps.m_rz_neut.resize(0);
268 m_maps.m_rz_chad.resize(0);
269
270 m_maps.m_full_rz_tid .resize(0);
271 m_maps.m_full_rz_eion.resize(0);
272 m_maps.m_full_rz_niel.resize(0);
273 m_maps.m_full_rz_h20 .resize(0);
274 m_maps.m_full_rz_neut.resize(0);
275 m_maps.m_full_rz_chad.resize(0);
276
277 m_maps.m_3d_tid .resize(0);
278 m_maps.m_3d_eion.resize(0);
279 m_maps.m_3d_niel.resize(0);
280 m_maps.m_3d_h20 .resize(0);
281 m_maps.m_3d_neut.resize(0);
282 m_maps.m_3d_chad.resize(0);
283
284 m_maps.m_rz_tid_time .resize(0);
285 m_maps.m_rz_ht_time .resize(0);
286 m_maps.m_full_rz_tid_time .resize(0);
287 m_maps.m_full_rz_ht_time .resize(0);
288
289 m_maps.m_rz_element .resize(0);
290 m_maps.m_full_rz_element .resize(0);
291
292 m_maps.m_rz_neut_spec .resize(0);
293 m_maps.m_full_rz_neut_spec.resize(0);
294 m_maps.m_rz_gamm_spec .resize(0);
295 m_maps.m_full_rz_gamm_spec.resize(0);
296 m_maps.m_rz_elec_spec .resize(0);
297 m_maps.m_full_rz_elec_spec.resize(0);
298 m_maps.m_rz_muon_spec .resize(0);
299 m_maps.m_full_rz_muon_spec.resize(0);
300 m_maps.m_rz_pion_spec .resize(0);
301 m_maps.m_full_rz_pion_spec.resize(0);
302 m_maps.m_rz_prot_spec .resize(0);
303 m_maps.m_full_rz_prot_spec.resize(0);
304 m_maps.m_rz_rest_spec .resize(0);
305 m_maps.m_full_rz_rest_spec.resize(0);
306
307 m_maps.m_theta_full_rz_neut_spec .resize(0);
308 m_maps.m_theta_full_rz_gamm_spec .resize(0);
309 m_maps.m_theta_full_rz_elec_spec .resize(0);
310 m_maps.m_theta_full_rz_muon_spec .resize(0);
311 m_maps.m_theta_full_rz_pion_spec .resize(0);
312 m_maps.m_theta_full_rz_prot_spec .resize(0);
313 m_maps.m_theta_full_rz_rchgd_spec.resize(0);
314 m_maps.m_theta_full_rz_rneut_spec.resize(0);
315
316 if (!m_config.materials.empty()) {
317 // need volume fraction only if particular materials is selected
318 m_maps.m_rz_vol .resize(0);
319 m_maps.m_rz_norm.resize(0);
320 m_maps.m_full_rz_vol .resize(0);
321 m_maps.m_full_rz_norm.resize(0);
322 m_maps.m_3d_vol .resize(0);
323 m_maps.m_3d_norm.resize(0);
324 }
325
326 // then resize to proper size and initialize with 0's
327
328 m_maps.m_rz_tid .resize(m_config.nBinsz*m_config.nBinsr,0.0);
329 m_maps.m_rz_eion.resize(m_config.nBinsz*m_config.nBinsr,0.0);
330 m_maps.m_rz_niel.resize(m_config.nBinsz*m_config.nBinsr,0.0);
331 m_maps.m_rz_h20 .resize(m_config.nBinsz*m_config.nBinsr,0.0);
332 m_maps.m_rz_neut.resize(m_config.nBinsz*m_config.nBinsr,0.0);
333 m_maps.m_rz_chad.resize(m_config.nBinsz*m_config.nBinsr,0.0);
334
335 m_maps.m_full_rz_tid .resize(m_config.nBinsz*m_config.nBinsr,0.0);
336 m_maps.m_full_rz_eion.resize(m_config.nBinsz*m_config.nBinsr,0.0);
337 m_maps.m_full_rz_niel.resize(m_config.nBinsz*m_config.nBinsr,0.0);
338 m_maps.m_full_rz_h20 .resize(m_config.nBinsz*m_config.nBinsr,0.0);
339 m_maps.m_full_rz_neut.resize(m_config.nBinsz*m_config.nBinsr,0.0);
340 m_maps.m_full_rz_chad.resize(m_config.nBinsz*m_config.nBinsr,0.0);
341
342 m_maps.m_3d_tid .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
343 m_maps.m_3d_eion.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
344 m_maps.m_3d_niel.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
345 m_maps.m_3d_h20 .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
346 m_maps.m_3d_neut.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
347 m_maps.m_3d_chad.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
348
349 m_maps.m_rz_tid_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
350 m_maps.m_rz_ht_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
351 m_maps.m_full_rz_tid_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
352 m_maps.m_full_rz_ht_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
353
354 m_maps.m_rz_element .resize(m_config.nBinsz*m_config.nBinsr*(m_config.elemZMax-m_config.elemZMin+1),0.0);
355 m_maps.m_full_rz_element .resize(m_config.nBinsz*m_config.nBinsr*(m_config.elemZMax-m_config.elemZMin+1),0.0);
356
357 m_maps.m_rz_neut_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
358 m_maps.m_full_rz_neut_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
359 m_maps.m_rz_gamm_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
360 m_maps.m_full_rz_gamm_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
361 m_maps.m_rz_elec_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
362 m_maps.m_full_rz_elec_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
363 m_maps.m_rz_muon_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
364 m_maps.m_full_rz_muon_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
365 m_maps.m_rz_pion_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
366 m_maps.m_full_rz_pion_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
367 m_maps.m_rz_prot_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
368 m_maps.m_full_rz_prot_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
369 m_maps.m_rz_rest_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
370 m_maps.m_full_rz_rest_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
371
372 m_maps.m_theta_full_rz_neut_spec .resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
373 m_maps.m_theta_full_rz_gamm_spec .resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
374 m_maps.m_theta_full_rz_elec_spec .resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
375 m_maps.m_theta_full_rz_muon_spec .resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
376 m_maps.m_theta_full_rz_pion_spec .resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
377 m_maps.m_theta_full_rz_prot_spec .resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
378 m_maps.m_theta_full_rz_rchgd_spec.resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
379 m_maps.m_theta_full_rz_rneut_spec.resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
380
381 if (!m_config.materials.empty()) {
382 // need volume fraction only if particular materials are selected
383 // 2d zoom
384 m_maps.m_rz_vol .resize(m_config.nBinsz*m_config.nBinsr,0.0);
385 m_maps.m_rz_norm.resize(m_config.nBinsz*m_config.nBinsr,0.0);
386 // 2d full
387 m_maps.m_full_rz_vol .resize(m_config.nBinsz*m_config.nBinsr,0.0);
388 m_maps.m_full_rz_norm.resize(m_config.nBinsz*m_config.nBinsr,0.0);
389 // 3d
390 m_maps.m_3d_vol .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
391 m_maps.m_3d_norm.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
392 }
393
400 std::string fpSiA = PathResolver::find_file("protons_Si_Gunnar_Summers.dat" ,"DATAPATH");//E_kin < 15 MeV
401 std::string fpSiB = PathResolver::find_file("protons_Si_Gunnar_Huhtinen.dat" ,"DATAPATH");//E_kin > 15 MeV
402 std::string fnSiA = PathResolver::find_file("neutrons_Si_Gunnar_Griffin.dat" ,"DATAPATH");//E_kin < 20 MeV
403 std::string fnSiB = PathResolver::find_file("neutrons_Si_Gunnar_Konobeyev.dat","DATAPATH");//20 MeV < E_kin < 800 MeV
404 std::string fnSiC = PathResolver::find_file("neutrons_Si_Gunnar_Huhtinen.dat" ,"DATAPATH");//E_kin > 800 MeV
405 std::string fpiSi = PathResolver::find_file("pions_Si_Gunnar_Huhtinen.dat" ,"DATAPATH");//E_kin > 15 MeV
406 std::string feSi = PathResolver::find_file("electrons_Si_Summers.dat" ,"DATAPATH");//E_kin > 0.1 MeV
407
408 if ( !m_tgpSiA )
409 m_tgpSiA = new TGraph(fpSiA.c_str()); //E_kin < 15 MeV
410 if ( !m_tgpSiB )
411 m_tgpSiB = new TGraph(fpSiB.c_str()); //E_kin > 15 MeV
412 if ( !m_tgnSiA )
413 m_tgnSiA = new TGraph(fnSiA.c_str()); //E_kin < 20 MeV
414 if ( !m_tgnSiB )
415 m_tgnSiB = new TGraph(fnSiB.c_str()); //20 MeV < E_kin < 800 MeV
416 if ( !m_tgnSiC )
417 m_tgnSiC = new TGraph(fnSiC.c_str()); //E_kin > 800 MeV
418 if ( !m_tgpiSi )
419 m_tgpiSi = new TGraph(fpiSi.c_str()); //E_kin > 15 MeV
420 if ( !m_tgeSi )
421 m_tgeSi = new TGraph(feSi.c_str()); //E_kin > 0.1 MeV
422
426
427 std::string fHn = PathResolver::find_file("NeutronDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
428 std::string fHg = PathResolver::find_file("PhotonDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
429 std::string fHp = PathResolver::find_file("ProtonDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
430 std::string fHem = PathResolver::find_file("ElectronDose_in_Humans_ICRP_116.dat","DATAPATH");
431 std::string fHep = PathResolver::find_file("PositronDose_in_Humans_ICRP_116.dat","DATAPATH");
432 std::string fHmm = PathResolver::find_file("MuMinusDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
433 std::string fHmp = PathResolver::find_file("MuPlusDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
434 std::string fHpm = PathResolver::find_file("PiMinusDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
435 std::string fHpp = PathResolver::find_file("PiPlusDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
436 std::string fHa = PathResolver::find_file("AlphaDose_in_Humans_ICRP_116.dat" ,"DATAPATH");
437
438 if ( !m_tgHn )
439 m_tgHn = new TGraph(fHn.c_str() ,"%lg %*lg %*lg %*lg %*lg %*lg %lg");
440 if ( !m_tgHg )
441 m_tgHg = new TGraph(fHg.c_str() ,"%lg %*lg %*lg %*lg %*lg %*lg %lg");
442 if ( !m_tgHp )
443 m_tgHp = new TGraph(fHp.c_str() ,"%lg %*lg %*lg %*lg %*lg %*lg %lg");
444 if ( !m_tgHem )
445 m_tgHem = new TGraph(fHem.c_str(),"%lg %*lg %*lg %lg");
446 if ( !m_tgHep )
447 m_tgHep = new TGraph(fHep.c_str(),"%lg %*lg %*lg %lg");
448 if ( !m_tgHmm )
449 m_tgHmm = new TGraph(fHmm.c_str(),"%lg %*lg %*lg %lg");
450 if ( !m_tgHmp )
451 m_tgHmp = new TGraph(fHmp.c_str(),"%lg %*lg %*lg %lg");
452 if ( !m_tgHpm )
453 m_tgHpm = new TGraph(fHpm.c_str(),"%lg %*lg %*lg %lg");
454 if ( !m_tgHpp )
455 m_tgHpp = new TGraph(fHpp.c_str(),"%lg %*lg %*lg %lg");
456 if ( !m_tgHa )
457 m_tgHa = new TGraph(fHa.c_str() ,"%lg %*lg %*lg %lg");
458 }
459
460 void RadiationMapsMaker::UserSteppingAction(const G4Step* aStep){
461
462 bool goodMaterial(false);
463
464 if (m_config.materials.empty()) {
465 // if no material is specified maps are made for all materials
466 goodMaterial = true;
467 }
468 else {
469 // if a list of materials is specified, maps are made just for those.
470 // Note that still all steps need to be treated to calculate the
471 // volume fraction of the target materials in each bin!
472 for (unsigned int i=0;i<m_config.materials.size();i++) {
473 if ( m_config.materials[i].compare(aStep->GetTrack()->GetMaterial()->GetName()) == 0) {
474 goodMaterial = true;
475 break;
476 }
477 }
478 }
479
480 if ( goodMaterial && m_activationFile.is_open() ) {
481 // print list of unstable secondaries independent of status and energy (to see full chain)
482 const std::vector<const G4Track*>* tv = aStep->GetSecondaryInCurrentStep();
483 //const G4TrackVector * tv = aStep->GetSecondary();
484 //std::cout << "GetSecondary" << std::endl;
485 if ( tv ) {
486 //std::cout << "tv" << std::endl;
487 for (unsigned int is=0;is<tv->size();is++) {
488 //std::cout << "is " << is << std::endl;
489 const G4ParticleDefinition * pd = (*tv)[is]->GetParticleDefinition ();
490 if ( pd == G4Triton::Definition() || ( pd->GetParticleTable()->GetIonTable()->IsIon(pd) && !pd->GetPDGStable() ) ) {
491 G4String svname("unknown");
492 G4String svmaterial("unknown");
493 if ( (*tv)[is]->GetVolume() ) {
494 svname = (*tv)[is]->GetVolume()->GetName();
495 if ( (*tv)[is]->GetVolume()->GetLogicalVolume() ) {
496 if ( (*tv)[is]->GetVolume()->GetLogicalVolume()->GetMaterial() ) {
497 svmaterial = (*tv)[is]->GetVolume()->GetLogicalVolume()->GetMaterial()->GetName();
498 }
499 }
500 }
501 m_activationFile.width(18);
502 m_activationFile << (*tv)[is]->GetDefinition()->GetParticleName()
503 << " Proc.: " << (*tv)[is]->GetCreatorProcess()->GetProcessType() << " Mother: ";
504 m_activationFile.width(18);
505 m_activationFile << aStep->GetTrack()->GetDefinition()->GetParticleName() << " (x,y,z,t): (";
506 m_activationFile.setf ( std::ios::scientific );
507 m_activationFile.precision(4);
508 m_activationFile.width(12);
509 m_activationFile << aStep->GetPostStepPoint()->GetPosition().x() << ",";
510 m_activationFile.width(12);
511 m_activationFile << aStep->GetPostStepPoint()->GetPosition().y() << ",";
512 m_activationFile.width(12);
513 m_activationFile << aStep->GetPostStepPoint()->GetPosition().z() << ",";
514 m_activationFile.width(12);
515 m_activationFile << aStep->GetPostStepPoint()->GetGlobalTime()<< ")" ;
516 m_activationFile.unsetf ( std::ios::scientific );
517 m_activationFile << " Mat.: ";
518 m_activationFile.width(40);
519 m_activationFile << svmaterial << " Vol.: ";
520 m_activationFile << svname << std::endl;
521 }
522 }
523 }
524 }
525
526 int pdgid = 10;
527 if (aStep->GetTrack()->GetDefinition()==G4Gamma::Definition()){
528 pdgid=0;
529 } else if (aStep->GetTrack()->GetDefinition()==G4Proton::Definition()){
530 pdgid=1;
531 } else if (aStep->GetTrack()->GetDefinition()==G4PionPlus::Definition() ||
532 aStep->GetTrack()->GetDefinition()==G4PionMinus::Definition()){
533 pdgid=2;
534 } else if(aStep->GetTrack()->GetDefinition()==G4MuonPlus::Definition() ||
535 aStep->GetTrack()->GetDefinition()==G4MuonMinus::Definition()){
536 pdgid=3;
537 } else if(aStep->GetTrack()->GetDefinition()==G4Electron::Definition() ||
538 aStep->GetTrack()->GetDefinition()==G4Positron::Definition()){
539 if (aStep->GetTrack()->GetKineticEnergy()>0.5){
540 pdgid=4;
541 } else {
542 pdgid=5;
543 }
544 } else if(aStep->GetTrack()->GetDefinition()==G4Neutron::Definition()){
545 if (aStep->GetTrack()->GetKineticEnergy()>10.){
546 pdgid=6;
547 } else {
548 pdgid=7;
549 }
550 } else if (aStep->GetTrack()->GetDefinition()==G4KaonPlus::Definition() ||
551 aStep->GetTrack()->GetDefinition()==G4KaonMinus::Definition() ||
552 aStep->GetTrack()->GetDefinition()==G4KaonZeroShort::Definition() ||
553 aStep->GetTrack()->GetDefinition()==G4KaonZeroLong::Definition() ||
554 aStep->GetTrack()->GetDefinition()==G4KaonZero::Definition() ||
555 aStep->GetTrack()->GetDefinition()==G4AntiKaonZero::Definition() ||
556 aStep->GetTrack()->GetDefinition()==G4PionZero::Definition()){
557 pdgid=8; // particles not charged pions treated as charged pions for NIEL
558
559 } else if (aStep->GetTrack()->GetDefinition()==G4AntiProton::Definition() ||
560 aStep->GetTrack()->GetDefinition()==G4AntiNeutron::Definition() || // note n-bar is treated as proton like in FLUKA ... o.k. for > 10 MeV ...
561 aStep->GetTrack()->GetDefinition()==G4Lambda::Definition() ||
562 aStep->GetTrack()->GetDefinition()==G4AntiLambda::Definition() ||
563 aStep->GetTrack()->GetDefinition()==G4SigmaPlus::Definition() ||
564 aStep->GetTrack()->GetDefinition()==G4AntiSigmaPlus::Definition() ||
565 aStep->GetTrack()->GetDefinition()==G4SigmaMinus::Definition() ||
566 aStep->GetTrack()->GetDefinition()==G4AntiSigmaMinus::Definition() ||
567 aStep->GetTrack()->GetDefinition()==G4SigmaZero::Definition() ||
568 aStep->GetTrack()->GetDefinition()==G4AntiSigmaZero::Definition() ||
569 aStep->GetTrack()->GetDefinition()==G4XiMinus::Definition() ||
570 aStep->GetTrack()->GetDefinition()==G4AntiXiMinus::Definition() ||
571 aStep->GetTrack()->GetDefinition()==G4XiZero::Definition() ||
572 aStep->GetTrack()->GetDefinition()==G4AntiXiZero::Definition() ||
573 aStep->GetTrack()->GetDefinition()==G4OmegaMinus::Definition() ||
574 aStep->GetTrack()->GetDefinition()==G4AntiOmegaMinus::Definition() ||
575 aStep->GetTrack()->GetDefinition()==G4Alpha::Definition() ||
576 aStep->GetTrack()->GetDefinition()==G4Deuteron::Definition() ||
577 aStep->GetTrack()->GetDefinition()==G4Triton::Definition() ||
578 aStep->GetTrack()->GetDefinition()==G4He3::Definition()){
579 pdgid=9; // particles not protons treated as protons for NIEL
580
581 } else if (aStep->GetTrack()->GetDefinition()==G4Geantino::Definition()) {
582 pdgid = 999; // geantinos are used for volume calculation
583 } else if (!aStep->GetTrack()->GetDefinition()->GetPDGCharge()) return; // Not one of those and not a primary?
584
585 if ( pdgid == 10 && aStep->GetTrack()->GetKineticEnergy() > 10 ) {
586 // check Z for heavy ions with more than 10 MeV to also treat those for NIEL
587 if ( aStep->GetTrack()->GetDefinition()->GetAtomicNumber() > 1 ) {
588 pdgid = 9;
589 }
590 }
591
592 // process spectra, NIEL, h20, Edep, NEUT and CHAD particles only
593
594 if ( pdgid == 0 || pdgid == 3 || /* spectra */
595 pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 || /* NIEL & h20*/
596 aStep->GetTotalEnergyDeposit() > 0 || pdgid == 999) {
597
598 double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
599
600 double absq = std::fabs(charge);
601
602 double rho = aStep->GetTrack()->GetMaterial()->GetDensity()/CLHEP::g*CLHEP::cm3;
603
604 const G4Material * theMaterial = aStep->GetTrack()->GetMaterial();
605
606 // get time of step as average between pre- and post-step point
607 G4StepPoint* pre_step_point = aStep->GetPreStepPoint();
608 G4StepPoint* post_step_point = aStep->GetPostStepPoint();
609 const G4ThreeVector& startPoint = pre_step_point->GetPosition();
610 const G4ThreeVector& endPoint = post_step_point->GetPosition();
611 G4ThreeVector p = (startPoint + endPoint) * 0.5;
612
613 // process upper hemisphere only in case PositiveYOnly is true
614 if ( m_config.posYOnly && p.y() < 0 ) return;
615
616 double timeOfFlight = (pre_step_point->GetGlobalTime() +
617 post_step_point->GetGlobalTime()) * 0.5;
618
619 // then compute time difference to a particle traveling with speed of light until this position
620 double deltatime = (timeOfFlight - p.mag()/CLHEP::c_light)/CLHEP::s;
621
622 // and use the log10 of that time difference in seconds to fill time-dependent histograms
623 // note that the upper bin boundary is the relevant cut. The first bin contains thus all times even shorter than the first lower bin boundary
624 double logt = (deltatime<0?m_config.logTMin-1:log10(deltatime));
625
626 // fluence dN/da = dl/dV
627 double x0 = aStep->GetPreStepPoint()->GetPosition().x()*0.1;
628 double x1 = aStep->GetPostStepPoint()->GetPosition().x()*0.1;
629 double y0 = aStep->GetPreStepPoint()->GetPosition().y()*0.1;
630 double y1 = aStep->GetPostStepPoint()->GetPosition().y()*0.1;
631 double z0 = aStep->GetPreStepPoint()->GetPosition().z()*0.1;
632 double z1 = aStep->GetPostStepPoint()->GetPosition().z()*0.1;
633
634 double l = sqrt(pow(x1-x0,2)+pow(y1-y0,2)+pow(z1-z0,2));
635 // make 5 mm steps but at least 1 step
636 double dl0 = 0.5;
637 unsigned int nStep = l/dl0+1;
638 double dx = (x1-x0)/nStep;
639 double dy = (y1-y0)/nStep;
640 double dz = (z1-z0)/nStep;
641 double dl = l/nStep;
642
643 double weight = 0; // weight for NIEL
644 double wHt = 0; // weight for effecitve dose in humans
645 double eKin = aStep->GetTrack()->GetKineticEnergy();
646 double theta = aStep->GetTrack()->GetMomentumDirection().theta()*180./M_PI;
647 // if theta range in configuration is 0 - 90 assume that 180-theta should
648 // be used for theta > 90; otherwise leave theta unchanged
649 theta = (theta > 90&&m_config.thetaMin==0&&m_config.thetaMax==90?180-theta:theta);
650 double dphi = (aStep->GetTrack()->GetMomentumDirection().phi()-p.phi())*180./M_PI;
651 while ( dphi > 360 ) dphi -= 360.;
652 while ( dphi < 0 ) dphi += 360.;
653
654 double logEKin = (eKin > 0?log10(eKin):(m_config.logEMinn<m_config.logEMino?m_config.logEMinn:m_config.logEMino)-1);
655
656 if ( pdgid == 1 || pdgid == 9 ) {
657 // H_T
658 if ( pdgid == 1 ) {
659 // protons
660 wHt = m_tgHp->Eval(eKin);
661 }
662 else {
663 // baryons and heavy ions treat as alpha
664 wHt = m_tgHa->Eval(eKin);
665 }
666 if ( eKin < 15 ) {
667 if ( eKin > 10 ) {
668 weight = m_tgpSiA->Eval(eKin);
669 }
670 }
671 else {
672 weight = m_tgpSiB->Eval(eKin);
673 }
674 }
675 else if ( pdgid == 2 || pdgid == 8 ) {
676 // H_T
677 // pions and light mesons
678 if ( charge < 0 ) {
679 wHt = m_tgHpm->Eval(eKin);
680 }
681 else {
682 wHt = m_tgHpp->Eval(eKin);
683 }
684 if ( eKin > 15 ) {
685 weight = m_tgpiSi->Eval(eKin);
686 }
687 }
688 else if ( pdgid == 6 || pdgid == 7 ) {
689 // H_T
690 // neutrons
691 wHt = m_tgHn->Eval(eKin);
692 if ( eKin < 20 ) {
693 weight = m_tgnSiA->Eval(eKin);
694 }
695 else if ( eKin < 800 ) {
696 weight = m_tgnSiB->Eval(eKin);
697 }
698 else {
699 weight = m_tgnSiC->Eval(eKin);
700 }
701 }
702 else if ( pdgid == 4 || pdgid == 5 ) {
703 // H_T
704 // electrons and positrons
705 if ( charge < 0 ) {
706 wHt = m_tgHem->Eval(eKin);
707 }
708 else {
709 wHt = m_tgHep->Eval(eKin);
710 }
711 if ( eKin > 0.1 ) {
712 weight = m_tgeSi->Eval(eKin);
713 }
714 }
715 else if ( pdgid == 0 ) {
716 // H_T
717 // photons
718 wHt = m_tgHg->Eval(eKin);
719 }
720 else if ( pdgid == 3 ) {
721 // H_T
722 // muons
723 if ( charge < 0 ) {
724 wHt = m_tgHmm->Eval(eKin);
725 }
726 else {
727 wHt = m_tgHmp->Eval(eKin);
728 }
729 }
730
731 double dE_TOT = aStep->GetTotalEnergyDeposit()/nStep;
732 double dE_NIEL = aStep->GetNonIonizingEnergyDeposit()/nStep;
733 double dE_ION = dE_TOT-dE_NIEL;
734 double dH_T = 1e-12*wHt; // H_T convert from pSv to Sv
735
736 for(unsigned int i=0;i<nStep;i++) {
737 // randomize position along current sub-step (flat fraction from 0-1)
738 double subpos = G4UniformRand();
739 double abszorz = z0+dz*(i+subpos);
740 // if |z| instead of z take abs value
741 if ( m_config.zMinFull >= 0 ) abszorz = std::fabs(abszorz);
742
743 double rr = sqrt(pow(x0+dx*(i+subpos),2)+
744 pow(y0+dy*(i+subpos),2));
745 double pphi = atan2(y0+dy*(i+subpos),x0+dx*(i+subpos))*180/M_PI;
746
747 int vBinZoom = -1;
748 int vBinFull = -1;
749 int vBin3d = -1;
750 int vBinZoomSpecn = -1;
751 int vBinFullSpecn = -1;
752 int vBinZoomSpeco = -1;
753 int vBinFullSpeco = -1;
754 int vBinZoomTime = -1;
755 int vBinFullTime = -1;
756 int vBinThetaFullSpecn = -1;
757 int vBinThetaFullSpeco = -1;
758
759 // zoom 2d
760 if ( m_config.zMinZoom < abszorz &&
761 m_config.zMaxZoom > abszorz ) {
762 int iz = (abszorz-m_config.zMinZoom)/(m_config.zMaxZoom-m_config.zMinZoom)*m_config.nBinsz;
763 if ( m_config.rMinZoom < rr &&
764 m_config.rMaxZoom > rr ) {
765 int ir = (rr-m_config.rMinZoom)/(m_config.rMaxZoom-m_config.rMinZoom)*m_config.nBinsr;
766 vBinZoom = m_config.nBinsr*iz+ir;
767 if ( m_config.logTMax > logt ){
768 int ilt = (logt < m_config.logTMin?0:(logt-m_config.logTMin)/(m_config.logTMax-m_config.logTMin)*m_config.nBinslogT);
769 vBinZoomTime = m_config.nBinsr*m_config.nBinslogT*iz+ir*m_config.nBinslogT+ilt;
770 }
771 if ( m_config.logEMinn < logEKin &&
772 m_config.logEMaxn > logEKin &&
773 (pdgid == 6 || pdgid == 7)) {
774 int ile = (logEKin-m_config.logEMinn)/(m_config.logEMaxn-m_config.logEMinn)*m_config.nBinslogEn;
775 vBinZoomSpecn = m_config.nBinsr*m_config.nBinslogEn*iz+ir*m_config.nBinslogEn+ile;
776 }
777 if ( m_config.logEMino < logEKin &&
778 m_config.logEMaxo > logEKin &&
779 pdgid != 6 && pdgid != 7) {
780 int ile = (logEKin-m_config.logEMino)/(m_config.logEMaxo-m_config.logEMino)*m_config.nBinslogEo;
781 vBinZoomSpeco = m_config.nBinsr*m_config.nBinslogEo*iz+ir*m_config.nBinslogEo+ile;
782 }
783 }
784 }
785
786 // full 2d
787 if ( m_config.zMinFull < abszorz &&
788 m_config.zMaxFull > abszorz ) {
789 int iz = (abszorz-m_config.zMinFull)/(m_config.zMaxFull-m_config.zMinFull)*m_config.nBinsz;
790 if ( m_config.rMinFull < rr &&
791 m_config.rMaxFull > rr ) {
792 int ir = (rr-m_config.rMinFull)/(m_config.rMaxFull-m_config.rMinFull)*m_config.nBinsr;
793 vBinFull = m_config.nBinsr*iz+ir;
794 if ( m_config.logTMax > logt ){
795 int ilt = (logt < m_config.logTMin?0:(logt-m_config.logTMin)/(m_config.logTMax-m_config.logTMin)*m_config.nBinslogT);
796 vBinFullTime = m_config.nBinsr*m_config.nBinslogT*iz+ir*m_config.nBinslogT+ilt;
797 }
798 if ( m_config.logEMinn < logEKin &&
799 m_config.logEMaxn > logEKin &&
800 (pdgid == 6 || pdgid == 7)) {
801 int ile = (logEKin-m_config.logEMinn)/(m_config.logEMaxn-m_config.logEMinn)*m_config.nBinslogEn;
802 vBinFullSpecn = m_config.nBinsr*m_config.nBinslogEn*iz+ir*m_config.nBinslogEn+ile;
803 if ( m_config.thetaMin < theta &&
804 m_config.thetaMax > theta ) {
805 int ith = (theta-m_config.thetaMin)/(m_config.thetaMax-m_config.thetaMin)*m_config.nBinstheta;
806 int idph = dphi/360.*m_config.nBinsdphi;
807 ith = ith*m_config.nBinsdphi+idph;
808 vBinThetaFullSpecn = m_config.nBinsr*m_config.nBinslogEn*m_config.nBinsdphi*m_config.nBinstheta*iz+ir*m_config.nBinslogEn*m_config.nBinsdphi*m_config.nBinstheta+ile*m_config.nBinsdphi*m_config.nBinstheta+ith;
809 }
810 }
811 if ( m_config.logEMino < logEKin &&
812 m_config.logEMaxo > logEKin &&
813 pdgid != 6 && pdgid != 7) {
814 int ile = (logEKin-m_config.logEMino)/(m_config.logEMaxo-m_config.logEMino)*m_config.nBinslogEo;
815 vBinFullSpeco = m_config.nBinsr*m_config.nBinslogEo*iz+ir*m_config.nBinslogEo+ile;
816 if ( m_config.thetaMin < theta &&
817 m_config.thetaMax > theta ) {
818 int ith = (theta-m_config.thetaMin)/(m_config.thetaMax-m_config.thetaMin)*m_config.nBinstheta;
819 int idph = dphi/360.*m_config.nBinsdphi;
820 ith = ith*m_config.nBinsdphi+idph;
821 vBinThetaFullSpeco = m_config.nBinsr*m_config.nBinslogEo*m_config.nBinsdphi*m_config.nBinstheta*iz+ir*m_config.nBinslogEo*m_config.nBinsdphi*m_config.nBinstheta+ile*m_config.nBinsdphi*m_config.nBinstheta+ith;
822 }
823 }
824 }
825 }
826
827 // zoom 3d
828 if ( m_config.zMinZoom < abszorz &&
829 m_config.zMaxZoom > abszorz ) {
830 int iz = (abszorz-m_config.zMinZoom)/(m_config.zMaxZoom-m_config.zMinZoom)*m_config.nBinsz3d;
831 if ( m_config.rMinZoom < rr &&
832 m_config.rMaxZoom > rr ) {
833 int ir = (rr-m_config.rMinZoom)/(m_config.rMaxZoom-m_config.rMinZoom)*m_config.nBinsr3d;
834 if ( m_config.phiMinZoom == 0) {
835 // assume that all phi should be mapped to the selected phi range
836 double phi_mapped = pphi;
837 // first use phi range from 0 - 360 degrees
838 if (phi_mapped < 0) {
839 phi_mapped = 360 + phi_mapped;
840 }
841 // then map to selected phi-range
842 int iphi = phi_mapped/m_config.phiMaxZoom;
843 phi_mapped -= iphi*m_config.phiMaxZoom;
844 iphi = phi_mapped/m_config.phiMaxZoom*m_config.nBinsphi3d;
845 vBin3d = m_config.nBinsr3d*m_config.nBinsphi3d*iz+m_config.nBinsphi3d*ir+iphi;
846 }
847 else if ( m_config.phiMinZoom < pphi &&
848 m_config.phiMaxZoom > pphi ) {
849 int iphi = (pphi-m_config.phiMinZoom)/(m_config.phiMaxZoom-m_config.phiMinZoom)*m_config.nBinsphi3d;
850 vBin3d = m_config.nBinsr3d*m_config.nBinsphi3d*iz+m_config.nBinsphi3d*ir+iphi;
851 }
852 }
853 }
854 // TID & EION
855 if ( goodMaterial && vBinZoom >=0 ) {
856 if ( pdgid == 999 ) {
857 m_maps.m_rz_tid [vBinZoom] += rr*dl;
858 m_maps.m_rz_eion[vBinZoom] += rho*rr*dl;
859 for (size_t ie=0;ie<theMaterial->GetNumberOfElements();ie++ ) {
860 const G4Element * theElement = theMaterial->GetElement (ie);
861 const double mFrac = theMaterial->GetFractionVector()[ie];
862 const int zElem = theElement->GetZ();
863 if ( zElem >= m_config.elemZMin &&
864 zElem <= m_config.elemZMax) {
865 m_maps.m_rz_element[vBinZoom*(m_config.elemZMax-m_config.elemZMin+1)+zElem-m_config.elemZMin] += rho*rr*dl*mFrac;
866 }
867 }
868 }
869 else {
870 m_maps.m_rz_tid [vBinZoom] += dE_ION/rho;
871 m_maps.m_rz_eion[vBinZoom] += dE_ION;
872 for (size_t ie=0;ie<theMaterial->GetNumberOfElements();ie++ ) {
873 const G4Element * theElement = theMaterial->GetElement (ie);
874 const double mFrac = theMaterial->GetFractionVector()[ie];
875 const int zElem = theElement->GetZ();
876 if ( zElem >= m_config.elemZMin &&
877 zElem <= m_config.elemZMax) {
878 m_maps.m_rz_element[vBinZoom*(m_config.elemZMax-m_config.elemZMin+1)+zElem-m_config.elemZMin] += dE_ION*mFrac;
879 }
880 }
881 }
882 }
883 if ( goodMaterial && vBinZoomTime >=0 ) {
884 m_maps.m_rz_tid_time[vBinZoomTime] += dE_ION/rho;
885 if ( dH_T > 0 ) {
886 m_maps.m_rz_ht_time[vBinZoomTime] += dH_T*dl;
887 }
888 }
889 if ( goodMaterial && vBinFull >=0 ) {
890 if ( pdgid == 999 ) {
891 m_maps.m_full_rz_tid [vBinFull] += rr*dl;
892 m_maps.m_full_rz_eion[vBinFull] += rho*rr*dl;
893 for (size_t ie=0;ie<theMaterial->GetNumberOfElements();ie++ ) {
894 const G4Element * theElement = theMaterial->GetElement (ie);
895 const double mFrac = theMaterial->GetFractionVector()[ie];
896 const int zElem = theElement->GetZ();
897 if ( zElem >= m_config.elemZMin &&
898 zElem <= m_config.elemZMax) {
899 m_maps.m_full_rz_element[vBinFull*(m_config.elemZMax-m_config.elemZMin+1)+zElem-m_config.elemZMin] += rho*rr*dl*mFrac;
900 }
901 }
902 }
903 else {
904 m_maps.m_full_rz_tid [vBinFull] += dE_ION/rho;
905 m_maps.m_full_rz_eion[vBinFull] += dE_ION;
906 for (size_t ie=0;ie<theMaterial->GetNumberOfElements();ie++ ) {
907 const G4Element * theElement = theMaterial->GetElement (ie);
908 const double mFrac = theMaterial->GetFractionVector()[ie];
909 const int zElem = theElement->GetZ();
910 if ( zElem >= m_config.elemZMin &&
911 zElem <= m_config.elemZMax) {
912 m_maps.m_full_rz_element[vBinFull*(m_config.elemZMax-m_config.elemZMin+1)+zElem-m_config.elemZMin] += dE_ION*mFrac;
913 }
914 }
915 }
916 }
917 if ( goodMaterial && vBinFullTime >=0 ) {
918 m_maps.m_full_rz_tid_time[vBinFullTime] += dE_ION/rho;
919 if ( dH_T > 0 ) {
920 m_maps.m_full_rz_ht_time[vBinFullTime] += dH_T*dl;
921 }
922 }
923 if ( goodMaterial && vBin3d >=0 ) {
924 if ( pdgid == 999 ) {
925 m_maps.m_3d_tid [vBin3d] += rr*dl;
926 m_maps.m_3d_eion[vBin3d] += rho*rr*dl;
927 }
928 else {
929 m_maps.m_3d_tid [vBin3d] += dE_ION/rho;
930 m_maps.m_3d_eion[vBin3d] += dE_ION;
931 }
932 }
933
934 if ( goodMaterial && (pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 )) {
935 // NIEL
936 if ( weight > 0 ) {
937 if ( vBinZoom >=0 ) {
938 m_maps.m_rz_niel [vBinZoom] += weight*dl;
939 }
940 if ( vBinFull >=0 ) {
941 m_maps.m_full_rz_niel [vBinFull] += weight*dl;
942 }
943 if ( vBin3d >=0 ) {
944 m_maps.m_3d_niel [vBin3d] += weight*dl;
945 }
946 }
947 // SEE
948 if ( eKin > 20 && (pdgid == 1 || pdgid == 2 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9) ) {
949 if ( vBinZoom >=0 ) {
950 m_maps.m_rz_h20 [vBinZoom] += dl;
951 }
952 if ( vBinFull >=0 ) {
953 m_maps.m_full_rz_h20 [vBinFull] += dl;
954 }
955 if ( vBin3d >=0 ) {
956 m_maps.m_3d_h20 [vBin3d] += dl;
957 }
958 }
959 // NEUT
960 if ( eKin > 0.1 && (pdgid == 6 || pdgid == 7 ) ) {
961 if ( vBinZoom >=0 ) {
962 m_maps.m_rz_neut [vBinZoom] += dl;
963 }
964 if ( vBinFull >=0 ) {
965 m_maps.m_full_rz_neut [vBinFull] += dl;
966 }
967 if ( vBin3d >=0 ) {
968 m_maps.m_3d_neut [vBin3d] += dl;
969 }
970 }
971 // CHAD
972 if ( absq > 0 && (pdgid == 1 || pdgid == 2 || pdgid == 8 || pdgid == 9 ) ) {
973 if ( vBinZoom >=0 ) {
974 m_maps.m_rz_chad [vBinZoom] += dl;
975 }
976 if ( vBinFull >=0 ) {
977 m_maps.m_full_rz_chad [vBinFull] += dl;
978 }
979 if ( vBin3d >=0 ) {
980 m_maps.m_3d_chad [vBin3d] += dl;
981 }
982 }
983 // Neutron Energy Spectra
984 if ( pdgid == 6 || pdgid == 7 ) {
985 if ( vBinZoomSpecn >=0 ) {
986 m_maps.m_rz_neut_spec [vBinZoomSpecn] += dl;
987 }
988 if ( vBinFullSpecn >=0 ) {
989 m_maps.m_full_rz_neut_spec [vBinFullSpecn] += dl;
990 }
991 if ( vBinThetaFullSpecn >=0 ) {
992 m_maps.m_theta_full_rz_neut_spec [vBinThetaFullSpecn] += dl;
993 }
994 }
995 }
996
997 if ( goodMaterial && (pdgid < 6 || pdgid > 7 )){
998 // Other particle Energy Spectra
999 if ( vBinZoomSpeco >=0 ) {
1000 if ( pdgid == 0 ) {
1001 m_maps.m_rz_gamm_spec [vBinZoomSpeco] += dl;
1002 }
1003 else if ( pdgid == 1 ) {
1004 m_maps.m_rz_prot_spec [vBinZoomSpeco] += dl;
1005 }
1006 else if ( pdgid == 2 ) {
1007 m_maps.m_rz_pion_spec [vBinZoomSpeco] += dl;
1008 }
1009 else if ( pdgid == 3 ) {
1010 m_maps.m_rz_muon_spec [vBinZoomSpeco] += dl;
1011 }
1012 else if ( pdgid == 4 || pdgid == 5 ) {
1013 m_maps.m_rz_elec_spec [vBinZoomSpeco] += dl;
1014 }
1015 else {
1016 m_maps.m_rz_rest_spec [vBinZoomSpeco] += dl;
1017 }
1018 }
1019 if ( vBinFullSpeco >=0 ) {
1020 if ( pdgid == 0 ) {
1021 m_maps.m_full_rz_gamm_spec [vBinFullSpeco] += dl;
1022 }
1023 else if ( pdgid == 1 ) {
1024 m_maps.m_full_rz_prot_spec [vBinFullSpeco] += dl;
1025 }
1026 else if ( pdgid == 2 ) {
1027 m_maps.m_full_rz_pion_spec [vBinFullSpeco] += dl;
1028 }
1029 else if ( pdgid == 3 ) {
1030 m_maps.m_full_rz_muon_spec [vBinFullSpeco] += dl;
1031 }
1032 else if ( pdgid == 4 || pdgid == 5 ) {
1033 m_maps.m_full_rz_elec_spec [vBinFullSpeco] += dl;
1034 }
1035 else {
1036 m_maps.m_full_rz_rest_spec [vBinFullSpeco] += dl;
1037 }
1038 }
1039 if ( vBinThetaFullSpeco >=0 ) {
1040 if ( pdgid == 0 ) {
1041 m_maps.m_theta_full_rz_gamm_spec [vBinThetaFullSpeco] += dl;
1042 }
1043 else if ( pdgid == 1 ) {
1044 m_maps.m_theta_full_rz_prot_spec [vBinThetaFullSpeco] += dl;
1045 }
1046 else if ( pdgid == 2 ) {
1047 m_maps.m_theta_full_rz_pion_spec [vBinThetaFullSpeco] += dl;
1048 }
1049 else if ( pdgid == 3 ) {
1050 m_maps.m_theta_full_rz_muon_spec [vBinThetaFullSpeco] += dl;
1051 }
1052 else if ( pdgid == 4 || pdgid == 5 ) {
1053 m_maps.m_theta_full_rz_elec_spec [vBinThetaFullSpeco] += dl;
1054 }
1055 else {
1056 if ( absq > 0 ) {
1057 m_maps.m_theta_full_rz_rchgd_spec [vBinThetaFullSpeco] += dl;
1058 }
1059 else {
1060 m_maps.m_theta_full_rz_rneut_spec [vBinThetaFullSpeco] += dl;
1061 }
1062 }
1063 }
1064 }
1065
1066 if (!m_config.materials.empty()) {
1067 // need volume fraction only if particular materials are selected
1068 if ( (eKin > 1 && (pdgid == 6 || pdgid == 7)) || pdgid == 999) {
1069 // count all neutron > 1 MeV track lengths weighted by r
1070 // to get norm for volume per bin. High energetic
1071 // neutrons are used because they travel far enough to
1072 // map entire bins and are not bent by magnetic fields.
1073 // dl is a measure of length inside the current bin.
1074 // The multiplication by r accounts for the larger
1075 // volume corresponding to larger r assuming that the
1076 // neutron flux is locally mainly from inside to
1077 // outside. In regions where the neutron flux differs
1078 // substantially from this cylindrical assumption a
1079 // cylindrical Geantino scan (vertex: flat in z, x=y=0;
1080 // momentum: pz=0, flat in phi) should be used to get
1081 // the correct volume fraction.
1082 if ( vBinZoom >=0 ) {
1083 m_maps.m_rz_norm[vBinZoom] += rr*dl;
1084 }
1085 if ( vBinFull >=0 ) {
1086 m_maps.m_full_rz_norm[vBinFull] += rr*dl;
1087 }
1088 if ( vBin3d >=0 ) {
1089 m_maps.m_3d_norm[vBin3d] += rr*dl;
1090 }
1091 if ( goodMaterial ) {
1092 // same but only inside the materials of interest.
1093 // The ratio vol/norm gives the volume fraction of the desired
1094 // materials inside the current bin
1095 if ( vBinZoom >=0 ) {
1096 m_maps.m_rz_vol [vBinZoom] += rr*dl;
1097 }
1098 if ( vBinFull >=0 ) {
1099 m_maps.m_full_rz_vol [vBinFull] += rr*dl;
1100 }
1101 if ( vBin3d >=0 ) {
1102 m_maps.m_3d_vol [vBin3d] += rr*dl;
1103 }
1104 }
1105 }
1106 }
1107 }
1108 }
1109 }
1110} // namespace G4UA
1111
rr
const boost::regex rr(r_r)
M_PI
#define M_PI
Definition ActiveFraction.h:11
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
charge
double charge(const T &p)
Definition AtlasPID.h:997
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
G4UA::RadiationMapsMaker::BeginOfRunAction
virtual void BeginOfRunAction(const G4Run *) override final
Definition RadiationMapsMaker.cxx:255
G4UA::RadiationMapsMaker::RadiationMapsMaker
RadiationMapsMaker(const Config &config)
Definition RadiationMapsMaker.cxx:62
G4UA::RadiationMapsMaker::UserSteppingAction
virtual void UserSteppingAction(const G4Step *) override final
Definition RadiationMapsMaker.cxx:460
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path)
Definition PathResolver.cxx:220
ir
int ir
counter of the current depth
Definition fastadd.cxx:49
G4UA::RadiationMapsMaker::Report
Simple struct for holding the radiation maps.
Definition RadiationMapsMaker.h:91
G4UA::RadiationMapsMaker::Report::m_3d_eion
std::vector< double > m_3d_eion
vector of ionizing energy density seen by thread in 3d
Definition RadiationMapsMaker.h:135
G4UA::RadiationMapsMaker::Report::m_full_rz_tid
std::vector< double > m_full_rz_tid
vector of tid seen by thread in full area
Definition RadiationMapsMaker.h:113
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_rneut_spec
std::vector< double > m_theta_full_rz_rneut_spec
vector of rest neutral spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:217
G4UA::RadiationMapsMaker::Report::m_3d_h20
std::vector< double > m_3d_h20
vector of >20 MeV hadron flux seen by thread in 3d
Definition RadiationMapsMaker.h:139
G4UA::RadiationMapsMaker::Report::m_rz_neut_spec
std::vector< double > m_rz_neut_spec
vector of neutron spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:157
G4UA::RadiationMapsMaker::Report::m_3d_chad
std::vector< double > m_3d_chad
vector of charged hadron flux seen by thread in 3d
Definition RadiationMapsMaker.h:143
G4UA::RadiationMapsMaker::Report::m_rz_vol
std::vector< double > m_rz_vol
next two vectors are used only in case maps are needed for particular materials instead of all
Definition RadiationMapsMaker.h:108
G4UA::RadiationMapsMaker::Report::m_full_rz_muon_spec
std::vector< double > m_full_rz_muon_spec
vector of mu^+/- spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:186
G4UA::RadiationMapsMaker::Report::m_3d_neut
std::vector< double > m_3d_neut
vector of >100 keV hadron flux seen by thread in 3d
Definition RadiationMapsMaker.h:141
G4UA::RadiationMapsMaker::Report::m_full_rz_chad
std::vector< double > m_full_rz_chad
vector of charged hadron flux seen by thread in full area
Definition RadiationMapsMaker.h:123
G4UA::RadiationMapsMaker::Report::m_3d_norm
std::vector< double > m_3d_norm
vector to normalize the volume fraction in 3d
Definition RadiationMapsMaker.h:150
G4UA::RadiationMapsMaker::Report::m_full_rz_tid_time
std::vector< double > m_full_rz_tid_time
vector of time dependent TID in full 2d grid
Definition RadiationMapsMaker.h:224
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_elec_spec
std::vector< double > m_theta_full_rz_elec_spec
vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:179
G4UA::RadiationMapsMaker::Report::m_rz_chad
std::vector< double > m_rz_chad
vector of charged hadron flux seen by thread in zoomed area
Definition RadiationMapsMaker.h:103
G4UA::RadiationMapsMaker::Report::m_3d_niel
std::vector< double > m_3d_niel
vector of 1 MeV neutron equivalent flux seen by thread in 3d
Definition RadiationMapsMaker.h:137
G4UA::RadiationMapsMaker::Report::m_3d_tid
std::vector< double > m_3d_tid
vector of tid seen by thread in 3d
Definition RadiationMapsMaker.h:133
G4UA::RadiationMapsMaker::Report::m_rz_h20
std::vector< double > m_rz_h20
vector of >20 MeV hadron flux seen by thread in zoomed area
Definition RadiationMapsMaker.h:99
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_gamm_spec
std::vector< double > m_theta_full_rz_gamm_spec
vector of gamma spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:170
G4UA::RadiationMapsMaker::Report::m_full_rz_element
std::vector< double > m_full_rz_element
vector of element fractions in full 2d grid
Definition RadiationMapsMaker.h:236
G4UA::RadiationMapsMaker::Report::m_rz_muon_spec
std::vector< double > m_rz_muon_spec
vector of mu^+/- spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:184
G4UA::RadiationMapsMaker::Report::m_full_rz_pion_spec
std::vector< double > m_full_rz_pion_spec
vector of pi^+/- spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:195
G4UA::RadiationMapsMaker::Report::m_full_rz_elec_spec
std::vector< double > m_full_rz_elec_spec
vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:177
G4UA::RadiationMapsMaker::Report::m_rz_tid
std::vector< double > m_rz_tid
vector of tid seen by thread in zoomed area
Definition RadiationMapsMaker.h:93
G4UA::RadiationMapsMaker::Report::m_full_rz_niel
std::vector< double > m_full_rz_niel
vector of 1 MeV neutron equivalent flux seen by thread in full area
Definition RadiationMapsMaker.h:117
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_muon_spec
std::vector< double > m_theta_full_rz_muon_spec
vector of mu^+/- spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:188
G4UA::RadiationMapsMaker::Report::m_rz_tid_time
std::vector< double > m_rz_tid_time
vector of time dependent TID in zoom 2d grid
Definition RadiationMapsMaker.h:222
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_prot_spec
std::vector< double > m_theta_full_rz_prot_spec
vector of proton spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:206
G4UA::RadiationMapsMaker::Report::m_rz_gamm_spec
std::vector< double > m_rz_gamm_spec
vector of gamma spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:166
G4UA::RadiationMapsMaker::Report::m_full_rz_gamm_spec
std::vector< double > m_full_rz_gamm_spec
vector of gamma spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:168
G4UA::RadiationMapsMaker::Report::m_full_rz_vol
std::vector< double > m_full_rz_vol
next two vectors are used only in case maps are needed for particular materials instead of all
Definition RadiationMapsMaker.h:128
G4UA::RadiationMapsMaker::Report::m_rz_norm
std::vector< double > m_rz_norm
vector to normalize the volume fraction in zoomed area
Definition RadiationMapsMaker.h:110
G4UA::RadiationMapsMaker::Report::m_rz_ht_time
std::vector< double > m_rz_ht_time
vector of time dependent H_T in zoom 2d grid
Definition RadiationMapsMaker.h:227
G4UA::RadiationMapsMaker::Report::m_full_rz_h20
std::vector< double > m_full_rz_h20
vector of >20 MeV hadron flux seen by thread in full area
Definition RadiationMapsMaker.h:119
G4UA::RadiationMapsMaker::Report::m_full_rz_norm
std::vector< double > m_full_rz_norm
vector to normalize the volume fraction in full area
Definition RadiationMapsMaker.h:130
G4UA::RadiationMapsMaker::Report::m_rz_prot_spec
std::vector< double > m_rz_prot_spec
vector of proton spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:202
G4UA::RadiationMapsMaker::Report::m_3d_vol
std::vector< double > m_3d_vol
next two vectors are used only in case maps are needed for particular materials instead of all
Definition RadiationMapsMaker.h:148
G4UA::RadiationMapsMaker::Report::m_rz_niel
std::vector< double > m_rz_niel
vector of 1 MeV neutron equivalent flux seen by thread in zoomed area
Definition RadiationMapsMaker.h:97
G4UA::RadiationMapsMaker::Report::m_full_rz_rest_spec
std::vector< double > m_full_rz_rest_spec
vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:213
G4UA::RadiationMapsMaker::Report::m_rz_neut
std::vector< double > m_rz_neut
vector of >100 keV hadron flux seen by thread in zoomed area
Definition RadiationMapsMaker.h:101
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_pion_spec
std::vector< double > m_theta_full_rz_pion_spec
vector of pi^+/- spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:197
G4UA::RadiationMapsMaker::Report::m_rz_pion_spec
std::vector< double > m_rz_pion_spec
vector of pi^+/- spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:193
G4UA::RadiationMapsMaker::Report::m_full_rz_neut
std::vector< double > m_full_rz_neut
vector of >100 keV hadron flux seen by thread in full area
Definition RadiationMapsMaker.h:121
G4UA::RadiationMapsMaker::Report::m_rz_eion
std::vector< double > m_rz_eion
vector of ionizing energy density seen by thread in zoomed area
Definition RadiationMapsMaker.h:95
G4UA::RadiationMapsMaker::Report::m_rz_elec_spec
std::vector< double > m_rz_elec_spec
vector of e^+/- spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:175
G4UA::RadiationMapsMaker::Report::m_full_rz_ht_time
std::vector< double > m_full_rz_ht_time
vector of time dependent H_T in full 2d grid
Definition RadiationMapsMaker.h:229
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_rchgd_spec
std::vector< double > m_theta_full_rz_rchgd_spec
vector of rest charged spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:215
G4UA::RadiationMapsMaker::Report::m_rz_rest_spec
std::vector< double > m_rz_rest_spec
vector of other particle spectra in log10(E/MeV) bins and the zoom 2d grid
Definition RadiationMapsMaker.h:211
G4UA::RadiationMapsMaker::Report::m_full_rz_prot_spec
std::vector< double > m_full_rz_prot_spec
vector of proton spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:204
G4UA::RadiationMapsMaker::Report::m_full_rz_eion
std::vector< double > m_full_rz_eion
vector of ionizing energy density seen by thread in full area
Definition RadiationMapsMaker.h:115
G4UA::RadiationMapsMaker::Report::m_full_rz_neut_spec
std::vector< double > m_full_rz_neut_spec
vector of neutron spectra in log10(E/MeV) bins and the full 2d grid
Definition RadiationMapsMaker.h:159
G4UA::RadiationMapsMaker::Report::m_rz_element
std::vector< double > m_rz_element
vector of element fractions in zoom 2d grid
Definition RadiationMapsMaker.h:234
G4UA::RadiationMapsMaker::Report::m_theta_full_rz_neut_spec
std::vector< double > m_theta_full_rz_neut_spec
vector of neutron spectra in log10(E/MeV) bins and the full 2d grid x theta bins
Definition RadiationMapsMaker.h:161
Generated on for ATLAS Offline Software by doxygen 1.14.0