From 007a98d6965e3558aa2ed46511155e2cfd5affd5 Mon Sep 17 00:00:00 2001
From: Sven Menke <sven.menke@cern.ch>
Date: Sat, 27 Apr 2019 10:51:45 +0200
Subject: [PATCH] latest 21.0 changes to RadiationMapsMaker (particle spectra)
ported to master
---
.../G4UserActions/RadiationMapsMaker.h | 61 +++
.../share/jobOptions.RadiationMapsMaker.py | 6 +
.../G4UserActions/src/RadiationMapsMaker.cxx | 454 +++++++++++-------
.../src/RadiationMapsMakerTool.cxx | 235 ++++++++-
4 files changed, 578 insertions(+), 178 deletions(-)
diff --git a/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h b/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
index 34dfd248d8b..f5e266af515 100644
--- a/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
+++ b/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
@@ -51,6 +51,16 @@ namespace G4UA
double phiMinZoom = -180.; // degrees
double phiMaxZoom = 180.; // degrees
+
+ // neutron spectra
+ int nBinslogEn = 90;
+ double logEMinn = -11.; // min log10(E_kin/MeV)
+ double logEMaxn = 7.; // max log10(E_kin/MeV)
+
+ // particle spectra for gamma,e^+/-,mu^+/-,pi^+/-,p,rest
+ int nBinslogEo = 45;
+ double logEMino = -2.; // min log10(E_kin/MeV)
+ double logEMaxo = 7.; // max log10(E_kin/MeV)
};
@@ -117,6 +127,57 @@ namespace G4UA
/// vector to normalize the volume fraction in 3d
std::vector<double> m_3d_norm;
+ // particle spectra
+
+ // neutrons
+
+ /// vector of neutron spectra in log10(E/MeV) bins and the zoom 2d grid
+ 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;
+
+ // 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
+ std::vector<double> m_full_rz_gamm_spec;
+
+ // e^+/-
+
+ /// vector of e^+/- spectra in log10(E/MeV) bins and the zoom 2d grid
+ 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;
+
+ // mu^+/-
+
+ /// vector of mu^+/- spectra in log10(E/MeV) bins and the zoom 2d grid
+ 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;
+
+ // pi^+/-
+
+ /// vector of pi^+/- spectra in log10(E/MeV) bins and the zoom 2d grid
+ 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;
+
+ // proton
+
+ /// vector of proton spectra in log10(E/MeV) bins and the zoom 2d grid
+ 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;
+
+ // rest
+
+ /// vector of other particle spectra in log10(E/MeV) bins and the zoom 2d grid
+ 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;
+
void merge(const Report& maps);
};
diff --git a/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py b/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
index dcc4ffcae05..8636a2680da 100644
--- a/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
+++ b/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
@@ -8,6 +8,8 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
# radmaptool.Material = "" # if left empty all materials are used (default)
# radmaptool.NBinsR = 120
# radmaptool.NBinsZ = 240
+# radmaptool.NBinsLogEn = 90
+# radmaptool.NBinsLogEo = 45
# radmaptool.NBinsR3D = 30
# radmaptool.NBinsZ3D = 60
# radmaptool.NBinsPhi3D = 32
@@ -21,6 +23,10 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
# radmaptool.ZMaxFull = 2400.0 # in cm
# radmaptool.PhiMinZoom = -180.0 # in degrees
# radmaptool.PhiMaxZoom = 180.0 # in degrees
+# radmaptool.LogEMinn = -11.0 # in log10(E/MeV)
+# radmaptool.LogEMaxn = 7.0 # in log10(E/MeV)
+# radmaptool.LogEMino = -2.0 # in log10(E/MeV)
+# radmaptool.LogEMaxo = 7.0 # in log10(E/MeV)
#
simFlags.OptionalUserActionList.addAction('G4UA::RadiationMapsMakerTool')
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
index b58aeb2b8fc..4e791a49b6e 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
@@ -110,6 +110,27 @@ namespace G4UA{
m_3d_vol [i] += maps.m_3d_vol [i];
m_3d_norm[i] += maps.m_3d_norm[i];
}
+
+ // neutron spectra have different size from all other particle's spectra
+ for(unsigned int i=0;i<maps.m_rz_neut_spec.size();i++) {
+ m_rz_neut_spec [i] += maps.m_rz_neut_spec [i];
+ m_full_rz_neut_spec[i] += maps.m_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];
+ m_full_rz_gamm_spec[i] += maps.m_full_rz_gamm_spec[i];
+ m_rz_elec_spec [i] += maps.m_rz_elec_spec [i];
+ m_full_rz_elec_spec[i] += maps.m_full_rz_elec_spec[i];
+ m_rz_muon_spec [i] += maps.m_rz_muon_spec [i];
+ m_full_rz_muon_spec[i] += maps.m_full_rz_muon_spec[i];
+ m_rz_pion_spec [i] += maps.m_rz_pion_spec [i];
+ m_full_rz_pion_spec[i] += maps.m_full_rz_pion_spec[i];
+ m_rz_prot_spec [i] += maps.m_rz_prot_spec [i];
+ m_full_rz_prot_spec[i] += maps.m_full_rz_prot_spec[i];
+ m_rz_rest_spec [i] += maps.m_rz_rest_spec [i];
+ m_full_rz_rest_spec[i] += maps.m_full_rz_rest_spec[i];
+ }
}
void RadiationMapsMaker::BeginOfRunAction(const G4Run*){
@@ -137,6 +158,21 @@ namespace G4UA{
m_maps.m_3d_neut.resize(0);
m_maps.m_3d_chad.resize(0);
+ m_maps.m_rz_neut_spec .resize(0);
+ m_maps.m_full_rz_neut_spec.resize(0);
+ m_maps.m_rz_gamm_spec .resize(0);
+ m_maps.m_full_rz_gamm_spec.resize(0);
+ m_maps.m_rz_elec_spec .resize(0);
+ m_maps.m_full_rz_elec_spec.resize(0);
+ m_maps.m_rz_muon_spec .resize(0);
+ m_maps.m_full_rz_muon_spec.resize(0);
+ m_maps.m_rz_pion_spec .resize(0);
+ m_maps.m_full_rz_pion_spec.resize(0);
+ m_maps.m_rz_prot_spec .resize(0);
+ 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);
+
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
m_maps.m_rz_vol .resize(0);
@@ -170,6 +206,21 @@ namespace G4UA{
m_maps.m_3d_neut.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
m_maps.m_3d_chad.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
+ m_maps.m_rz_neut_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
+ m_maps.m_full_rz_neut_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEn,0.0);
+ m_maps.m_rz_gamm_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_full_rz_gamm_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_rz_elec_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_full_rz_elec_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_rz_muon_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_full_rz_muon_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_rz_pion_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_full_rz_pion_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_rz_prot_spec .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ m_maps.m_full_rz_prot_spec.resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogEo,0.0);
+ 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);
+
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// 2d zoom
@@ -276,14 +327,15 @@ namespace G4UA{
pdgid = 9;
}
}
- // process NIEL, h20, Edep, NEUT and CHAD particles only
- if ( pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 || /* NIEL & h20*/
+ // process spectra, NIEL, h20, Edep, NEUT and CHAD particles only
+
+ if ( pdgid == 0 || pdgid == 3 || /* spectra */
+ pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 || /* NIEL & h20*/
aStep->GetTotalEnergyDeposit() > 0 || pdgid == 999) {
-
double absq = fabs(aStep->GetTrack()->GetDefinition()->GetPDGCharge());
-
+
double rho = aStep->GetTrack()->GetMaterial()->GetDensity()/CLHEP::g*CLHEP::cm3;
bool goodMaterial(false);
@@ -318,6 +370,8 @@ namespace G4UA{
double weight = 0; // weight for NIEL
double eKin = aStep->GetTrack()->GetKineticEnergy();
+ double logEKin = (eKin > 0?log10(eKin):(m_config.logEMinn<m_config.logEMino?m_config.logEMinn:m_config.logEMino)-1);
+
if ( pdgid == 1 || pdgid == 9 ) {
if ( eKin < 15 ) {
if ( eKin > 10 ) {
@@ -346,195 +400,271 @@ namespace G4UA{
weight = m_tgeSi->Eval(eKin);
}
}
-
double dE_TOT = aStep->GetTotalEnergyDeposit()/nStep;
double dE_NIEL = aStep->GetNonIonizingEnergyDeposit()/nStep;
double dE_ION = dE_TOT-dE_NIEL;
-
- // NIEL SEE TID NEUT CHAD geantino
- if ( weight > 0 || eKin > 20 || dE_TOT > 0 || ((pdgid == 6 || pdgid == 7) && eKin > 0.1) || (absq>0 && (pdgid == 1 || pdgid == 2 || pdgid == 8 || pdgid == 9)) || pdgid == 999 ) {
-
- for(unsigned int i=0;i<nStep;i++) {
- double absz = fabs(z0+dz*(i+0.5));
- double rr = sqrt(pow(x0+dx*(i+0.5),2)+
- pow(y0+dy*(i+0.5),2));
- double pphi = atan2(y0+dy*(i+0.5),x0+dx*(i+0.5))*180/M_PI;
-
- int vBinZoom = -1;
- int vBinFull = -1;
- int vBin3d = -1;
-
- // zoom 2d
- if ( m_config.zMinZoom < absz &&
- m_config.zMaxZoom > absz ) {
- int iz = (absz-m_config.zMinZoom)/(m_config.zMaxZoom-m_config.zMinZoom)*m_config.nBinsz;
- if ( m_config.rMinZoom < rr &&
- m_config.rMaxZoom > rr ) {
- int ir = (rr-m_config.rMinZoom)/(m_config.rMaxZoom-m_config.rMinZoom)*m_config.nBinsr;
- vBinZoom = m_config.nBinsr*iz+ir;
+
+ for(unsigned int i=0;i<nStep;i++) {
+ double absz = fabs(z0+dz*(i+0.5));
+ double rr = sqrt(pow(x0+dx*(i+0.5),2)+
+ pow(y0+dy*(i+0.5),2));
+ double pphi = atan2(y0+dy*(i+0.5),x0+dx*(i+0.5))*180/M_PI;
+
+ int vBinZoom = -1;
+ int vBinFull = -1;
+ int vBin3d = -1;
+ int vBinZoomSpecn = -1;
+ int vBinFullSpecn = -1;
+ int vBinZoomSpeco = -1;
+ int vBinFullSpeco = -1;
+
+ // zoom 2d
+ if ( m_config.zMinZoom < absz &&
+ m_config.zMaxZoom > absz ) {
+ int iz = (absz-m_config.zMinZoom)/(m_config.zMaxZoom-m_config.zMinZoom)*m_config.nBinsz;
+ if ( m_config.rMinZoom < rr &&
+ m_config.rMaxZoom > rr ) {
+ int ir = (rr-m_config.rMinZoom)/(m_config.rMaxZoom-m_config.rMinZoom)*m_config.nBinsr;
+ vBinZoom = m_config.nBinsr*iz+ir;
+ if ( m_config.logEMinn < logEKin &&
+ m_config.logEMaxn > logEKin &&
+ (pdgid == 6 || pdgid == 7)) {
+ int ile = (logEKin-m_config.logEMinn)/(m_config.logEMaxn-m_config.logEMinn)*m_config.nBinslogEn;
+ vBinZoomSpecn = m_config.nBinsr*m_config.nBinslogEn*iz+ir*m_config.nBinslogEn+ile;
+ }
+ 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;
+ vBinZoomSpeco = m_config.nBinsr*m_config.nBinslogEo*iz+ir*m_config.nBinslogEo+ile;
}
}
-
- // full 2d
- if ( m_config.zMinFull < absz &&
- m_config.zMaxFull > absz ) {
- int iz = (absz-m_config.zMinFull)/(m_config.zMaxFull-m_config.zMinFull)*m_config.nBinsz;
- if ( m_config.rMinFull < rr &&
- m_config.rMaxFull > rr ) {
- int ir = (rr-m_config.rMinFull)/(m_config.rMaxFull-m_config.rMinFull)*m_config.nBinsr;
- vBinFull = m_config.nBinsr*iz+ir;
+ }
+
+ // full 2d
+ if ( m_config.zMinFull < absz &&
+ m_config.zMaxFull > absz ) {
+ int iz = (absz-m_config.zMinFull)/(m_config.zMaxFull-m_config.zMinFull)*m_config.nBinsz;
+ if ( m_config.rMinFull < rr &&
+ m_config.rMaxFull > rr ) {
+ int ir = (rr-m_config.rMinFull)/(m_config.rMaxFull-m_config.rMinFull)*m_config.nBinsr;
+ vBinFull = m_config.nBinsr*iz+ir;
+ if ( m_config.logEMinn < logEKin &&
+ m_config.logEMaxn > logEKin &&
+ (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.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;
}
}
-
- // zoom 3d
- if ( m_config.zMinZoom < absz &&
- m_config.zMaxZoom > absz ) {
- int iz = (absz-m_config.zMinZoom)/(m_config.zMaxZoom-m_config.zMinZoom)*m_config.nBinsz3d;
- if ( m_config.rMinZoom < rr &&
- m_config.rMaxZoom > rr ) {
- int ir = (rr-m_config.rMinZoom)/(m_config.rMaxZoom-m_config.rMinZoom)*m_config.nBinsr3d;
- if ( m_config.phiMinZoom == 0) {
- // assume that all phi should be mapped to the selected phi range
- double phi_mapped = pphi;
- // first use phi range from 0 - 360 degrees
- if (phi_mapped < 0) {
- phi_mapped = 360 + phi_mapped;
- }
- // then map to selected phi-range
- int iphi = phi_mapped/m_config.phiMaxZoom;
- phi_mapped -= iphi*m_config.phiMaxZoom;
- iphi = phi_mapped/m_config.phiMaxZoom*m_config.nBinsphi3d;
- vBin3d = m_config.nBinsr3d*m_config.nBinsphi3d*iz+m_config.nBinsphi3d*ir+iphi;
- }
- else if ( m_config.phiMinZoom < pphi &&
- m_config.phiMaxZoom > pphi ) {
- int iphi = (pphi-m_config.phiMinZoom)/(m_config.phiMaxZoom-m_config.phiMinZoom)*m_config.nBinsphi3d;
- vBin3d = m_config.nBinsr3d*m_config.nBinsphi3d*iz+m_config.nBinsphi3d*ir+iphi;
+ }
+
+ // zoom 3d
+ if ( m_config.zMinZoom < absz &&
+ m_config.zMaxZoom > absz ) {
+ int iz = (absz-m_config.zMinZoom)/(m_config.zMaxZoom-m_config.zMinZoom)*m_config.nBinsz3d;
+ if ( m_config.rMinZoom < rr &&
+ m_config.rMaxZoom > rr ) {
+ int ir = (rr-m_config.rMinZoom)/(m_config.rMaxZoom-m_config.rMinZoom)*m_config.nBinsr3d;
+ if ( m_config.phiMinZoom == 0) {
+ // assume that all phi should be mapped to the selected phi range
+ double phi_mapped = pphi;
+ // first use phi range from 0 - 360 degrees
+ if (phi_mapped < 0) {
+ phi_mapped = 360 + phi_mapped;
}
+ // then map to selected phi-range
+ int iphi = phi_mapped/m_config.phiMaxZoom;
+ phi_mapped -= iphi*m_config.phiMaxZoom;
+ iphi = phi_mapped/m_config.phiMaxZoom*m_config.nBinsphi3d;
+ vBin3d = m_config.nBinsr3d*m_config.nBinsphi3d*iz+m_config.nBinsphi3d*ir+iphi;
+ }
+ else if ( m_config.phiMinZoom < pphi &&
+ m_config.phiMaxZoom > pphi ) {
+ int iphi = (pphi-m_config.phiMinZoom)/(m_config.phiMaxZoom-m_config.phiMinZoom)*m_config.nBinsphi3d;
+ vBin3d = m_config.nBinsr3d*m_config.nBinsphi3d*iz+m_config.nBinsphi3d*ir+iphi;
}
}
-
- // TID & EION
- if ( goodMaterial && vBinZoom >=0 ) {
- if ( pdgid == 999 ) {
- m_maps.m_rz_tid [vBinZoom] += dl;
- m_maps.m_rz_eion[vBinZoom] += rho*dl;
+ }
+ // TID & EION
+ if ( goodMaterial && vBinZoom >=0 ) {
+ if ( pdgid == 999 ) {
+ m_maps.m_rz_tid [vBinZoom] += dl;
+ m_maps.m_rz_eion[vBinZoom] += rho*dl;
+ }
+ else {
+ m_maps.m_rz_tid [vBinZoom] += dE_ION/rho;
+ m_maps.m_rz_eion[vBinZoom] += dE_ION;
+ }
+ }
+ if ( goodMaterial && vBinFull >=0 ) {
+ if ( pdgid == 999 ) {
+ m_maps.m_full_rz_tid [vBinFull] += dl;
+ m_maps.m_full_rz_eion[vBinFull] += rho*dl;
+ }
+ else {
+ m_maps.m_full_rz_tid [vBinFull] += dE_ION/rho;
+ m_maps.m_full_rz_eion[vBinFull] += dE_ION;
+ }
+ }
+ if ( goodMaterial && vBin3d >=0 ) {
+ if ( pdgid == 999 ) {
+ m_maps.m_3d_tid [vBin3d] += dl;
+ m_maps.m_3d_eion[vBin3d] += rho*dl;
+ }
+ else {
+ m_maps.m_3d_tid [vBin3d] += dE_ION/rho;
+ m_maps.m_3d_eion[vBin3d] += dE_ION;
+ }
+ }
+
+ if ( goodMaterial && (pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 )) {
+ // NIEL
+ if ( weight > 0 ) {
+ if ( vBinZoom >=0 ) {
+ m_maps.m_rz_niel [vBinZoom] += weight*dl;
}
- else {
- m_maps.m_rz_tid [vBinZoom] += dE_ION/rho;
- m_maps.m_rz_eion[vBinZoom] += dE_ION;
+ if ( vBinFull >=0 ) {
+ m_maps.m_full_rz_niel [vBinFull] += weight*dl;
+ }
+ if ( vBin3d >=0 ) {
+ m_maps.m_3d_niel [vBin3d] += weight*dl;
}
}
- if ( goodMaterial && vBinFull >=0 ) {
- if ( pdgid == 999 ) {
- m_maps.m_full_rz_tid [vBinFull] += dl;
- m_maps.m_full_rz_eion[vBinFull] += rho*dl;
+ // SEE
+ if ( eKin > 20 && (pdgid == 1 || pdgid == 2 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9) ) {
+ if ( vBinZoom >=0 ) {
+ m_maps.m_rz_h20 [vBinZoom] += dl;
}
- else {
- m_maps.m_full_rz_tid [vBinFull] += dE_ION/rho;
- m_maps.m_full_rz_eion[vBinFull] += dE_ION;
+ if ( vBinFull >=0 ) {
+ m_maps.m_full_rz_h20 [vBinFull] += dl;
+ }
+ if ( vBin3d >=0 ) {
+ m_maps.m_3d_h20 [vBin3d] += dl;
+ }
+ }
+ // NEUT
+ if ( eKin > 0.1 && (pdgid == 6 || pdgid == 7 ) ) {
+ if ( vBinZoom >=0 ) {
+ m_maps.m_rz_neut [vBinZoom] += dl;
+ }
+ if ( vBinFull >=0 ) {
+ m_maps.m_full_rz_neut [vBinFull] += dl;
+ }
+ if ( vBin3d >=0 ) {
+ m_maps.m_3d_neut [vBin3d] += dl;
+ }
+ }
+ // CHAD
+ if ( absq > 0 && (pdgid == 1 || pdgid == 2 || pdgid == 8 || pdgid == 9 ) ) {
+ if ( vBinZoom >=0 ) {
+ m_maps.m_rz_chad [vBinZoom] += dl;
+ }
+ if ( vBinFull >=0 ) {
+ m_maps.m_full_rz_chad [vBinFull] += dl;
+ }
+ if ( vBin3d >=0 ) {
+ m_maps.m_3d_chad [vBin3d] += dl;
+ }
+ }
+ // Neutron Energy Spectra
+ if ( pdgid == 6 || pdgid == 7 ) {
+ if ( vBinZoomSpecn >=0 ) {
+ m_maps.m_rz_neut_spec [vBinZoomSpecn] += dl;
+ }
+ if ( vBinFullSpecn >=0 ) {
+ m_maps.m_full_rz_neut_spec [vBinFullSpecn] += dl;
}
}
- if ( goodMaterial && vBin3d >=0 ) {
- if ( pdgid == 999 ) {
- m_maps.m_3d_tid [vBin3d] += dl;
- m_maps.m_3d_eion[vBin3d] += rho*dl;
+ }
+
+ if ( goodMaterial && (pdgid < 6 || pdgid > 7 )){
+ // Other particle Energy Spectra
+ if ( vBinZoomSpeco >=0 ) {
+ if ( pdgid == 0 ) {
+ m_maps.m_rz_gamm_spec [vBinZoomSpeco] += dl;
+ }
+ else if ( pdgid == 1 ) {
+ m_maps.m_rz_prot_spec [vBinZoomSpeco] += dl;
+ }
+ else if ( pdgid == 2 ) {
+ m_maps.m_rz_pion_spec [vBinZoomSpeco] += dl;
+ }
+ else if ( pdgid == 3 ) {
+ m_maps.m_rz_muon_spec [vBinZoomSpeco] += dl;
+ }
+ else if ( pdgid == 4 || pdgid == 5 ) {
+ m_maps.m_rz_elec_spec [vBinZoomSpeco] += dl;
}
else {
- m_maps.m_3d_tid [vBin3d] += dE_ION/rho;
- m_maps.m_3d_eion[vBin3d] += dE_ION;
+ m_maps.m_rz_rest_spec [vBinZoomSpeco] += dl;
}
}
-
- if ( goodMaterial && (pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 )) {
- // NIEL
- if ( weight > 0 ) {
- if ( vBinZoom >=0 ) {
- m_maps.m_rz_niel [vBinZoom] += weight*dl;
- }
- if ( vBinFull >=0 ) {
- m_maps.m_full_rz_niel [vBinFull] += weight*dl;
- }
- if ( vBin3d >=0 ) {
- m_maps.m_3d_niel [vBin3d] += weight*dl;
- }
+ if ( vBinFullSpeco >=0 ) {
+ if ( pdgid == 0 ) {
+ m_maps.m_full_rz_gamm_spec [vBinFullSpeco] += dl;
}
- // SEE
- if ( eKin > 20 && (pdgid == 1 || pdgid == 2 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9) ) {
- if ( vBinZoom >=0 ) {
- m_maps.m_rz_h20 [vBinZoom] += dl;
- }
- if ( vBinFull >=0 ) {
- m_maps.m_full_rz_h20 [vBinFull] += dl;
- }
- if ( vBin3d >=0 ) {
- m_maps.m_3d_h20 [vBin3d] += dl;
- }
+ else if ( pdgid == 1 ) {
+ m_maps.m_full_rz_prot_spec [vBinFullSpeco] += dl;
}
- // NEUT
- if ( eKin > 0.1 && (pdgid == 6 || pdgid == 7 ) ) {
- if ( vBinZoom >=0 ) {
- m_maps.m_rz_neut [vBinZoom] += dl;
- }
- if ( vBinFull >=0 ) {
- m_maps.m_full_rz_neut [vBinFull] += dl;
- }
- if ( vBin3d >=0 ) {
- m_maps.m_3d_neut [vBin3d] += dl;
- }
+ else if ( pdgid == 2 ) {
+ m_maps.m_full_rz_pion_spec [vBinFullSpeco] += dl;
}
- // CHAD
- if ( absq > 0 && (pdgid == 1 || pdgid == 2 || pdgid == 8 || pdgid == 9 ) ) {
- if ( vBinZoom >=0 ) {
- m_maps.m_rz_chad [vBinZoom] += dl;
- }
- if ( vBinFull >=0 ) {
- m_maps.m_full_rz_chad [vBinFull] += dl;
- }
- if ( vBin3d >=0 ) {
- m_maps.m_3d_chad [vBin3d] += dl;
- }
+ else if ( pdgid == 3 ) {
+ m_maps.m_full_rz_muon_spec [vBinFullSpeco] += dl;
+ }
+ else if ( pdgid == 4 || pdgid == 5 ) {
+ m_maps.m_full_rz_elec_spec [vBinFullSpeco] += dl;
+ }
+ else {
+ m_maps.m_full_rz_rest_spec [vBinFullSpeco] += dl;
}
}
- if (!m_config.material.empty()) {
- // need volume fraction only if particular material is selected
- if ( (eKin > 1 && (pdgid == 6 || pdgid == 7)) || pdgid == 999) {
- // count all neutron > 1 MeV track lengths weighted by r
- // to get norm for volume per bin. High energetic
- // neutrons are used because they travel far enough to
- // map entire bins and are not bent by magnetic fields.
- // dl is a measure of length inside the current bin.
- // The multiplication by r accounts for the larger
- // volume corresponding to larger r assuming that the
- // neutron flux is locally mainly from inside to
- // outside. In regions where the neutron flux differs
- // substantially from this cylindrical assumption a
- // cylindrical Geantino scan (vertex: flat in z, x=y=0;
- // momentum: pz=0, flat in phi) should be used to get
- // the correct volume fraction.
+ }
+
+ if (!m_config.material.empty()) {
+ // need volume fraction only if particular material is selected
+ if ( (eKin > 1 && (pdgid == 6 || pdgid == 7)) || pdgid == 999) {
+ // count all neutron > 1 MeV track lengths weighted by r
+ // to get norm for volume per bin. High energetic
+ // neutrons are used because they travel far enough to
+ // map entire bins and are not bent by magnetic fields.
+ // dl is a measure of length inside the current bin.
+ // The multiplication by r accounts for the larger
+ // volume corresponding to larger r assuming that the
+ // neutron flux is locally mainly from inside to
+ // outside. In regions where the neutron flux differs
+ // substantially from this cylindrical assumption a
+ // cylindrical Geantino scan (vertex: flat in z, x=y=0;
+ // momentum: pz=0, flat in phi) should be used to get
+ // the correct volume fraction.
+ if ( vBinZoom >=0 ) {
+ m_maps.m_rz_norm[vBinZoom] += rr*dl;
+ }
+ if ( vBinFull >=0 ) {
+ m_maps.m_full_rz_norm[vBinFull] += rr*dl;
+ }
+ if ( vBin3d >=0 ) {
+ m_maps.m_3d_norm[vBin3d] += rr*dl;
+ }
+ if ( goodMaterial ) {
+ // same but only inside the material of interest.
+ // The ratio vol/norm gives the volume fraction of the desired
+ // material inside the current bin
if ( vBinZoom >=0 ) {
- m_maps.m_rz_norm[vBinZoom] += rr*dl;
+ m_maps.m_rz_vol [vBinZoom] += rr*dl;
}
if ( vBinFull >=0 ) {
- m_maps.m_full_rz_norm[vBinFull] += rr*dl;
+ m_maps.m_full_rz_vol [vBinFull] += rr*dl;
}
if ( vBin3d >=0 ) {
- m_maps.m_3d_norm[vBin3d] += rr*dl;
- }
- if ( goodMaterial ) {
- // same but only inside the material of interest.
- // The ratio vol/norm gives the volume fraction of the desired
- // material inside the current bin
- if ( vBinZoom >=0 ) {
- m_maps.m_rz_vol [vBinZoom] += rr*dl;
- }
- if ( vBinFull >=0 ) {
- m_maps.m_full_rz_vol [vBinFull] += rr*dl;
- }
- if ( vBin3d >=0 ) {
- m_maps.m_3d_vol [vBin3d] += rr*dl;
- }
+ m_maps.m_3d_vol [vBin3d] += rr*dl;
}
}
}
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
index de6c9f17a0f..9015b2f8c1f 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
@@ -22,6 +22,10 @@ namespace G4UA{
/// number of bins in r and z for all 2D maps
declareProperty("NBinsR" , m_config.nBinsr);
declareProperty("NBinsZ" , m_config.nBinsz);
+ /// number of bins in logE for energy spectra of neutrons in 2D grids
+ 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 r, z and phi for all 3D maps
declareProperty("NBinsR3D" , m_config.nBinsr3d);
declareProperty("NBinsZ3D" , m_config.nBinsz3d);
@@ -40,6 +44,12 @@ namespace G4UA{
/// for Zoomed area in 3D
declareProperty("PhiMinZoom", m_config.phiMinZoom);
declareProperty("PhiMaxZoom", m_config.phiMaxZoom);
+ /// for logE of neutrons in 2D spectra
+ declareProperty("LogEMinn" , m_config.logEMinn);
+ declareProperty("LogEMaxn" , m_config.logEMaxn);
+ /// for logE of other particles in 2D spectra
+ declareProperty("LogEMino" , m_config.logEMino);
+ declareProperty("LogEMaxo" , m_config.logEMaxo);
}
//---------------------------------------------------------------------------
@@ -47,21 +57,25 @@ namespace G4UA{
//---------------------------------------------------------------------------
StatusCode RadiationMapsMakerTool::initialize()
{
- ATH_MSG_INFO( "Initializing " << name() << "\n" <<
- "OutputFile: " << m_radMapsFileName << "\n" <<
- "Material: " << m_config.material << "\n" <<
- "2D Maps: " << m_config.nBinsz << " |z|-bins, " <<
- m_config.nBinsr << " r-bins" << "\n" <<
- "Zoom: " << m_config.zMinZoom << " < |z|/cm < " << m_config.zMaxZoom << ", " <<
- m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << "\n" <<
- "Full: " << m_config.zMinFull << " < |z|/cm < " << m_config.zMaxFull << ", " <<
- m_config.rMinFull << " < r/cm < " << m_config.rMaxFull << "\n" <<
- "3D Maps: " << m_config.nBinsz3d << " |z|-bins, " <<
- m_config.nBinsr3d << " r-bins, " <<
- m_config.nBinsphi3d << " phi-bins" << "\n" <<
- "Zoom: " << m_config.zMinZoom << " < |z|/cm < " << m_config.zMaxZoom << ", " <<
- m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << ", " <<
- m_config.phiMinZoom << " < phi/degrees < " << m_config.phiMaxZoom );
+ ATH_MSG_INFO( "Initializing " << name() << "\n" <<
+ "OutputFile: " << m_radMapsFileName << "\n" <<
+ "Material: " << m_config.material << "\n" <<
+ "2D Maps: " << m_config.nBinsz << " |z|-bins, " <<
+ m_config.nBinsr << " r-bins" << "\n" <<
+ "Zoom: " << m_config.zMinZoom << " < |z|/cm < " << m_config.zMaxZoom << ", " <<
+ m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << "\n" <<
+ "Full: " << m_config.zMinFull << " < |z|/cm < " << m_config.zMaxFull << ", " <<
+ m_config.rMinFull << " < r/cm < " << m_config.rMaxFull << "\n" <<
+ "Neutron Spectra: " << m_config.nBinslogEn << " log10E-bins" << ", " <<
+ 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" <<
+ "3D Maps: " << m_config.nBinsz3d << " |z|-bins, " <<
+ m_config.nBinsr3d << " r-bins, " <<
+ m_config.nBinsphi3d << " phi-bins" << "\n" <<
+ "Zoom: " << m_config.zMinZoom << " < |z|/cm < " << m_config.zMaxZoom << ", " <<
+ m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << ", " <<
+ m_config.phiMinZoom << " < phi/degrees < " << m_config.phiMaxZoom );
return StatusCode::SUCCESS;
}
@@ -98,6 +112,21 @@ namespace G4UA{
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);
+
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// 2d zoom
@@ -134,6 +163,20 @@ namespace G4UA{
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);
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// 2d zoom
@@ -245,6 +288,104 @@ namespace G4UA{
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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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("|z| [cm]");
+ 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}]");
+
TH2D * h_rz_vol = 0;
TH2D * h_rz_norm = 0;
TH2D * h_full_rz_vol = 0;
@@ -322,6 +463,30 @@ namespace G4UA{
// 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);
+ }
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// VOL
@@ -339,7 +504,14 @@ namespace G4UA{
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();
+
// 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++) {
@@ -369,6 +541,30 @@ namespace G4UA{
// 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);
+ }
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// VOL
@@ -386,6 +582,13 @@ namespace G4UA{
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();
// normalize to volume element per bin
for(int i=0;i<h_3d_tid->GetNbinsX();i++) { /* |z| */
--
GitLab