ATLAS Offline Software
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"
55 
56 namespace G4UA{
57 
58 
59  //----------------------------------------------------------------------------
60  // Constructor
61  //----------------------------------------------------------------------------
63  : m_config(config)
64  {
65  }
66 
67  //----------------------------------------------------------------------------
68  // Destructor
69  //----------------------------------------------------------------------------
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  //---------------------------------------------------------------------------
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
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];
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];
223  }
224  // x theta bins
225  for(unsigned int i=0;i<maps.m_theta_full_rz_neut_spec.size();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];
232  m_rz_elec_spec [i] += maps.m_rz_elec_spec [i];
234  m_rz_muon_spec [i] += maps.m_rz_muon_spec [i];
236  m_rz_pion_spec [i] += maps.m_rz_pion_spec [i];
238  m_rz_prot_spec [i] += maps.m_rz_prot_spec [i];
240  m_rz_rest_spec [i] += maps.m_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++) {
252  }
253  }
254 
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 
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 
334 
341 
348 
353 
356 
371 
380 
381  if (!m_config.materials.empty()) {
382  // need volume fraction only if particular materials are selected
383  // 2d zoom
386  // 2d full
389  // 3d
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 ) {
763  if ( m_config.rMinZoom < rr &&
764  m_config.rMaxZoom > rr ) {
766  vBinZoom = m_config.nBinsr*iz+ir;
767  if ( m_config.logTMax > logt ){
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)) {
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) {
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 ) {
790  if ( m_config.rMinFull < rr &&
791  m_config.rMaxFull > rr ) {
793  vBinFull = m_config.nBinsr*iz+ir;
794  if ( m_config.logTMax > logt ){
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)) {
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 ) {
806  int idph = dphi/360.*m_config.nBinsdphi;
807  ith = ith*m_config.nBinsdphi+idph;
809  }
810  }
811  if ( m_config.logEMino < logEKin &&
812  m_config.logEMaxo > logEKin &&
813  pdgid != 6 && pdgid != 7) {
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 ) {
819  int idph = dphi/360.*m_config.nBinsdphi;
820  ith = ith*m_config.nBinsdphi+idph;
822  }
823  }
824  }
825  }
826 
827  // zoom 3d
828  if ( m_config.zMinZoom < abszorz &&
829  m_config.zMaxZoom > abszorz ) {
831  if ( m_config.rMinZoom < rr &&
832  m_config.rMaxZoom > rr ) {
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;
846  }
847  else if ( m_config.phiMinZoom < pphi &&
848  m_config.phiMaxZoom > pphi ) {
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) {
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) {
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 
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_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_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::Config::nBinslogT
int nBinslogT
Definition: RadiationMapsMaker.h:79
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_element
std::vector< double > m_rz_element
vector of element fractions in zoom 2d grid
Definition: RadiationMapsMaker.h:234
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
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
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::m_tgeSi
TGraph * m_tgeSi
Definition: RadiationMapsMaker.h:267
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::m_activationFile
std::ofstream m_activationFile
Definition: RadiationMapsMaker.h:280
TestSUSYToolsAlg.dl
dl
Definition: TestSUSYToolsAlg.py:83
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
G4UA::RadiationMapsMaker::m_tgHa
TGraph * m_tgHa
Definition: RadiationMapsMaker.h:278
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
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
G4UA::RadiationMapsMaker::Config::logEMino
double logEMino
Definition: RadiationMapsMaker.h:75
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path, SearchType search_type=LocalSearch)
Definition: PathResolver.cxx:251
G4UA::RadiationMapsMaker::m_tgHmp
TGraph * m_tgHmp
Definition: RadiationMapsMaker.h:275
G4UA::RadiationMapsMaker::m_tgnSiA
TGraph * m_tgnSiA
Definition: RadiationMapsMaker.h:263
G4UA::RadiationMapsMaker::Config::nBinslogEn
int nBinslogEn
Definition: RadiationMapsMaker.h:69
G4UA
for nSW
Definition: CalibrationDefaultProcessing.h:19
G4UA::RadiationMapsMaker::RadiationMapsMaker
RadiationMapsMaker(const Config &config)
Definition: RadiationMapsMaker.cxx:62
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_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::Config::logEMinn
double logEMinn
Definition: RadiationMapsMaker.h:70
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::m_tgnSiC
TGraph * m_tgnSiC
Definition: RadiationMapsMaker.h:265
G4UA::RadiationMapsMaker::Config::activationFileName
std::string activationFileName
Definition: RadiationMapsMaker.h:34
G4UA::RadiationMapsMaker::Config::rMinZoom
double rMinZoom
Definition: RadiationMapsMaker.h:41
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:71
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
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
G4UA::RadiationMapsMaker::Config::logTMax
double logTMax
Definition: RadiationMapsMaker.h:81
cm3
#define cm3
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
M_PI
#define M_PI
Definition: ActiveFraction.h:11
G4UA::RadiationMapsMaker::m_tgHpm
TGraph * m_tgHpm
Definition: RadiationMapsMaker.h:276
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
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_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_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::Config::nBinstheta
int nBinstheta
Definition: RadiationMapsMaker.h:63
G4UA::RadiationMapsMaker::BeginOfRunAction
virtual void BeginOfRunAction(const G4Run *) override final
Definition: RadiationMapsMaker.cxx:255
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_eion
std::vector< double > m_rz_eion
vector of ionizing energy density seen by thread in zoomed area
Definition: RadiationMapsMaker.h:95
G4UA::RadiationMapsMaker::Config::nBinsz3d
int nBinsz3d
Definition: RadiationMapsMaker.h:54
G4UA::RadiationMapsMaker::m_tgnSiB
TGraph * m_tgnSiB
Definition: RadiationMapsMaker.h:264
G4UA::RadiationMapsMaker::Config::zMinZoom
double zMinZoom
Definition: RadiationMapsMaker.h:47
G4UA::RadiationMapsMaker::Report::m_3d_norm
std::vector< double > m_3d_norm
vector to normalize the volume fraction in 3d
Definition: RadiationMapsMaker.h:150
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
G4UA::RadiationMapsMaker::m_tgpSiB
TGraph * m_tgpSiB
Definition: RadiationMapsMaker.h:262
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
PlotCalibFromCool.ie
ie
Definition: PlotCalibFromCool.py:420
G4UA::RadiationMapsMaker::Config::materials
std::vector< std::string > materials
bin sizes and ranges match the requirements for the Radiation Estimate Web tool for the default value...
Definition: RadiationMapsMaker.h:32
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:200
G4UA::RadiationMapsMaker::Config::thetaMin
double thetaMin
Definition: RadiationMapsMaker.h:65
G4UA::RadiationMapsMaker::Config::rMaxFull
double rMaxFull
Definition: RadiationMapsMaker.h:45
G4UA::RadiationMapsMaker::m_tgHg
TGraph * m_tgHg
Definition: RadiationMapsMaker.h:270
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_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::Config::zMaxFull
double zMaxFull
Definition: RadiationMapsMaker.h:51
G4UA::RadiationMapsMaker::Config
Definition: RadiationMapsMaker.h:27
lumiFormat.i
int i
Definition: lumiFormat.py:92
G4UA::RadiationMapsMaker::UserSteppingAction
virtual void UserSteppingAction(const G4Step *) override final
Definition: RadiationMapsMaker.cxx:460
G4UA::RadiationMapsMaker::m_tgpSiA
TGraph * m_tgpSiA
Definition: RadiationMapsMaker.h:261
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
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
G4UA::RadiationMapsMaker::m_tgHem
TGraph * m_tgHem
Definition: RadiationMapsMaker.h:272
G4UA::RadiationMapsMaker::Config::nBinsz
int nBinsz
Definition: RadiationMapsMaker.h:39
G4UA::RadiationMapsMaker::Config::zMinFull
double zMinFull
Definition: RadiationMapsMaker.h:48
G4UA::RadiationMapsMaker::Config::phiMaxZoom
double phiMaxZoom
Definition: RadiationMapsMaker.h:58
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::Config::nBinsr
int nBinsr
Definition: RadiationMapsMaker.h:38
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::~RadiationMapsMaker
~RadiationMapsMaker()
Definition: RadiationMapsMaker.cxx:70
G4UA::RadiationMapsMaker::Report::merge
void merge(const Report &maps)
Definition: RadiationMapsMaker.cxx:155
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::m_tgHn
TGraph * m_tgHn
Definition: RadiationMapsMaker.h:269
G4UA::RadiationMapsMaker::m_tgHpp
TGraph * m_tgHpp
Definition: RadiationMapsMaker.h:277
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
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
G4UA::RadiationMapsMaker::Config::elemZMax
int elemZMax
Definition: RadiationMapsMaker.h:85
G4UA::RadiationMapsMaker::Report
Simple struct for holding the radiation maps.
Definition: RadiationMapsMaker.h:91
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::m_tgpiSi
TGraph * m_tgpiSi
Definition: RadiationMapsMaker.h:266
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_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_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::m_tgHp
TGraph * m_tgHp
Definition: RadiationMapsMaker.h:271
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
RadiationMapsMaker.h
PathResolver.h
G4UA::RadiationMapsMaker::Config::rMinFull
double rMinFull
Definition: RadiationMapsMaker.h:42
G4UA::RadiationMapsMaker::Config::logEMaxo
double logEMaxo
Definition: RadiationMapsMaker.h:76
G4UA::RadiationMapsMaker::Config::rMaxZoom
double rMaxZoom
Definition: RadiationMapsMaker.h:44
G4UA::RadiationMapsMaker::Config::logEMaxn
double logEMaxn
Definition: RadiationMapsMaker.h:71
G4UA::RadiationMapsMaker::Config::thetaMax
double thetaMax
Definition: RadiationMapsMaker.h:66
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
charge
double charge(const T &p)
Definition: AtlasPID.h:494
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
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_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_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::m_config
Config m_config
Definition: RadiationMapsMaker.h:257
G4UA::RadiationMapsMaker::m_tgHep
TGraph * m_tgHep
Definition: RadiationMapsMaker.h:273
ir
int ir
counter of the current depth
Definition: fastadd.cxx:49
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
G4UA::RadiationMapsMaker::Config::logTMin
double logTMin
Definition: RadiationMapsMaker.h:80
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
G4UA::RadiationMapsMaker::Config::nBinslogEo
int nBinslogEo
Definition: RadiationMapsMaker.h:74
G4UA::RadiationMapsMaker::m_maps
Report m_maps
Definition: RadiationMapsMaker.h:259
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_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_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_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_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_neut
std::vector< double > m_rz_neut
vector of >100 keV hadron flux seen by thread in zoomed area
Definition: RadiationMapsMaker.h:101
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
G4UA::RadiationMapsMaker::Config::posYOnly
bool posYOnly
Definition: RadiationMapsMaker.h:36
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::Config::nBinsr3d
int nBinsr3d
Definition: RadiationMapsMaker.h:53
G4UA::RadiationMapsMaker::Config::phiMinZoom
double phiMinZoom
Definition: RadiationMapsMaker.h:57
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
rr
const boost::regex rr(r_r)
G4UA::RadiationMapsMaker::Config::nBinsphi3d
int nBinsphi3d
Definition: RadiationMapsMaker.h:55
G4UA::RadiationMapsMaker::Config::zMaxZoom
double zMaxZoom
Definition: RadiationMapsMaker.h:50
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::m_tgHmm
TGraph * m_tgHmm
Definition: RadiationMapsMaker.h:274
G4UA::RadiationMapsMaker::Config::nBinsdphi
int nBinsdphi
Definition: RadiationMapsMaker.h:62
G4UA::RadiationMapsMaker::Config::elemZMin
int elemZMin
Definition: RadiationMapsMaker.h:84
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_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_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_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_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
fitman.rho
rho
Definition: fitman.py:532
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