ATLAS Offline Software
Loading...
Searching...
No Matches
Classes | Public Member Functions | Private Attributes | List of all members
G4UA::RadiationMapsMaker Class Reference

#include <RadiationMapsMaker.h>

Inheritance diagram for G4UA::RadiationMapsMaker:
Collaboration diagram for G4UA::RadiationMapsMaker:

Classes

struct  Config
struct  Report
 Simple struct for holding the radiation maps. More...

Public Member Functions

 RadiationMapsMaker (const Config &config)
 ~RadiationMapsMaker ()
virtual void BeginOfRunAction (const G4Run *) override final
virtual void UserSteppingAction (const G4Step *) override final
const Report & getReport () const
 Retrieve my maps.

Private Attributes

Config m_config
Report m_maps
TGraph * m_tgpSiA = 0
TGraph * m_tgpSiB = 0
TGraph * m_tgnSiA = 0
TGraph * m_tgnSiB = 0
TGraph * m_tgnSiC = 0
TGraph * m_tgpiSi = 0
TGraph * m_tgeSi = 0
TGraph * m_tgHn = 0
TGraph * m_tgHg = 0
TGraph * m_tgHp = 0
TGraph * m_tgHem = 0
TGraph * m_tgHep = 0
TGraph * m_tgHmm = 0
TGraph * m_tgHmp = 0
TGraph * m_tgHpm = 0
TGraph * m_tgHpp = 0
TGraph * m_tgHa = 0
std::ofstream m_activationFile

Detailed Description

Definition at line 19 of file RadiationMapsMaker.h.

Constructor & Destructor Documentation

◆ RadiationMapsMaker()

G4UA::RadiationMapsMaker::RadiationMapsMaker ( const Config & config)

Definition at line 62 of file RadiationMapsMaker.cxx.

63 : m_config(config)
64 {
65 }

◆ ~RadiationMapsMaker()

G4UA::RadiationMapsMaker::~RadiationMapsMaker ( )

Definition at line 70 of file RadiationMapsMaker.cxx.

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 }

Member Function Documentation

◆ BeginOfRunAction()

void G4UA::RadiationMapsMaker::BeginOfRunAction ( const G4Run * )
finaloverridevirtual

data files for NIEL damage factors in Si

The data resides in the share directory of the package and should not be changed. To prevent modification of the damage factors by accident the file names are not configurable

data files for effective dose in Humans from ICRP 116

Definition at line 255 of file RadiationMapsMaker.cxx.

255 {
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 }
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path)
Definition PathResolver.cxx:220

◆ getReport()

const Report & G4UA::RadiationMapsMaker::getReport ( ) const
inline

Retrieve my maps.

Definition at line 252 of file RadiationMapsMaker.h.

253 { return m_maps; }

◆ UserSteppingAction()

void G4UA::RadiationMapsMaker::UserSteppingAction ( const G4Step * aStep)
finaloverridevirtual

Definition at line 460 of file RadiationMapsMaker.cxx.

460 {
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 }
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
ir
int ir
counter of the current depth
Definition fastadd.cxx:49
PrintTrkAnaSummary.l
l
Printing final latex table to .tex output file.
Definition PrintTrkAnaSummary.py:370
dqt_zlumi_pandas.weight
int weight
Definition dqt_zlumi_pandas.py:190
makeDTCalibBlob_pickPhase.pd
pd
Definition makeDTCalibBlob_pickPhase.py:342

Member Data Documentation

◆ m_activationFile

std::ofstream G4UA::RadiationMapsMaker::m_activationFile
private

Definition at line 280 of file RadiationMapsMaker.h.

◆ m_config

Config G4UA::RadiationMapsMaker::m_config
private

Definition at line 257 of file RadiationMapsMaker.h.

◆ m_maps

Report G4UA::RadiationMapsMaker::m_maps
private

Definition at line 259 of file RadiationMapsMaker.h.

◆ m_tgeSi

TGraph* G4UA::RadiationMapsMaker::m_tgeSi = 0
private

Definition at line 267 of file RadiationMapsMaker.h.

◆ m_tgHa

TGraph* G4UA::RadiationMapsMaker::m_tgHa = 0
private

Definition at line 278 of file RadiationMapsMaker.h.

◆ m_tgHem

TGraph* G4UA::RadiationMapsMaker::m_tgHem = 0
private

Definition at line 272 of file RadiationMapsMaker.h.

◆ m_tgHep

TGraph* G4UA::RadiationMapsMaker::m_tgHep = 0
private

Definition at line 273 of file RadiationMapsMaker.h.

◆ m_tgHg

TGraph* G4UA::RadiationMapsMaker::m_tgHg = 0
private

Definition at line 270 of file RadiationMapsMaker.h.

◆ m_tgHmm

TGraph* G4UA::RadiationMapsMaker::m_tgHmm = 0
private

Definition at line 274 of file RadiationMapsMaker.h.

◆ m_tgHmp

TGraph* G4UA::RadiationMapsMaker::m_tgHmp = 0
private

Definition at line 275 of file RadiationMapsMaker.h.

◆ m_tgHn

TGraph* G4UA::RadiationMapsMaker::m_tgHn = 0
private

Definition at line 269 of file RadiationMapsMaker.h.

◆ m_tgHp

TGraph* G4UA::RadiationMapsMaker::m_tgHp = 0
private

Definition at line 271 of file RadiationMapsMaker.h.

◆ m_tgHpm

TGraph* G4UA::RadiationMapsMaker::m_tgHpm = 0
private

Definition at line 276 of file RadiationMapsMaker.h.

◆ m_tgHpp

TGraph* G4UA::RadiationMapsMaker::m_tgHpp = 0
private

Definition at line 277 of file RadiationMapsMaker.h.

◆ m_tgnSiA

TGraph* G4UA::RadiationMapsMaker::m_tgnSiA = 0
private

Definition at line 263 of file RadiationMapsMaker.h.

◆ m_tgnSiB

TGraph* G4UA::RadiationMapsMaker::m_tgnSiB = 0
private

Definition at line 264 of file RadiationMapsMaker.h.

◆ m_tgnSiC

TGraph* G4UA::RadiationMapsMaker::m_tgnSiC = 0
private

Definition at line 265 of file RadiationMapsMaker.h.

◆ m_tgpiSi

TGraph* G4UA::RadiationMapsMaker::m_tgpiSi = 0
private

Definition at line 266 of file RadiationMapsMaker.h.

◆ m_tgpSiA

TGraph* G4UA::RadiationMapsMaker::m_tgpSiA = 0
private

Definition at line 261 of file RadiationMapsMaker.h.

◆ m_tgpSiB

TGraph* G4UA::RadiationMapsMaker::m_tgpSiB = 0
private

Definition at line 262 of file RadiationMapsMaker.h.

The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.14.0