From 621bea712a09e88956e8ca4561c459b0316a57e2 Mon Sep 17 00:00:00 2001
From: Sven Menke <sven.menke@cern.ch>
Date: Mon, 29 Jun 2020 12:09:37 +0000
Subject: [PATCH] RadiationMapsMaker: Direction info for spectra and memory
footprint improvement in finalize
---
.../share/jobOptions.RadiationMapsMaker.py | 4 +
.../G4UserActions/src/RadiationMapsMaker.cxx | 86 +-
.../G4UserActions/src/RadiationMapsMaker.h | 26 +-
.../src/RadiationMapsMakerTool.cxx | 984 ++++++++----------
.../src/RadiationMapsMakerTool.h | 2 +-
5 files changed, 566 insertions(+), 536 deletions(-)
diff --git a/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py b/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
index c72dd9e6702..d8dc202a3bf 100644
--- a/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
+++ b/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
@@ -15,6 +15,8 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
# radmaptool.NBinsLogEn = 90
# radmaptool.NBinsLogEo = 45
# radmaptool.NBinsLogTimeCut = 20
+# radmaptool.NBinsDPhi = 18
+# radmaptool.NBinsTheta = 9
# radmaptool.NBinsR3D = 30
# radmaptool.NBinsZ3D = 60
# radmaptool.NBinsPhi3D = 32
@@ -34,6 +36,8 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
# radmaptool.LogEMaxo = 7.0 # in log10(E/MeV)
# radmaptool.LogTMin = -9.0 # in log10(t_cut/s)
# radmaptool.LogTMax = 11.0 # in log10(t_cut/s)
+# radmaptool.ThetaMin = 0.0 # in degrees
+# radmaptool.ThetaMin = 90.0 # in degrees
# radmaptool.ElemZMin = 1 # Atomic number
# radmaptool.ElemZMax = 92 # Atomic number
#
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
index 8bf2474349b..b309896e414 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
@@ -131,6 +131,10 @@ namespace G4UA{
m_rz_neut_spec [i] += maps.m_rz_neut_spec [i];
m_full_rz_neut_spec[i] += maps.m_full_rz_neut_spec[i];
}
+ // x theta bins
+ for(unsigned int i=0;i<maps.m_theta_full_rz_neut_spec.size();i++) {
+ m_theta_full_rz_neut_spec[i] += maps.m_theta_full_rz_neut_spec[i];
+ }
// all other particle's spectra have same size
for(unsigned int i=0;i<maps.m_rz_gamm_spec.size();i++) {
m_rz_gamm_spec [i] += maps.m_rz_gamm_spec [i];
@@ -146,6 +150,16 @@ namespace G4UA{
m_rz_rest_spec [i] += maps.m_rz_rest_spec [i];
m_full_rz_rest_spec[i] += maps.m_full_rz_rest_spec[i];
}
+ // x theta bins
+ for(unsigned int i=0;i<maps.m_theta_full_rz_gamm_spec.size();i++) {
+ m_theta_full_rz_gamm_spec[i] += maps.m_theta_full_rz_gamm_spec[i];
+ m_theta_full_rz_elec_spec[i] += maps.m_theta_full_rz_elec_spec[i];
+ m_theta_full_rz_muon_spec[i] += maps.m_theta_full_rz_muon_spec[i];
+ m_theta_full_rz_pion_spec[i] += maps.m_theta_full_rz_pion_spec[i];
+ m_theta_full_rz_prot_spec[i] += maps.m_theta_full_rz_prot_spec[i];
+ m_theta_full_rz_rchgd_spec[i] += maps.m_theta_full_rz_rchgd_spec[i];
+ m_theta_full_rz_rneut_spec[i] += maps.m_theta_full_rz_rneut_spec[i];
+ }
}
void RadiationMapsMaker::BeginOfRunAction(const G4Run*){
@@ -193,6 +207,15 @@ namespace G4UA{
m_maps.m_full_rz_prot_spec.resize(0);
m_maps.m_rz_rest_spec .resize(0);
m_maps.m_full_rz_rest_spec.resize(0);
+
+ m_maps.m_theta_full_rz_neut_spec .resize(0);
+ m_maps.m_theta_full_rz_gamm_spec .resize(0);
+ m_maps.m_theta_full_rz_elec_spec .resize(0);
+ m_maps.m_theta_full_rz_muon_spec .resize(0);
+ m_maps.m_theta_full_rz_pion_spec .resize(0);
+ m_maps.m_theta_full_rz_prot_spec .resize(0);
+ m_maps.m_theta_full_rz_rchgd_spec.resize(0);
+ m_maps.m_theta_full_rz_rneut_spec.resize(0);
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
@@ -248,6 +271,15 @@ namespace G4UA{
m_maps.m_rz_rest_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
m_maps.m_full_rz_rest_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// 2d zoom
@@ -407,7 +439,7 @@ namespace G4UA{
double y1 = aStep->GetPostStepPoint()->GetPosition().y()*0.1;
double z0 = aStep->GetPreStepPoint()->GetPosition().z()*0.1;
double z1 = aStep->GetPostStepPoint()->GetPosition().z()*0.1;
-
+
double l = sqrt(pow(x1-x0,2)+pow(y1-y0,2)+pow(z1-z0,2));
// make 1 mm steps but at least 1 step
double dl0 = 0.1;
@@ -419,6 +451,14 @@ namespace G4UA{
double weight = 0; // weight for NIEL
double eKin = aStep->GetTrack()->GetKineticEnergy();
+ double theta = aStep->GetTrack()->GetMomentumDirection().theta()*180./M_PI;
+ // if theta range in configuration is 0 - 90 assume that 180-theta should
+ // be used for theta > 90; otherwise leave theta unchanged
+ theta = (theta > 90&&m_config.thetaMin==0&&m_config.thetaMax==90?180-theta:theta);
+ double dphi = (aStep->GetTrack()->GetMomentumDirection().phi()-p.phi())*180./M_PI;
+ while ( dphi > 360 ) dphi -= 360.;
+ while ( dphi < 0 ) dphi += 360.;
+
double logEKin = (eKin > 0?log10(eKin):(m_config.logEMinn<m_config.logEMino?m_config.logEMinn:m_config.logEMino)-1);
if ( pdgid == 1 || pdgid == 9 ) {
@@ -472,6 +512,8 @@ namespace G4UA{
int vBinFullSpeco = -1;
int vBinZoomTime = -1;
int vBinFullTime = -1;
+ int vBinThetaFullSpecn = -1;
+ int vBinThetaFullSpeco = -1;
// zoom 2d
if ( m_config.zMinZoom < abszorz &&
@@ -517,12 +559,26 @@ namespace G4UA{
(pdgid == 6 || pdgid == 7)) {
int ile = (logEKin-m_config.logEMinn)/(m_config.logEMaxn-m_config.logEMinn)*m_config.nBinslogEn;
vBinFullSpecn = m_config.nBinsr*m_config.nBinslogEn*iz+ir*m_config.nBinslogEn+ile;
+ if ( m_config.thetaMin < theta &&
+ m_config.thetaMax > theta ) {
+ int ith = (theta-m_config.thetaMin)/(m_config.thetaMax-m_config.thetaMin)*m_config.nBinstheta;
+ int idph = dphi/360.*m_config.nBinsdphi;
+ ith = ith*m_config.nBinsdphi+idph;
+ 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;
+ }
}
if ( m_config.logEMino < logEKin &&
m_config.logEMaxo > logEKin &&
pdgid != 6 && pdgid != 7) {
int ile = (logEKin-m_config.logEMino)/(m_config.logEMaxo-m_config.logEMino)*m_config.nBinslogEo;
vBinFullSpeco = m_config.nBinsr*m_config.nBinslogEo*iz+ir*m_config.nBinslogEo+ile;
+ if ( m_config.thetaMin < theta &&
+ m_config.thetaMax > theta ) {
+ int ith = (theta-m_config.thetaMin)/(m_config.thetaMax-m_config.thetaMin)*m_config.nBinstheta;
+ int idph = dphi/360.*m_config.nBinsdphi;
+ ith = ith*m_config.nBinsdphi+idph;
+ 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;
+ }
}
}
}
@@ -685,6 +741,9 @@ namespace G4UA{
if ( vBinFullSpecn >=0 ) {
m_maps.m_full_rz_neut_spec [vBinFullSpecn] += dl;
}
+ if ( vBinThetaFullSpecn >=0 ) {
+ m_maps.m_theta_full_rz_neut_spec [vBinThetaFullSpecn] += dl;
+ }
}
}
@@ -730,6 +789,31 @@ namespace G4UA{
m_maps.m_full_rz_rest_spec [vBinFullSpeco] += dl;
}
}
+ if ( vBinThetaFullSpeco >=0 ) {
+ if ( pdgid == 0 ) {
+ m_maps.m_theta_full_rz_gamm_spec [vBinThetaFullSpeco] += dl;
+ }
+ else if ( pdgid == 1 ) {
+ m_maps.m_theta_full_rz_prot_spec [vBinThetaFullSpeco] += dl;
+ }
+ else if ( pdgid == 2 ) {
+ m_maps.m_theta_full_rz_pion_spec [vBinThetaFullSpeco] += dl;
+ }
+ else if ( pdgid == 3 ) {
+ m_maps.m_theta_full_rz_muon_spec [vBinThetaFullSpeco] += dl;
+ }
+ else if ( pdgid == 4 || pdgid == 5 ) {
+ m_maps.m_theta_full_rz_elec_spec [vBinThetaFullSpeco] += dl;
+ }
+ else {
+ if ( absq > 0 ) {
+ m_maps.m_theta_full_rz_rchgd_spec [vBinThetaFullSpeco] += dl;
+ }
+ else {
+ m_maps.m_theta_full_rz_rneut_spec [vBinThetaFullSpeco] += dl;
+ }
+ }
+ }
}
if (!m_config.material.empty()) {
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.h b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.h
index 682ba45e496..b718b060a5a 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.h
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.h
@@ -54,6 +54,14 @@ namespace G4UA
double phiMinZoom = -180.; // degrees
double phiMaxZoom = 180.; // degrees
+ // theta x dphi bins are used in the theta-spectra
+
+ int nBinsdphi = 18; // 0 degrees <= dphi < 360 degrees
+ int nBinstheta = 9;
+
+ double thetaMin = 0.; // degrees
+ double thetaMax = 90.; // degrees
+
// neutron spectra
int nBinslogEn = 90;
double logEMinn = -11.; // min log10(E_kin/MeV)
@@ -146,13 +154,17 @@ namespace G4UA
std::vector<double> m_rz_neut_spec;
/// vector of neutron spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_neut_spec;
+ /// vector of neutron spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_neut_spec;
// gamma
/// vector of gamma spectra in log10(E/MeV) bins and the zoom 2d grid
std::vector<double> m_rz_gamm_spec;
- /// vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid
+ /// vector of gamma spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_gamm_spec;
+ /// vector of gamma spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_gamm_spec;
// e^+/-
@@ -160,6 +172,8 @@ namespace G4UA
std::vector<double> m_rz_elec_spec;
/// vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_elec_spec;
+ /// vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_elec_spec;
// mu^+/-
@@ -167,6 +181,8 @@ namespace G4UA
std::vector<double> m_rz_muon_spec;
/// vector of mu^+/- spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_muon_spec;
+ /// vector of mu^+/- spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_muon_spec;
// pi^+/-
@@ -174,6 +190,8 @@ namespace G4UA
std::vector<double> m_rz_pion_spec;
/// vector of pi^+/- spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_pion_spec;
+ /// vector of pi^+/- spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_pion_spec;
// proton
@@ -181,6 +199,8 @@ namespace G4UA
std::vector<double> m_rz_prot_spec;
/// vector of proton spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_prot_spec;
+ /// vector of proton spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_prot_spec;
// rest
@@ -188,6 +208,10 @@ namespace G4UA
std::vector<double> m_rz_rest_spec;
/// vector of e^+/- spectra in log10(E/MeV) bins and the full 2d grid
std::vector<double> m_full_rz_rest_spec;
+ /// vector of rest charged spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_rchgd_spec;
+ /// vector of rest neutral spectra in log10(E/MeV) bins and the full 2d grid x theta bins
+ std::vector<double> m_theta_full_rz_rneut_spec;
// time dependent maps
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
index 0940d3c2f84..b1cd1fd0755 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
@@ -18,6 +18,7 @@ namespace G4UA{
declareProperty("RadMapsFileName", m_radMapsFileName);
/// Name of the material to make radiation maps for (take all if empty)
declareProperty("Material" , m_config.material);
+ /// If true consider hits with y>0 only -- useful for shafts
declareProperty("PositiveYOnly" , m_config.posYOnly);
/// map granularities
/// number of bins in r and z for all 2D maps
@@ -27,6 +28,10 @@ namespace G4UA{
declareProperty("NBinsLogEn" , m_config.nBinslogEn);
/// number of bins in logE for energy spectra of other particles in 2D grids
declareProperty("NBinsLogEo" , m_config.nBinslogEo);
+ /// number of bins in dphi for dphi x theta dependent energy spectra
+ declareProperty("NBinsDPhi" , m_config.nBinsdphi);
+ /// number of bins in theta for dphi x theta dependent energy spectra
+ declareProperty("NBinsTheta" , m_config.nBinstheta);
/// number of bins in r, z and phi for all 3D maps
declareProperty("NBinsR3D" , m_config.nBinsr3d);
declareProperty("NBinsZ3D" , m_config.nBinsz3d);
@@ -53,6 +58,9 @@ namespace G4UA{
/// for logE of other particles in 2D spectra
declareProperty("LogEMino" , m_config.logEMino);
declareProperty("LogEMaxo" , m_config.logEMaxo);
+ /// for Theta in 2D spectra
+ declareProperty("ThetaMin" , m_config.thetaMin);
+ declareProperty("ThetaMax" , m_config.thetaMax);
/// for logT in time-dependent TID 2D maps
declareProperty("LogTMin" , m_config.logTMin);
declareProperty("LogTMax" , m_config.logTMax);
@@ -80,6 +88,9 @@ namespace G4UA{
m_config.logEMinn << " < log10(E/MeV) < " << m_config.logEMaxn << "\n" <<
"Other Spectra: " << m_config.nBinslogEo << " log10E-bins" << ", " <<
m_config.logEMino << " < log10(E/MeV) < " << m_config.logEMaxo << "\n" <<
+ "DPhi-Bins: " << m_config.nBinsdphi << " dphi-bins, 0 < dphi < 360 \n" <<
+ "Theta-Bins: " << m_config.nBinstheta << " theta-bins" << ", " <<
+ m_config.thetaMin << " < theta < " << m_config.thetaMax << "\n" <<
"3D Maps: " << m_config.nBinsz3d << (m_config.zMinFull<0?" z-bins, ":" |z|-bins, ") <<
m_config.nBinsr3d << " r-bins, " <<
m_config.nBinsphi3d << " phi-bins" << "\n" <<
@@ -101,46 +112,72 @@ namespace G4UA{
{
ATH_MSG_DEBUG( "Finalizing " << name() );
- // first make sure the vectors are empty
-
RadiationMapsMaker::Report maps;
- maps.m_rz_tid .resize(0);
- maps.m_rz_eion.resize(0);
- maps.m_rz_niel.resize(0);
- maps.m_rz_h20 .resize(0);
- maps.m_rz_neut.resize(0);
- maps.m_rz_chad.resize(0);
-
- maps.m_full_rz_tid .resize(0);
- maps.m_full_rz_eion.resize(0);
- maps.m_full_rz_niel.resize(0);
- maps.m_full_rz_h20 .resize(0);
- maps.m_full_rz_neut.resize(0);
- maps.m_full_rz_chad.resize(0);
-
- maps.m_3d_tid .resize(0);
- maps.m_3d_eion.resize(0);
- maps.m_3d_niel.resize(0);
- maps.m_3d_h20 .resize(0);
- maps.m_3d_neut.resize(0);
- maps.m_3d_chad.resize(0);
-
- maps.m_rz_neut_spec .resize(0);
- maps.m_full_rz_neut_spec.resize(0);
- maps.m_rz_gamm_spec .resize(0);
- maps.m_full_rz_gamm_spec.resize(0);
- maps.m_rz_elec_spec .resize(0);
- maps.m_full_rz_elec_spec.resize(0);
- maps.m_rz_muon_spec .resize(0);
- maps.m_full_rz_muon_spec.resize(0);
- maps.m_rz_pion_spec .resize(0);
- maps.m_full_rz_pion_spec.resize(0);
- maps.m_rz_prot_spec .resize(0);
- maps.m_full_rz_prot_spec.resize(0);
- maps.m_rz_rest_spec .resize(0);
- maps.m_full_rz_rest_spec.resize(0);
+ // vector of pointers to vectors of double to save space
+ std::vector<std::vector<double> *> pVal = {
+ &(maps.m_rz_tid),&(maps.m_rz_eion),&(maps.m_rz_niel),&(maps.m_rz_h20),&(maps.m_rz_neut),&(maps.m_rz_chad),
+ &(maps.m_full_rz_tid),&(maps.m_full_rz_eion),&(maps.m_full_rz_niel),&(maps.m_full_rz_h20),&(maps.m_full_rz_neut),&(maps.m_full_rz_chad)
+ };
+
+ std::vector<std::vector<double> *> pVal3d = {
+ &(maps.m_3d_tid),&(maps.m_3d_eion),&(maps.m_3d_niel),&(maps.m_3d_h20),&(maps.m_3d_neut),&(maps.m_3d_chad)
+ };
+
+ std::vector<std::vector<double> *> pValSpec = {
+ &(maps.m_rz_neut_spec),&(maps.m_rz_gamm_spec),&(maps.m_rz_elec_spec),&(maps.m_rz_muon_spec),&(maps.m_rz_pion_spec),&(maps.m_rz_prot_spec),&(maps.m_rz_rest_spec),
+ &(maps.m_full_rz_neut_spec),&(maps.m_full_rz_gamm_spec),&(maps.m_full_rz_elec_spec),&(maps.m_full_rz_muon_spec),&(maps.m_full_rz_pion_spec),&(maps.m_full_rz_prot_spec),&(maps.m_full_rz_rest_spec)
+ };
+
+ std::vector<int> nBinslogE = {
+ m_config.nBinslogEn,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,
+ m_config.nBinslogEn,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo
+ };
+
+ std::vector<std::vector<double> *> pValThetaSpec = {
+ &(maps.m_theta_full_rz_neut_spec),&(maps.m_theta_full_rz_gamm_spec),&(maps.m_theta_full_rz_elec_spec),&(maps.m_theta_full_rz_muon_spec),&(maps.m_theta_full_rz_pion_spec),&(maps.m_theta_full_rz_prot_spec),&(maps.m_theta_full_rz_rchgd_spec),&(maps.m_theta_full_rz_rneut_spec)
+ };
+
+ std::vector<int> nBinsThetalogE = {
+ m_config.nBinslogEn,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo,m_config.nBinslogEo
+ };
+
+ for(unsigned int hi=0;hi<pVal.size();hi++) {
+
+ // first make sure the vectors are empty
+
+ (pVal[hi])->resize(0);
+
+ // then resize to proper size and initialize with 0's
+ // all maps are needed for the merge - so have to do resize
+ // for all before merging any ...
+
+ (pVal[hi])->resize(m_config.nBinsz*m_config.nBinsr,0.0);
+ }
+
+ // same for 3d
+
+ for(unsigned int hi=0;hi<pVal3d.size();hi++) {
+ (pVal3d[hi])->resize(0);
+ (pVal3d[hi])->resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
+ }
+
+ // same for spectra
+ for(unsigned int hi=0;hi<pValSpec.size();hi++) {
+ (pValSpec[hi])->resize(0);
+ (pValSpec[hi])->resize(m_config.nBinsz*m_config.nBinsr*nBinslogE[hi],0.0);
+ }
+
+ // same for theta x dphi spectra
+
+ for(unsigned int hi=0;hi<pValThetaSpec.size();hi++) {
+ (pValThetaSpec[hi])->resize(0);
+ (pValThetaSpec[hi])->resize(m_config.nBinsdphi*m_config.nBinstheta*m_config.nBinsz*m_config.nBinsr*nBinsThetalogE[hi],0.0);
+ }
+
+ // same for misc individual maps
+
maps.m_rz_tid_time .resize(0);
maps.m_full_rz_tid_time .resize(0);
@@ -161,42 +198,8 @@ namespace G4UA{
}
// then resize to proper size and initialize with 0's
-
- maps.m_rz_tid .resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_rz_eion.resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_rz_niel.resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_rz_h20 .resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_rz_neut.resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_rz_chad.resize(m_config.nBinsz*m_config.nBinsr,0.0);
-
- maps.m_full_rz_tid .resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_full_rz_eion.resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_full_rz_niel.resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_full_rz_h20 .resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_full_rz_neut.resize(m_config.nBinsz*m_config.nBinsr,0.0);
- maps.m_full_rz_chad.resize(m_config.nBinsz*m_config.nBinsr,0.0);
-
- maps.m_3d_tid .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
- maps.m_3d_eion.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
- maps.m_3d_niel.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
- maps.m_3d_h20 .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
- maps.m_3d_neut.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
- maps.m_3d_chad.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
-
- maps.m_rz_neut_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
- maps.m_full_rz_neut_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
- maps.m_rz_gamm_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_full_rz_gamm_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_rz_elec_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_full_rz_elec_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_rz_muon_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_full_rz_muon_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_rz_pion_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_full_rz_pion_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_rz_prot_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_full_rz_prot_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_rz_rest_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
- maps.m_full_rz_rest_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ // all maps are needed for the merge - so have to do resize
+ // for all first ...
maps.m_rz_tid_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
maps.m_full_rz_tid_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
@@ -224,321 +227,378 @@ namespace G4UA{
TFile * f = new TFile(m_radMapsFileName.c_str(),"RECREATE");
- TH2D * h_rz_tid = new TH2D("rz_tid" ,"rz_tid" ,m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- TH2D * h_rz_eion = new TH2D("rz_eion","rz_eion",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- TH2D * h_rz_niel = new TH2D("rz_niel","rz_niel",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- TH2D * h_rz_h20 = new TH2D("rz_h20" ,"rz_h20" ,m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- TH2D * h_rz_neut = new TH2D("rz_neut","rz_neut",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- TH2D * h_rz_chad = new TH2D("rz_chad","rz_chad",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
+ // histograms can be booked, filled, written and deleted in a loop
+ // in order to save memory. Need vectors to loop over them ...
+
+ std::vector<const char *> h2dNames = {
+ "rz_tid","rz_eion","rz_niel","rz_h20","rz_neut","rz_chad",
+ "full_rz_tid","full_rz_eion","full_rz_niel","full_rz_h20","full_rz_neut","full_rz_chad"
+ };
+
+ std::vector<const char *> h2dZTitles = {
+ "TID [Gy]","E_{ion}/V [MeV/cm^{3}]","NIEL [n_{eq}/cm^{2}]","SEE [h_{>20 MeV}/cm^{2}]","FLUX [n_{>100 keV}/cm^{2}]","FLUX [h_{charged}/cm^{2}]",
+ "TID [Gy]","E_{ion}/V [MeV/cm^{3}]","NIEL [n_{eq}/cm^{2}]","SEE [h_{>20 MeV}/cm^{2}]","FLUX [n_{>100 keV}/cm^{2}]","FLUX [h_{charged}/cm^{2}]"
+ };
+
+ std::vector<double> zMin = {
+ m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,
+ m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull
+ };
+
+ std::vector<double> zMax = {
+ m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,
+ m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull
+ };
+
+ std::vector<double> rMin = {
+ m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,
+ m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull
+ };
+
+ std::vector<double> rMax = {
+ m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,
+ m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull
+ };
const char * xtitle = (m_config.zMinFull<0?"z [cm]":"|z| [cm]");
- h_rz_tid ->SetXTitle(xtitle);
- h_rz_eion ->SetXTitle(xtitle);
- h_rz_niel ->SetXTitle(xtitle);
- h_rz_h20 ->SetXTitle(xtitle);
- h_rz_neut ->SetXTitle(xtitle);
- h_rz_chad ->SetXTitle(xtitle);
-
- h_rz_tid ->SetYTitle("r [cm]");
- h_rz_eion ->SetYTitle("r [cm]");
- h_rz_niel ->SetYTitle("r [cm]");
- h_rz_h20 ->SetYTitle("r [cm]");
- h_rz_neut ->SetYTitle("r [cm]");
- h_rz_chad ->SetYTitle("r [cm]");
-
- h_rz_tid ->SetZTitle("TID [Gy]");
- h_rz_eion ->SetZTitle("E_{ion}/V [MeV/cm^{3}]");
- h_rz_niel ->SetZTitle("NIEL [n_{eq}/cm^{2}]");
- h_rz_h20 ->SetZTitle("SEE [h_{>20 MeV}/cm^{2}]");
- h_rz_neut ->SetZTitle("FLUX [n_{>100 keV}/cm^{2}]");
- h_rz_chad ->SetZTitle("FLUX [h_{charged}/cm^{2}]");
-
- TH2D *h_full_rz_tid = new TH2D("full_rz_tid" ,"full_rz_tid" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- TH2D *h_full_rz_eion = new TH2D("full_rz_eion","full_rz_eion",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- TH2D *h_full_rz_niel = new TH2D("full_rz_niel","full_rz_niel",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- TH2D *h_full_rz_h20 = new TH2D("full_rz_h20" ,"full_rz_h20" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- TH2D *h_full_rz_neut = new TH2D("full_rz_neut","full_rz_neut",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- TH2D *h_full_rz_chad = new TH2D("full_rz_chad","full_rz_chad",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
-
- h_full_rz_tid ->SetXTitle(xtitle);
- h_full_rz_eion ->SetXTitle(xtitle);
- h_full_rz_niel ->SetXTitle(xtitle);
- h_full_rz_h20 ->SetXTitle(xtitle);
- h_full_rz_neut ->SetXTitle(xtitle);
- h_full_rz_chad ->SetXTitle(xtitle);
-
- h_full_rz_tid ->SetYTitle("r [cm]");
- h_full_rz_eion ->SetYTitle("r [cm]");
- h_full_rz_niel ->SetYTitle("r [cm]");
- h_full_rz_h20 ->SetYTitle("r [cm]");
- h_full_rz_neut ->SetYTitle("r [cm]");
- h_full_rz_chad ->SetYTitle("r [cm]");
-
- h_full_rz_tid ->SetZTitle("TID [Gy]");
- h_full_rz_eion ->SetZTitle("E_{ion}/V [MeV/cm^{3}]");
- h_full_rz_niel ->SetZTitle("NIEL [n_{eq}/cm^{2}]");
- h_full_rz_h20 ->SetZTitle("SEE [h_{>20 MeV}/cm^{2}]");
- h_full_rz_neut ->SetZTitle("FLUX [n_{>100 keV}/cm^{2}]");
- h_full_rz_chad ->SetZTitle("FLUX [h_{charged}/cm^{2}]");
-
- TH3D * h_3d_tid = new TH3D("h3d_tid" ,"h3d_tid" ,m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
- TH3D * h_3d_eion = new TH3D("h3d_eion","h3d_eion",m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
- TH3D * h_3d_niel = new TH3D("h3d_niel","h3d_niel",m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
- TH3D * h_3d_h20 = new TH3D("h3d_h20" ,"h3d_h20" ,m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
- TH3D * h_3d_neut = new TH3D("h3d_neut","h3d_neut",m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
- TH3D * h_3d_chad = new TH3D("h3d_chad","h3d_chad",m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
-
- h_3d_tid ->SetXTitle(xtitle);
- h_3d_eion ->SetXTitle(xtitle);
- h_3d_niel ->SetXTitle(xtitle);
- h_3d_h20 ->SetXTitle(xtitle);
- h_3d_neut ->SetXTitle(xtitle);
- h_3d_chad ->SetXTitle(xtitle);
-
- h_3d_tid ->SetYTitle("r [cm]");
- h_3d_eion ->SetYTitle("r [cm]");
- h_3d_niel ->SetYTitle("r [cm]");
- h_3d_h20 ->SetYTitle("r [cm]");
- h_3d_neut ->SetYTitle("r [cm]");
- h_3d_chad ->SetYTitle("r [cm]");
-
- h_3d_tid ->SetZTitle("#phi [#circ]");
- h_3d_eion ->SetZTitle("#phi [#circ]");
- h_3d_niel ->SetZTitle("#phi [#circ]");
- h_3d_h20 ->SetZTitle("#phi [#circ]");
- h_3d_neut ->SetZTitle("#phi [#circ]");
- h_3d_chad ->SetZTitle("#phi [#circ]");
-
- h_3d_tid ->SetTitle("TID [Gy]");
- h_3d_eion ->SetTitle("E_{ion}/V [MeV/cm^{3}]");
- h_3d_niel ->SetTitle("NIEL [n_{eq}/cm^{2}]");
- h_3d_h20 ->SetTitle("SEE [h_{>20 MeV}/cm^{2}]");
- h_3d_neut ->SetTitle("FLUX [n_{>100 keV}/cm^{2}]");
- h_3d_chad ->SetTitle("FLUX [h_{charged}/cm^{2}]");
-
- // neutron spectra
- TH3D *h_full_rz_neut_spec = new TH3D("full_rz_neut_spec","full_rz_neut_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEn,m_config.logEMinn,m_config.logEMaxn);
- TH3D *h_rz_neut_spec = new TH3D("rz_neut_spec","rz_neut_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEn,m_config.logEMinn,m_config.logEMaxn);
-
- h_rz_neut_spec ->SetXTitle(xtitle);
- h_rz_neut_spec ->SetYTitle("r [cm]");
- h_rz_neut_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_neut_spec ->SetTitle("FLUX [n(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_neut_spec ->SetXTitle(xtitle);
- h_full_rz_neut_spec ->SetYTitle("r [cm]");
- h_full_rz_neut_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_neut_spec ->SetTitle("FLUX [n(log_{10}(E/MeV))/cm^{2}]");
-
- // gamma spectra
- TH3D *h_full_rz_gamm_spec = new TH3D("full_rz_gamm_spec","full_rz_gamm_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
- TH3D *h_rz_gamm_spec = new TH3D("rz_gamm_spec","rz_gamm_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
-
- h_rz_gamm_spec ->SetXTitle(xtitle);
- h_rz_gamm_spec ->SetYTitle("r [cm]");
- h_rz_gamm_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_gamm_spec ->SetTitle("FLUX [#gamma(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_gamm_spec ->SetXTitle(xtitle);
- h_full_rz_gamm_spec ->SetYTitle("r [cm]");
- h_full_rz_gamm_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_gamm_spec ->SetTitle("FLUX [#gamma(log_{10}(E/MeV))/cm^{2}]");
-
- // e^+/- spectra
- TH3D *h_full_rz_elec_spec = new TH3D("full_rz_elec_spec","full_rz_elec_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
- TH3D *h_rz_elec_spec = new TH3D("rz_elec_spec","rz_elec_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
-
- h_rz_elec_spec ->SetXTitle(xtitle);
- h_rz_elec_spec ->SetYTitle("r [cm]");
- h_rz_elec_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_elec_spec ->SetTitle("FLUX [e^{#pm}(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_elec_spec ->SetXTitle(xtitle);
- h_full_rz_elec_spec ->SetYTitle("r [cm]");
- h_full_rz_elec_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_elec_spec ->SetTitle("FLUX [e^{#pm}(log_{10}(E/MeV))/cm^{2}]");
-
- // mu^+/- spectra
- TH3D *h_full_rz_muon_spec = new TH3D("full_rz_muon_spec","full_rz_muon_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
- TH3D *h_rz_muon_spec = new TH3D("rz_muon_spec","rz_muon_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
-
- h_rz_muon_spec ->SetXTitle(xtitle);
- h_rz_muon_spec ->SetYTitle("r [cm]");
- h_rz_muon_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_muon_spec ->SetTitle("FLUX [#mu^{#pm}(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_muon_spec ->SetXTitle(xtitle);
- h_full_rz_muon_spec ->SetYTitle("r [cm]");
- h_full_rz_muon_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_muon_spec ->SetTitle("FLUX [#mu^{#pm}(log_{10}(E/MeV))/cm^{2}]");
-
- // pi^+/- spectra
- TH3D *h_full_rz_pion_spec = new TH3D("full_rz_pion_spec","full_rz_pion_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
- TH3D *h_rz_pion_spec = new TH3D("rz_pion_spec","rz_pion_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
-
- h_rz_pion_spec ->SetXTitle(xtitle);
- h_rz_pion_spec ->SetYTitle("r [cm]");
- h_rz_pion_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_pion_spec ->SetTitle("FLUX [#pi^{#pm}(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_pion_spec ->SetXTitle(xtitle);
- h_full_rz_pion_spec ->SetYTitle("r [cm]");
- h_full_rz_pion_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_pion_spec ->SetTitle("FLUX [#pi^{#pm}(log_{10}(E/MeV))/cm^{2}]");
-
- // proton spectra
- TH3D *h_full_rz_prot_spec = new TH3D("full_rz_prot_spec","full_rz_prot_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
- TH3D *h_rz_prot_spec = new TH3D("rz_prot_spec","rz_prot_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
-
- h_rz_prot_spec ->SetXTitle(xtitle);
- h_rz_prot_spec ->SetYTitle("r [cm]");
- h_rz_prot_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_prot_spec ->SetTitle("FLUX [p(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_prot_spec ->SetXTitle(xtitle);
- h_full_rz_prot_spec ->SetYTitle("r [cm]");
- h_full_rz_prot_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_prot_spec ->SetTitle("FLUX [p(log_{10}(E/MeV))/cm^{2}]");
-
- // rest spectra
- TH3D *h_full_rz_rest_spec = new TH3D("full_rz_rest_spec","full_rz_rest_spec",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
- TH3D *h_rz_rest_spec = new TH3D("rz_rest_spec","rz_rest_spec",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogEo,m_config.logEMino,m_config.logEMaxo);
-
- h_rz_rest_spec ->SetXTitle(xtitle);
- h_rz_rest_spec ->SetYTitle("r [cm]");
- h_rz_rest_spec ->SetZTitle("log_{10}(E/MeV)");
- h_rz_rest_spec ->SetTitle("FLUX [rest(log_{10}(E/MeV))/cm^{2}]");
-
- h_full_rz_rest_spec ->SetXTitle(xtitle);
- h_full_rz_rest_spec ->SetYTitle("r [cm]");
- h_full_rz_rest_spec ->SetZTitle("log_{10}(E/MeV)");
- h_full_rz_rest_spec ->SetTitle("FLUX [rest(log_{10}(E/MeV))/cm^{2}]");
-
- // time dependent TID maps
+ TH2D * h_rz_vol = 0;
+ TH2D * h_rz_norm = 0;
+ TH2D * h_full_rz_vol = 0;
+ TH2D * h_full_rz_norm = 0;
+
+ for(unsigned int hi=0;hi<h2dNames.size();hi++) {
+ TH2D *h2d = new TH2D(h2dNames[hi],h2dNames[hi],m_config.nBinsz,zMin[hi],zMax[hi],m_config.nBinsr,rMin[hi],rMax[hi]);
+ h2d->SetXTitle(xtitle);
+ h2d->SetYTitle("r [cm]");
+ h2d->SetZTitle(h2dZTitles[hi]);
+ if (hi==0 && !m_config.material.empty()) {
+ // need volume fraction only if particular material is selected
+ //
+ // the maps for TID, NIEL and SEE need to be divided by the ratio of (vol/norm) in order to get
+ // the proper estimate per volume bin for the selected material.
+ // This is *not* done in the tool directly and left to the user after having summed the histograms
+ // from many individual jobs.
+ //
+ h_rz_vol = new TH2D("rz_vol" ,"rz_vol" ,m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
+ h_rz_norm = new TH2D("rz_norm","rz_norm",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
+ h_full_rz_vol = new TH2D("full_rz_vol" ,"full_rz_vol" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
+ h_full_rz_norm = new TH2D("full_rz_norm","full_rz_norm",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
+
+ h_rz_vol ->SetXTitle(xtitle);
+ h_rz_norm ->SetXTitle(xtitle);
+ h_rz_vol ->SetYTitle("r [cm]");
+ h_rz_norm ->SetYTitle("r [cm]");
+ std::string hname("Volume fraction of ");
+ hname += m_config.material;
+ h_rz_vol ->SetZTitle(hname.data());
+ h_rz_norm ->SetZTitle("Volume norm");
+
+ h_full_rz_vol ->SetXTitle(xtitle);
+ h_full_rz_norm ->SetXTitle(xtitle);
+ h_full_rz_vol ->SetYTitle("r [cm]");
+ h_full_rz_norm ->SetYTitle("r [cm]");
+ h_full_rz_vol ->SetZTitle(hname.data());
+ h_full_rz_norm ->SetZTitle("Volume norm");
+ }
+
+ // normalize to volume element per bin
+ for(int i=0;i<h2d->GetNbinsX();i++) {
+ for(int j=0;j<h2d->GetNbinsY();j++) {
+ int iBin = h2d->GetBin(i+1,j+1);
+ int vBin = m_config.nBinsr*i+j;
+ double r0=h2d->GetYaxis()->GetBinLowEdge(j+1);
+ double r1=h2d->GetYaxis()->GetBinUpEdge(j+1);
+ double z0=h2d->GetXaxis()->GetBinLowEdge(i+1);
+ double z1=h2d->GetXaxis()->GetBinUpEdge(i+1);
+ double vol=(z1-z0)*M_PI*(r1*r1-r0*r0);
+ // if |z| instead of z double the volume
+ if ( m_config.zMinFull >= 0 ) vol *= 2;
+ // if positive y hemisphere is used only -- half the volume
+ if ( m_config.posYOnly ) vol *= 0.5;
+ double val = (*(pVal[hi]))[vBin];
+ h2d->SetBinContent(iBin,val/vol);
+ if (hi ==0 && !m_config.material.empty()) {
+ // need volume fraction only if particular material is selected
+ // VOL
+ val =maps.m_rz_vol[vBin];
+ h_rz_vol->SetBinContent(iBin,val/vol);
+ // NORM
+ val =maps.m_rz_norm[vBin];
+ h_rz_norm->SetBinContent(iBin,val/vol);
+ }
+ if (hi ==h2dNames.size()/2 && !m_config.material.empty()) {
+ // need volume fraction only if particular material is selected
+ // VOL
+ val =maps.m_full_rz_vol[vBin];
+ h_full_rz_vol->SetBinContent(iBin,val/vol);
+ // NORM
+ val =maps.m_full_rz_norm[vBin];
+ h_full_rz_norm->SetBinContent(iBin,val/vol);
+ }
+ }
+ }
+ h2d->Write();
+ h2d->Delete();
+ if (hi ==h2dNames.size()/2 && !m_config.material.empty()) {
+ h_rz_vol->Write();
+ h_rz_vol->Delete();
+ h_rz_norm->Write();
+ h_rz_norm->Delete();
+ h_full_rz_vol->Write();
+ h_full_rz_vol->Delete();
+ h_full_rz_norm->Write();
+ h_full_rz_norm->Delete();
+ }
+ }
+
+ std::vector<const char *> h3dNames = {
+ "h3d_tid","h3d_eion","h3d_niel","h3d_h20","h3d_neut","h3d_chad"
+ };
+
+ std::vector<const char *> h3dTitles = {
+ "TID [Gy]","E_{ion}/V [MeV/cm^{3}]","NIEL [n_{eq}/cm^{2}]","SEE [h_{>20 MeV}/cm^{2}]","FLUX [n_{>100 keV}/cm^{2}]","FLUX [h_{charged}/cm^{2}]"
+ };
+
+ TH3D * h_3d_vol = 0;
+ TH3D * h_3d_norm = 0;
+
+ for(unsigned int hi=0;hi<h3dNames.size();hi++) {
+ TH3D *h3d = new TH3D(h3dNames[hi],h3dNames[hi],m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
+ h3d->SetXTitle(xtitle);
+ h3d->SetYTitle("r [cm]");
+ h3d->SetZTitle("#phi [#circ]");
+ h3d->SetTitle(h3dTitles[hi]);
+ if (hi == 0 && !m_config.material.empty()) {
+ h_3d_vol = new TH3D("h3d_vol" ,"h3d_vol" ,m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
+ h_3d_norm = new TH3D("h3d_norm","h3d_norm",m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
+ h_3d_vol ->SetXTitle(xtitle);
+ h_3d_norm ->SetXTitle(xtitle);
+ h_3d_vol ->SetYTitle("r [cm]");
+ h_3d_norm ->SetYTitle("r [cm]");
+ h_3d_vol ->SetZTitle("#phi [#circ]");
+ h_3d_norm ->SetZTitle("#phi [#circ]");
+ std::string hname("Volume fraction of ");
+ hname += m_config.material;
+ h_3d_vol ->SetTitle(hname.data());
+ h_3d_norm ->SetTitle("Volume norm");
+ }
+ // normalize to volume element per bin
+ for(int i=0;i<h3d->GetNbinsX();i++) {
+ for(int j=0;j<h3d->GetNbinsY();j++) {
+ for(int k=0;k<h3d->GetNbinsZ();k++) {
+ int vBin = m_config.nBinsr3d*m_config.nBinsphi3d*i+m_config.nBinsphi3d*j+k;
+ int iBin = h3d->GetBin(i+1,j+1,k+1);
+ double phi0=h3d->GetZaxis()->GetBinLowEdge(k+1);
+ double phi1=h3d->GetZaxis()->GetBinUpEdge(k+1);
+ double r0=h3d->GetYaxis()->GetBinLowEdge(j+1);
+ double r1=h3d->GetYaxis()->GetBinUpEdge(j+1);
+ double z0=h3d->GetXaxis()->GetBinLowEdge(i+1);
+ double z1=h3d->GetXaxis()->GetBinUpEdge(i+1);
+ double vol=(z1-z0)*M_PI*(r1*r1-r0*r0)*(phi1-phi0)/360.;
+ // if |z| instead of z double the volume
+ if ( m_config.zMinFull >= 0 ) vol *= 2;
+ // assume that phi-range corresponds to full 360 degrees in case
+ // lower phi boundary is 0 - i.e. all phi-segments mapped to first
+ if ( m_config.phiMinZoom == 0 ) {
+ vol *= 360./m_config.phiMaxZoom;
+ // if positive y hemisphere is used only -- half the volume
+ if ( m_config.posYOnly )
+ vol *= 0.5;
+ }
+ double val = (*(pVal3d[hi]))[vBin];
+ h3d->SetBinContent(iBin,val/vol);
+ if (hi == 0 && !m_config.material.empty()) {
+ // VOL
+ val =maps.m_3d_vol[vBin];
+ h_3d_vol->SetBinContent(iBin,val/vol);
+ // NORM
+ val =maps.m_3d_norm[vBin];
+ h_3d_norm->SetBinContent(iBin,val/vol);
+ }
+ }
+ }
+ }
+ h3d->Write();
+ h3d->Delete();
+ if (hi == 0 && !m_config.material.empty()) {
+ h_3d_vol->Write();
+ h_3d_vol->Delete();
+ h_3d_norm->Write();
+ h_3d_norm->Delete();
+ }
+ }
+
+ // spectra
+
+ std::vector<const char *> hSpecNames = {
+ "rz_neut_spec","rz_gamm_spec","rz_elec_spec","rz_muon_spec","rz_pion_spec","rz_prot_spec","rz_rest_spec",
+ "full_rz_neut_spec","full_rz_gamm_spec","full_rz_elec_spec","full_rz_muon_spec","full_rz_pion_spec","full_rz_prot_spec","full_rz_rest_spec"
+ };
+
+ std::vector<const char *> hSpecTitles = {
+ "FLUX [n(log_{10}(E/MeV))/cm^{2}]","FLUX [#gamma(log_{10}(E/MeV))/cm^{2}]","FLUX [e^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [#mu^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [#pi^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [p(log_{10}(E/MeV))/cm^{2}]","FLUX [rest(log_{10}(E/MeV))/cm^{2}]",
+ "FLUX [n(log_{10}(E/MeV))/cm^{2}]","FLUX [#gamma(log_{10}(E/MeV))/cm^{2}]","FLUX [e^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [#mu^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [#pi^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [p(log_{10}(E/MeV))/cm^{2}]","FLUX [rest(log_{10}(E/MeV))/cm^{2}]"
+ };
+
+ std::vector<double> zMinSpec = {
+ m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,m_config.zMinZoom,
+ m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull,m_config.zMinFull
+ };
+
+ std::vector<double> zMaxSpec = {
+ m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,m_config.zMaxZoom,
+ m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull,m_config.zMaxFull
+ };
+
+ std::vector<double> rMinSpec = {
+ m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,m_config.rMinZoom,
+ m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull,m_config.rMinFull
+ };
+
+ std::vector<double> rMaxSpec = {
+ m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,m_config.rMaxZoom,
+ m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull,m_config.rMaxFull
+ };
+
+ std::vector<double> logEMin = {
+ m_config.logEMinn,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,
+ m_config.logEMinn,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino
+ };
+
+ std::vector<double> logEMax = {
+ m_config.logEMaxn,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,
+ m_config.logEMaxn,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo
+ };
+
+ for(unsigned int hi=0;hi<hSpecNames.size();hi++) {
+ TH3D *hSpec = new TH3D(hSpecNames[hi],hSpecNames[hi],m_config.nBinsz,zMinSpec[hi],zMaxSpec[hi],m_config.nBinsr,rMinSpec[hi],rMaxSpec[hi],nBinslogE[hi],logEMin[hi],logEMax[hi]);
+ hSpec->SetXTitle(xtitle);
+ hSpec->SetYTitle("r [cm]");
+ hSpec->SetZTitle("log_{10}(E/MeV)");
+ hSpec->SetTitle(hSpecTitles[hi]);
+ // normalize to volume element per bin
+ for(int i=0;i<hSpec->GetNbinsX();i++) {
+ for(int j=0;j<hSpec->GetNbinsY();j++) {
+ double r0=hSpec->GetYaxis()->GetBinLowEdge(j+1);
+ double r1=hSpec->GetYaxis()->GetBinUpEdge(j+1);
+ double z0=hSpec->GetXaxis()->GetBinLowEdge(i+1);
+ double z1=hSpec->GetXaxis()->GetBinUpEdge(i+1);
+ double vol=(z1-z0)*M_PI*(r1*r1-r0*r0);
+ // if |z| instead of z double the volume
+ if ( m_config.zMinFull >= 0 ) vol *= 2;
+ // if positive y hemisphere is used only -- half the volume
+ if ( m_config.posYOnly ) vol *= 0.5;
+ double val;
+ // Energy Spectra
+ for(int k=0;k<hSpec->GetNbinsZ();k++) {
+ int kBin = hSpec->GetBin(i+1,j+1,k+1);
+ int vBinE = m_config.nBinsr*nBinslogE[hi]*i+j*nBinslogE[hi]+k;
+ val = (*(pValSpec[hi]))[vBinE];
+ hSpec->SetBinContent(kBin,val/vol);
+ }
+ }
+ }
+ hSpec->Write();
+ hSpec->Delete();
+ }
+
+ std::vector<const char *> hThetaSpecNames = {
+ "theta_dphi_full_rz_neut_spec","theta_dphi_full_rz_gamm_spec","theta_dphi_full_rz_elec_spec","theta_dphi_full_rz_muon_spec","theta_dphi_full_rz_pion_spec","theta_dphi_full_rz_prot_spec","theta_dphi_full_rz_rchgd_spec","theta_dphi_full_rz_rneut_spec"
+ };
+
+ std::vector<const char *> hThetaSpecTitles = {
+ "FLUX [n(log_{10}(E/MeV))/cm^{2}]","FLUX [#gamma(log_{10}(E/MeV))/cm^{2}]","FLUX [e^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [#mu^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [#pi^{#pm}(log_{10}(E/MeV))/cm^{2}]","FLUX [p(log_{10}(E/MeV))/cm^{2}]","FLUX [rest^{chgd}(log_{10}(E/MeV))/cm^{2}]","FLUX [rest^{neut}(log_{10}(E/MeV))/cm^{2}]"
+ };
+
+ std::vector<double> thetalogEMin = {
+ m_config.logEMinn,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino,m_config.logEMino
+ };
+
+ std::vector<double> thetalogEMax = {
+ m_config.logEMaxn,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo,m_config.logEMaxo
+ };
+
+ for(unsigned int hi=0;hi<hThetaSpecNames.size();hi++) {
+ for(int i=0;i<m_config.nBinstheta;i++) {
+ for(int j=0;j<m_config.nBinsdphi;j++) {
+ TString hname(hThetaSpecNames[hi]);
+ hname += "_";
+ hname += (int)(m_config.thetaMin + i *(m_config.thetaMax-m_config.thetaMin)/m_config.nBinstheta);
+ hname += "_";
+ hname += (int)(m_config.thetaMin + (i+1)*(m_config.thetaMax-m_config.thetaMin)/m_config.nBinstheta);
+ hname += "_";
+ hname += (int)(j*360./m_config.nBinsdphi);
+ hname += "_";
+ hname += (int)((j+1)*360./m_config.nBinsdphi);
+ TH3D * hThetaSpec = new TH3D(hname.Data(),hname.Data(),m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,nBinsThetalogE[hi],thetalogEMin[hi],thetalogEMax[hi]);
+
+ hThetaSpec->SetXTitle(xtitle);
+ hThetaSpec->SetYTitle("r [cm]");
+ hThetaSpec->SetZTitle("log_{10}(E/MeV)");
+ hname = TString(hThetaSpecTitles[hi]);
+ hname += " Theta: ";
+ hname += (int)(m_config.thetaMin + i *(m_config.thetaMax-m_config.thetaMin)/m_config.nBinstheta);
+ hname += " - ";
+ hname += (int)(m_config.thetaMin + (i+1)*(m_config.thetaMax-m_config.thetaMin)/m_config.nBinstheta);
+ hname += " DPhi: ";
+ hname += (int)(j*360./m_config.nBinsdphi);
+ hname += " - ";
+ hname += (int)((j+1)*360./m_config.nBinsdphi);
+ hThetaSpec->SetTitle(hname.Data());
+ // normalize to volume element per bin
+ for(int k=0;k<hThetaSpec->GetNbinsX();k++) {
+ for(int l=0;l<hThetaSpec->GetNbinsY();l++) {
+ double r0=hThetaSpec->GetYaxis()->GetBinLowEdge(l+1);
+ double r1=hThetaSpec->GetYaxis()->GetBinUpEdge(l+1);
+ double z0=hThetaSpec->GetXaxis()->GetBinLowEdge(k+1);
+ double z1=hThetaSpec->GetXaxis()->GetBinUpEdge(k+1);
+ double vol=(z1-z0)*M_PI*(r1*r1-r0*r0);
+ // if |z| instead of z double the volume
+ if ( m_config.zMinFull >= 0 ) vol *= 2;
+ // if positive y hemisphere is used only -- half the volume
+ if ( m_config.posYOnly ) vol *= 0.5;
+ double val;
+ // Energy Spectra
+ for(int m=0;m<hThetaSpec->GetNbinsZ();m++) {
+ int mBin = hThetaSpec->GetBin(k+1,l+1,m+1);
+ int vBinETh = m_config.nBinsr*nBinsThetalogE[hi]*m_config.nBinsdphi*m_config.nBinstheta*k+l*nBinsThetalogE[hi]*m_config.nBinsdphi*m_config.nBinstheta+m*m_config.nBinsdphi*m_config.nBinstheta+i*m_config.nBinsdphi+j;
+ val = (*(pValThetaSpec[hi]))[vBinETh];
+ hThetaSpec->SetBinContent(mBin,val/vol);
+ }
+ }
+ }
+ hThetaSpec->Write();
+ hThetaSpec->Delete();
+ }
+ }
+ }
+
+ // time dependent TID maps zoom
TH3D * h_rz_tid_time = new TH3D("rz_tid_time" ,"rz_tid_time" ,m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinslogT,m_config.logTMin,m_config.logTMax);
- TH3D * h_full_rz_tid_time = new TH3D("full_rz_tid_time" ,"full_rz_tid_time" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogT,m_config.logTMin,m_config.logTMax);
h_rz_tid_time ->SetXTitle(xtitle);
- h_full_rz_tid_time->SetXTitle(xtitle);
h_rz_tid_time ->SetYTitle("r [cm]");
- h_full_rz_tid_time->SetYTitle("r [cm]");
h_rz_tid_time ->SetZTitle("log_{10}(t_{cut}/s)");
- h_full_rz_tid_time->SetZTitle("log_{10}(t_{cut}/s)");
- // element mass fraction maps
+ // element mass fraction maps zoom
TH3D * h_rz_element = new TH3D("rz_element" ,"rz_element" ,m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom,m_config.elemZMax-m_config.elemZMin+1,m_config.elemZMin-0.5,m_config.elemZMax+0.5);
- TH3D * h_full_rz_element = new TH3D("full_rz_element" ,"full_rz_element" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.elemZMax-m_config.elemZMin+1,m_config.elemZMin-0.5,m_config.elemZMax+0.5);
h_rz_element ->SetXTitle(xtitle);
- h_full_rz_element->SetXTitle(xtitle);
h_rz_element ->SetYTitle("r [cm]");
- h_full_rz_element->SetYTitle("r [cm]");
h_rz_element ->SetZTitle("Z");
- h_full_rz_element->SetZTitle("Z");
-
- TH2D * h_rz_vol = 0;
- TH2D * h_rz_norm = 0;
- TH2D * h_full_rz_vol = 0;
- TH2D * h_full_rz_norm = 0;
- TH3D * h_3d_vol = 0;
- TH3D * h_3d_norm = 0;
- if (!m_config.material.empty()) {
- // need volume fraction only if particular material is selected
- //
- // the maps for TID, NIEL and SEE need to be divided by the ratio of (vol/norm) in order to get
- // the proper estimate per volume bin for the selected material.
- // This is *not* done in the tool directly and left to the user after having summed the histograms
- // from many individual jobs.
- //
- h_rz_vol = new TH2D("rz_vol" ,"rz_vol" ,m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- h_rz_norm = new TH2D("rz_norm","rz_norm",m_config.nBinsz,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr,m_config.rMinZoom,m_config.rMaxZoom);
- h_full_rz_vol = new TH2D("full_rz_vol" ,"full_rz_vol" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- h_full_rz_norm = new TH2D("full_rz_norm","full_rz_norm",m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull);
- h_3d_vol = new TH3D("h3d_vol" ,"h3d_vol" ,m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
- h_3d_norm = new TH3D("h3d_norm","h3d_norm",m_config.nBinsz3d,m_config.zMinZoom,m_config.zMaxZoom,m_config.nBinsr3d,m_config.rMinZoom,m_config.rMaxZoom,m_config.nBinsphi3d,m_config.phiMinZoom,m_config.phiMaxZoom);
-
- h_rz_vol ->SetXTitle(xtitle);
- h_rz_norm ->SetXTitle(xtitle);
- h_rz_vol ->SetYTitle("r [cm]");
- h_rz_norm ->SetYTitle("r [cm]");
- std::string hname("Volume fraction of ");
- hname += m_config.material;
- h_rz_vol ->SetZTitle(hname.data());
- h_rz_norm ->SetZTitle("Volume norm");
-
- h_full_rz_vol ->SetXTitle(xtitle);
- h_full_rz_norm ->SetXTitle(xtitle);
- h_full_rz_vol ->SetYTitle("r [cm]");
- h_full_rz_norm ->SetYTitle("r [cm]");
- h_full_rz_vol ->SetZTitle(hname.data());
- h_full_rz_norm ->SetZTitle("Volume norm");
-
- h_3d_vol ->SetXTitle(xtitle);
- h_3d_norm ->SetXTitle(xtitle);
- h_3d_vol ->SetYTitle("r [cm]");
- h_3d_norm ->SetYTitle("r [cm]");
- h_3d_vol ->SetZTitle("#phi [#circ]");
- h_3d_norm ->SetZTitle("#phi [#circ]");
- h_3d_vol ->SetTitle(hname.data());
- h_3d_norm ->SetTitle("Volume norm");
- }
// normalize to volume element per bin
- for(int i=0;i<h_rz_tid->GetNbinsX();i++) {
- for(int j=0;j<h_rz_tid->GetNbinsY();j++) {
- int iBin = h_rz_tid->GetBin(i+1,j+1);
- int vBin = m_config.nBinsr*i+j;
- double r0=h_rz_tid->GetYaxis()->GetBinLowEdge(j+1);
- double r1=h_rz_tid->GetYaxis()->GetBinUpEdge(j+1);
- double z0=h_rz_tid->GetXaxis()->GetBinLowEdge(i+1);
- double z1=h_rz_tid->GetXaxis()->GetBinUpEdge(i+1);
+ for(int i=0;i<h_rz_tid_time->GetNbinsX();i++) {
+ for(int j=0;j<h_rz_tid_time->GetNbinsY();j++) {
+ double r0=h_rz_tid_time->GetYaxis()->GetBinLowEdge(j+1);
+ double r1=h_rz_tid_time->GetYaxis()->GetBinUpEdge(j+1);
+ double z0=h_rz_tid_time->GetXaxis()->GetBinLowEdge(i+1);
+ double z1=h_rz_tid_time->GetXaxis()->GetBinUpEdge(i+1);
double vol=(z1-z0)*M_PI*(r1*r1-r0*r0);
// if |z| instead of z double the volume
if ( m_config.zMinFull >= 0 ) vol *= 2;
// if positive y hemisphere is used only -- half the volume
if ( m_config.posYOnly ) vol *= 0.5;
double val;
- // TID
- val =maps.m_rz_tid[vBin];
- h_rz_tid->SetBinContent(iBin,val/vol);
- // EION
- val =maps.m_rz_eion[vBin];
- h_rz_eion->SetBinContent(iBin,val/vol);
- // NIEL
- val =maps.m_rz_niel[vBin];
- h_rz_niel->SetBinContent(iBin,val/vol);
- // SEE
- val =maps.m_rz_h20[vBin];
- h_rz_h20->SetBinContent(iBin,val/vol);
- // NEUT
- val =maps.m_rz_neut[vBin];
- h_rz_neut->SetBinContent(iBin,val/vol);
- // CHAD
- val =maps.m_rz_chad[vBin];
- h_rz_chad->SetBinContent(iBin,val/vol);
- // Neutron Energy Spectra
- for(int k=0;k<h_rz_neut_spec->GetNbinsZ();k++) {
- int kBin = h_rz_neut_spec->GetBin(i+1,j+1,k+1);
- int vBinE = m_config.nBinsr*m_config.nBinslogEn*i+j*m_config.nBinslogEn+k;
- val =maps.m_rz_neut_spec[vBinE];
- h_rz_neut_spec->SetBinContent(kBin,val/vol);
- }
- // All other particle's Energy Spectra
- for(int k=0;k<h_rz_gamm_spec->GetNbinsZ();k++) {
- int kBin = h_rz_gamm_spec->GetBin(i+1,j+1,k+1);
- int vBinE = m_config.nBinsr*m_config.nBinslogEo*i+j*m_config.nBinslogEo+k;
- val =maps.m_rz_gamm_spec[vBinE];
- h_rz_gamm_spec->SetBinContent(kBin,val/vol);
- val =maps.m_rz_elec_spec[vBinE];
- h_rz_elec_spec->SetBinContent(kBin,val/vol);
- val =maps.m_rz_muon_spec[vBinE];
- h_rz_muon_spec->SetBinContent(kBin,val/vol);
- val =maps.m_rz_pion_spec[vBinE];
- h_rz_pion_spec->SetBinContent(kBin,val/vol);
- val =maps.m_rz_prot_spec[vBinE];
- h_rz_prot_spec->SetBinContent(kBin,val/vol);
- val =maps.m_rz_rest_spec[vBinE];
- h_rz_rest_spec->SetBinContent(kBin,val/vol);
- }
// Time dependent TID maps
for(int k=0;k<h_rz_tid_time->GetNbinsZ();k++) {
int kBin = h_rz_tid_time->GetBin(i+1,j+1,k+1);
@@ -553,90 +613,39 @@ namespace G4UA{
val =maps.m_rz_element[vBinElem];
h_rz_element->SetBinContent(kBin,val/vol);
}
- if (!m_config.material.empty()) {
- // need volume fraction only if particular material is selected
- // VOL
- val =maps.m_rz_vol[vBin];
- h_rz_vol->SetBinContent(iBin,val/vol);
- // NORM
- val =maps.m_rz_norm[vBin];
- h_rz_norm->SetBinContent(iBin,val/vol);
- }
}
}
- h_rz_tid->Write();
- h_rz_eion->Write();
- h_rz_niel->Write();
- h_rz_h20->Write();
- h_rz_neut->Write();
- h_rz_chad->Write();
- h_rz_neut_spec->Write();
- h_rz_gamm_spec->Write();
- h_rz_elec_spec->Write();
- h_rz_muon_spec->Write();
- h_rz_pion_spec->Write();
- h_rz_prot_spec->Write();
- h_rz_rest_spec->Write();
+
h_rz_tid_time->Write();
+ h_rz_tid_time->Delete();
h_rz_element->Write();
+ h_rz_element->Delete();
+ // time dependent TID maps full
+ TH3D * h_full_rz_tid_time = new TH3D("full_rz_tid_time" ,"full_rz_tid_time" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.nBinslogT,m_config.logTMin,m_config.logTMax);
+ h_full_rz_tid_time->SetXTitle(xtitle);
+ h_full_rz_tid_time->SetYTitle("r [cm]");
+ h_full_rz_tid_time->SetZTitle("log_{10}(t_{cut}/s)");
+
+ // element mass fraction maps full
+ TH3D * h_full_rz_element = new TH3D("full_rz_element" ,"full_rz_element" ,m_config.nBinsz,m_config.zMinFull,m_config.zMaxFull,m_config.nBinsr,m_config.rMinFull,m_config.rMaxFull,m_config.elemZMax-m_config.elemZMin+1,m_config.elemZMin-0.5,m_config.elemZMax+0.5);
+ h_full_rz_element->SetXTitle(xtitle);
+ h_full_rz_element->SetYTitle("r [cm]");
+ h_full_rz_element->SetZTitle("Z");
+
// normalize to volume element per bin
- for(int i=0;i<h_full_rz_tid->GetNbinsX();i++) {
- for(int j=0;j<h_full_rz_tid->GetNbinsY();j++) {
- int iBin = h_full_rz_tid->GetBin(i+1,j+1);
- int vBin = m_config.nBinsr*i+j;
- double r0=h_full_rz_tid->GetYaxis()->GetBinLowEdge(j+1);
- double r1=h_full_rz_tid->GetYaxis()->GetBinUpEdge(j+1);
- double z0=h_full_rz_tid->GetXaxis()->GetBinLowEdge(i+1);
- double z1=h_full_rz_tid->GetXaxis()->GetBinUpEdge(i+1);
+ for(int i=0;i<h_full_rz_tid_time->GetNbinsX();i++) {
+ for(int j=0;j<h_full_rz_tid_time->GetNbinsY();j++) {
+ double r0=h_full_rz_tid_time->GetYaxis()->GetBinLowEdge(j+1);
+ double r1=h_full_rz_tid_time->GetYaxis()->GetBinUpEdge(j+1);
+ double z0=h_full_rz_tid_time->GetXaxis()->GetBinLowEdge(i+1);
+ double z1=h_full_rz_tid_time->GetXaxis()->GetBinUpEdge(i+1);
double vol=(z1-z0)*M_PI*(r1*r1-r0*r0);
// if |z| instead of z double the volume
if ( m_config.zMinFull >= 0 ) vol *= 2;
// if positive y hemisphere is used only -- half the volume
if ( m_config.posYOnly ) vol *= 0.5;
double val;
- // TID
- val =maps.m_full_rz_tid[vBin];
- h_full_rz_tid->SetBinContent(iBin,val/vol);
- // EION
- val =maps.m_full_rz_eion[vBin];
- h_full_rz_eion->SetBinContent(iBin,val/vol);
- // NIEL
- val =maps.m_full_rz_niel[vBin];
- h_full_rz_niel->SetBinContent(iBin,val/vol);
- // SEE
- val =maps.m_full_rz_h20[vBin];
- h_full_rz_h20->SetBinContent(iBin,val/vol);
- // NEUT
- val =maps.m_full_rz_neut[vBin];
- h_full_rz_neut->SetBinContent(iBin,val/vol);
- // CHAD
- val =maps.m_full_rz_chad[vBin];
- h_full_rz_chad->SetBinContent(iBin,val/vol);
- // Neutron Energy Spectra
- for(int k=0;k<h_full_rz_neut_spec->GetNbinsZ();k++) {
- int kBin = h_full_rz_neut_spec->GetBin(i+1,j+1,k+1);
- int vBinE = m_config.nBinsr*m_config.nBinslogEn*i+j*m_config.nBinslogEn+k;
- val =maps.m_full_rz_neut_spec[vBinE];
- h_full_rz_neut_spec->SetBinContent(kBin,val/vol);
- }
- // All other particle's Energy Spectra
- for(int k=0;k<h_full_rz_gamm_spec->GetNbinsZ();k++) {
- int kBin = h_full_rz_gamm_spec->GetBin(i+1,j+1,k+1);
- int vBinE = m_config.nBinsr*m_config.nBinslogEo*i+j*m_config.nBinslogEo+k;
- val =maps.m_full_rz_gamm_spec[vBinE];
- h_full_rz_gamm_spec->SetBinContent(kBin,val/vol);
- val =maps.m_full_rz_elec_spec[vBinE];
- h_full_rz_elec_spec->SetBinContent(kBin,val/vol);
- val =maps.m_full_rz_muon_spec[vBinE];
- h_full_rz_muon_spec->SetBinContent(kBin,val/vol);
- val =maps.m_full_rz_pion_spec[vBinE];
- h_full_rz_pion_spec->SetBinContent(kBin,val/vol);
- val =maps.m_full_rz_prot_spec[vBinE];
- h_full_rz_prot_spec->SetBinContent(kBin,val/vol);
- val =maps.m_full_rz_rest_spec[vBinE];
- h_full_rz_rest_spec->SetBinContent(kBin,val/vol);
- }
// Time dependent TID maps
for(int k=0;k<h_full_rz_tid_time->GetNbinsZ();k++) {
int kBin = h_full_rz_tid_time->GetBin(i+1,j+1,k+1);
@@ -651,103 +660,12 @@ namespace G4UA{
val =maps.m_full_rz_element[vBinElem];
h_full_rz_element->SetBinContent(kBin,val/vol);
}
- if (!m_config.material.empty()) {
- // need volume fraction only if particular material is selected
- // VOL
- val =maps.m_full_rz_vol[vBin];
- h_full_rz_vol->SetBinContent(iBin,val/vol);
- // NORM
- val =maps.m_full_rz_norm[vBin];
- h_full_rz_norm->SetBinContent(iBin,val/vol);
- }
}
}
- h_full_rz_tid->Write();
- h_full_rz_eion->Write();
- h_full_rz_niel->Write();
- h_full_rz_h20->Write();
- h_full_rz_neut->Write();
- h_full_rz_chad->Write();
- h_full_rz_neut_spec->Write();
- h_full_rz_gamm_spec->Write();
- h_full_rz_elec_spec->Write();
- h_full_rz_muon_spec->Write();
- h_full_rz_pion_spec->Write();
- h_full_rz_prot_spec->Write();
- h_full_rz_rest_spec->Write();
h_full_rz_tid_time->Write();
+ h_full_rz_tid_time->Delete();
h_full_rz_element->Write();
-
- // normalize to volume element per bin
- for(int i=0;i<h_3d_tid->GetNbinsX();i++) { /* |z| */
- for(int j=0;j<h_3d_tid->GetNbinsY();j++) { /* r */
- for(int k=0;k<h_3d_tid->GetNbinsZ();k++) { /* phi */
- int vBin = m_config.nBinsr3d*m_config.nBinsphi3d*i+m_config.nBinsphi3d*j+k;
- int iBin = h_3d_tid->GetBin(i+1,j+1,k+1);
- double phi0=h_3d_tid->GetZaxis()->GetBinLowEdge(k+1);
- double phi1=h_3d_tid->GetZaxis()->GetBinUpEdge(k+1);
- double r0=h_3d_tid->GetYaxis()->GetBinLowEdge(j+1);
- double r1=h_3d_tid->GetYaxis()->GetBinUpEdge(j+1);
- double z0=h_3d_tid->GetXaxis()->GetBinLowEdge(i+1);
- double z1=h_3d_tid->GetXaxis()->GetBinUpEdge(i+1);
- double vol=(z1-z0)*M_PI*(r1*r1-r0*r0)*(phi1-phi0)/360.;
- // if |z| instead of z double the volume
- if ( m_config.zMinFull >= 0 ) vol *= 2;
- // assume that phi-range corresponds to full 360 degrees in case
- // lower phi boundary is 0 - i.e. all phi-segments mapped to first
- if ( m_config.phiMinZoom == 0 ) {
- vol *= 360./m_config.phiMaxZoom;
- // if positive y hemisphere is used only -- half the volume
- if ( m_config.posYOnly )
- vol *= 0.5;
- }
- double val;
- // TID
- val =maps.m_3d_tid[vBin];
- h_3d_tid->SetBinContent(iBin,val/vol);
- // EION
- val =maps.m_3d_eion[vBin];
- h_3d_eion->SetBinContent(iBin,val/vol);
- // NIEL
- val =maps.m_3d_niel[vBin];
- h_3d_niel->SetBinContent(iBin,val/vol);
- // SEE
- val =maps.m_3d_h20[vBin];
- h_3d_h20->SetBinContent(iBin,val/vol);
- // NEUT
- val =maps.m_3d_neut[vBin];
- h_3d_neut->SetBinContent(iBin,val/vol);
- // CHAD
- val =maps.m_3d_chad[vBin];
- h_3d_chad->SetBinContent(iBin,val/vol);
- if (!m_config.material.empty()) {
- // need volume fraction only if particular material is selected
- // VOL
- val =maps.m_3d_vol[vBin];
- h_3d_vol->SetBinContent(iBin,val/vol);
- // NORM
- val =maps.m_3d_norm[vBin];
- h_3d_norm->SetBinContent(iBin,val/vol);
- }
- }
- }
- }
- h_3d_tid->Write();
- h_3d_eion->Write();
- h_3d_niel->Write();
- h_3d_h20->Write();
- h_3d_neut->Write();
- h_3d_chad->Write();
-
- if (!m_config.material.empty()) {
- // need volume fraction only if particular material is selected
- h_rz_vol->Write();
- h_rz_norm->Write();
- h_full_rz_vol->Write();
- h_full_rz_norm->Write();
- h_3d_vol->Write();
- h_3d_norm->Write();
- }
+ h_full_rz_element->Delete();
f->Close();
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.h b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.h
index 6313fafd75a..8d0cb87f4d2 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.h
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#ifndef G4USERACTIONS_G4UA__RADIATIONMAPSMAKERTOOL_H
--
GitLab