diff --git a/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h b/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
index f5e266af515b43a3c06b823c37246da7b4079d45..9d3361bc1e7525f0fa28208311234278a4d9cb01 100644
--- a/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
+++ b/Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
@@ -61,6 +61,11 @@ namespace G4UA
int nBinslogEo = 45;
double logEMino = -2.; // min log10(E_kin/MeV)
double logEMaxo = 7.; // max log10(E_kin/MeV)
+
+ // time dependent TID maps
+ int nBinslogT = 20;
+ double logTMin = -9.; // log10(t_cut/s); first bin for t < 1 ns
+ double logTMax = 11.; // log10(t_cut/s); last bin for t < 3169 a
};
@@ -177,7 +182,14 @@ 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;
-
+
+ // time dependent maps
+
+ /// vector of time dependent TID in zoom 2d grid
+ std::vector<double> m_rz_tid_time;
+ /// vector of time dependent TID in full 2d grid
+ std::vector<double> m_full_rz_tid_time;
+
void merge(const Report& maps);
};
diff --git a/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py b/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
index 8636a2680dafb712ffc7df1d6381f272401cfa77..ed68e98eadd1028f912cb2f0175165abdfbebd0a 100644
--- a/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
+++ b/Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
@@ -10,6 +10,7 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
# radmaptool.NBinsZ = 240
# radmaptool.NBinsLogEn = 90
# radmaptool.NBinsLogEo = 45
+# radmaptool.NBinsLogTimeCut = 20
# radmaptool.NBinsR3D = 30
# radmaptool.NBinsZ3D = 60
# radmaptool.NBinsPhi3D = 32
@@ -27,6 +28,8 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
# radmaptool.LogEMaxn = 7.0 # in log10(E/MeV)
# radmaptool.LogEMino = -2.0 # in log10(E/MeV)
# 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)
#
simFlags.OptionalUserActionList.addAction('G4UA::RadiationMapsMakerTool')
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
index 79cd139b50f1d427495af877c347b617f6a83d8f..1451e193c46f940e30508a4da7c6e0ca448d4983 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
@@ -105,6 +105,15 @@ namespace G4UA{
m_3d_chad[i] += maps.m_3d_chad[i];
}
+ // all time 2d vectors have same size
+ for(unsigned int i=0;i<maps.m_rz_tid_time.size();i++) {
+ // zoom 2d
+ m_rz_tid_time [i] += maps.m_rz_tid_time [i];
+
+ // full 2d
+ m_full_rz_tid_time [i] += maps.m_full_rz_tid_time [i];
+ }
+
for(unsigned int i=0;i<maps.m_3d_vol.size();i++) {
// need volume fraction only if particular material is selected
m_3d_vol [i] += maps.m_3d_vol [i];
@@ -158,6 +167,9 @@ namespace G4UA{
m_maps.m_3d_neut.resize(0);
m_maps.m_3d_chad.resize(0);
+ m_maps.m_rz_tid_time .resize(0);
+ m_maps.m_full_rz_tid_time .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);
@@ -206,6 +218,9 @@ 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_tid_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,0.0);
+ m_maps.m_full_rz_tid_time .resize(m_config.nBinsz*m_config.nBinsr*m_config.nBinslogT,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);
@@ -338,6 +353,23 @@ namespace G4UA{
double rho = aStep->GetTrack()->GetMaterial()->GetDensity()/CLHEP::g*CLHEP::cm3;
+ // get time of step as average between pre- and post-step point
+ G4StepPoint* pre_step_point = aStep->GetPreStepPoint();
+ G4StepPoint* post_step_point = aStep->GetPostStepPoint();
+ G4ThreeVector startPoint = pre_step_point->GetPosition();
+ G4ThreeVector endPoint = post_step_point->GetPosition();
+ G4ThreeVector p = (startPoint + endPoint) * 0.5;
+
+ double timeOfFlight = (pre_step_point->GetGlobalTime() +
+ post_step_point->GetGlobalTime()) * 0.5;
+
+ // then compute time difference to a particle traveling with speed of light until this position
+ double deltatime = (timeOfFlight - p.mag()/CLHEP::c_light)/CLHEP::s;
+
+ // and use the log10 of that time difference in seconds to fill time-dependent histograms
+ // note that the upper bin boundary is the relevant cut. The first bin contains thus all times even shorter than the first lower bin boundary
+ double logt = (deltatime<0?m_config.logTMin-1:log10(deltatime));
+
bool goodMaterial(false);
if (m_config.material.empty()) {
@@ -421,6 +453,8 @@ namespace G4UA{
int vBinFullSpecn = -1;
int vBinZoomSpeco = -1;
int vBinFullSpeco = -1;
+ int vBinZoomTime = -1;
+ int vBinFullTime = -1;
// zoom 2d
if ( m_config.zMinZoom < abszorz &&
@@ -430,6 +464,10 @@ namespace G4UA{
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.logTMax > logt ){
+ int ilt = (logt < m_config.logTMin?0:(logt-m_config.logTMin)/(m_config.logTMax-m_config.logTMin)*m_config.nBinslogT);
+ vBinZoomTime = m_config.nBinsr*m_config.nBinslogT*iz+ir*m_config.nBinslogT+ilt;
+ }
if ( m_config.logEMinn < logEKin &&
m_config.logEMaxn > logEKin &&
(pdgid == 6 || pdgid == 7)) {
@@ -453,6 +491,10 @@ namespace G4UA{
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.logTMax > logt ){
+ int ilt = (logt < m_config.logTMin?0:(logt-m_config.logTMin)/(m_config.logTMax-m_config.logTMin)*m_config.nBinslogT);
+ vBinFullTime = m_config.nBinsr*m_config.nBinslogT*iz+ir*m_config.nBinslogT+ilt;
+ }
if ( m_config.logEMinn < logEKin &&
m_config.logEMaxn > logEKin &&
(pdgid == 6 || pdgid == 7)) {
@@ -506,6 +548,9 @@ namespace G4UA{
m_maps.m_rz_eion[vBinZoom] += dE_ION;
}
}
+ if ( goodMaterial && vBinZoomTime >=0 ) {
+ m_maps.m_rz_tid_time[vBinZoomTime] += dE_ION/rho;
+ }
if ( goodMaterial && vBinFull >=0 ) {
if ( pdgid == 999 ) {
m_maps.m_full_rz_tid [vBinFull] += dl;
@@ -516,6 +561,9 @@ namespace G4UA{
m_maps.m_full_rz_eion[vBinFull] += dE_ION;
}
}
+ if ( goodMaterial && vBinFullTime >=0 ) {
+ m_maps.m_full_rz_tid_time[vBinFullTime] += dE_ION/rho;
+ }
if ( goodMaterial && vBin3d >=0 ) {
if ( pdgid == 999 ) {
m_maps.m_3d_tid [vBin3d] += dl;
diff --git a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
index 9ca7095103fe9f2cbf6993b3e2877084e4b12b59..aff0fefb21baebb2d63ef4b2f303d7a1f46a29c1 100644
--- a/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
+++ b/Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
@@ -17,39 +17,44 @@ namespace G4UA{
/// Output Filename for the Radiation Maps
declareProperty("RadMapsFileName", m_radMapsFileName);
/// Name of the material to make radiation maps for (take all if empty)
- declareProperty("Material" , m_config.material);
+ declareProperty("Material" , m_config.material);
/// map granularities
/// number of bins in r and z for all 2D maps
- declareProperty("NBinsR" , m_config.nBinsr);
- declareProperty("NBinsZ" , m_config.nBinsz);
+ 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);
+ declareProperty("NBinsLogEn" , m_config.nBinslogEn);
/// number of bins in logE for energy spectra of other particles in 2D grids
- declareProperty("NBinsLogEo", m_config.nBinslogEo);
+ 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);
- declareProperty("NBinsPhi3D", m_config.nBinsphi3d);
+ declareProperty("NBinsR3D" , m_config.nBinsr3d);
+ declareProperty("NBinsZ3D" , m_config.nBinsz3d);
+ declareProperty("NBinsPhi3D" , m_config.nBinsphi3d);
+ /// number of bins in logTimeCut for time dependent TID 2D maps
+ declareProperty("NBinsLogTimeCut", m_config.nBinslogT);
/// map ranges
/// for Zoomed area in 2D and 3D
- declareProperty("RMinZoom" , m_config.rMinZoom);
- declareProperty("RMaxZoom" , m_config.rMaxZoom);
- declareProperty("ZMinZoom" , m_config.zMinZoom);
- declareProperty("ZMaxZoom" , m_config.zMaxZoom);
+ declareProperty("RMinZoom" , m_config.rMinZoom);
+ declareProperty("RMaxZoom" , m_config.rMaxZoom);
+ declareProperty("ZMinZoom" , m_config.zMinZoom);
+ declareProperty("ZMaxZoom" , m_config.zMaxZoom);
/// for Full detector in 2D
- declareProperty("RMinFull" , m_config.rMinFull);
- declareProperty("RMaxFull" , m_config.rMaxFull);
- declareProperty("ZMinFull" , m_config.zMinFull);
- declareProperty("ZMaxFull" , m_config.zMaxFull);
+ declareProperty("RMinFull" , m_config.rMinFull);
+ declareProperty("RMaxFull" , m_config.rMaxFull);
+ declareProperty("ZMinFull" , m_config.zMinFull);
+ declareProperty("ZMaxFull" , m_config.zMaxFull);
/// for Zoomed area in 3D
- declareProperty("PhiMinZoom", m_config.phiMinZoom);
- declareProperty("PhiMaxZoom", m_config.phiMaxZoom);
+ 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);
+ 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);
+ declareProperty("LogEMino" , m_config.logEMino);
+ declareProperty("LogEMaxo" , m_config.logEMaxo);
+ /// for logT in time-dependent TID 2D maps
+ declareProperty("LogTMin" , m_config.logTMin);
+ declareProperty("LogTMax" , m_config.logTMax);
}
//---------------------------------------------------------------------------
@@ -57,25 +62,27 @@ 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 << (m_config.zMinFull<0?" z-bins, ":" |z|-bins, ") <<
- m_config.nBinsr << " r-bins" << "\n" <<
- "Zoom: " << m_config.zMinZoom << (m_config.zMinFull<0?" < z/cm < ":" < |z|/cm < ") << m_config.zMaxZoom << ", " <<
- m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << "\n" <<
- "Full: " << m_config.zMinFull << (m_config.zMinFull<0?" < z/cm < ":" < |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 << (m_config.zMinFull<0?" z-bins, ":" |z|-bins, ") <<
- m_config.nBinsr3d << " r-bins, " <<
- m_config.nBinsphi3d << " phi-bins" << "\n" <<
- "Zoom: " << m_config.zMinZoom << (m_config.zMinFull<0?" < z/cm < ":" < |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 << (m_config.zMinFull<0?" z-bins, ":" |z|-bins, ") <<
+ m_config.nBinsr << " r-bins" << "\n" <<
+ "Zoom: " << m_config.zMinZoom << (m_config.zMinFull<0?" < z/cm < ":" < |z|/cm < ") << m_config.zMaxZoom << ", " <<
+ m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << "\n" <<
+ "Full: " << m_config.zMinFull << (m_config.zMinFull<0?" < z/cm < ":" < |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 << (m_config.zMinFull<0?" z-bins, ":" |z|-bins, ") <<
+ m_config.nBinsr3d << " r-bins, " <<
+ m_config.nBinsphi3d << " phi-bins" << "\n" <<
+ "Zoom: " << m_config.zMinZoom << (m_config.zMinFull<0?" < z/cm < ":" < |z|/cm < ") << m_config.zMaxZoom << ", " <<
+ m_config.rMinZoom << " < r/cm < " << m_config.rMaxZoom << ", " <<
+ m_config.phiMinZoom << " < phi/degrees < " << m_config.phiMaxZoom << "\n" <<
+ "Time TID Maps: " << m_config.nBinslogT << " Time-cut bins, " <<
+ m_config.logTMin << " < log10(t_cut/s) < " << m_config.logTMax );
return StatusCode::SUCCESS;
}
@@ -127,6 +134,9 @@ namespace G4UA{
maps.m_rz_rest_spec .resize(0);
maps.m_full_rz_rest_spec.resize(0);
+ maps.m_rz_tid_time .resize(0);
+ maps.m_full_rz_tid_time .resize(0);
+
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// 2d zoom
@@ -177,6 +187,10 @@ namespace G4UA{
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);
+
+ 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);
+
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// 2d zoom
@@ -388,6 +402,16 @@ namespace G4UA{
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
+ 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)");
+
TH2D * h_rz_vol = 0;
TH2D * h_rz_norm = 0;
TH2D * h_full_rz_vol = 0;
@@ -435,7 +459,6 @@ namespace G4UA{
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++) {
@@ -491,6 +514,13 @@ namespace G4UA{
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);
+ int vBinT = m_config.nBinsr*m_config.nBinslogT*i+j*m_config.nBinslogT+k;
+ val =maps.m_rz_tid_time[vBinT];
+ h_rz_tid_time->SetBinContent(kBin,val/vol);
+ }
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// VOL
@@ -515,6 +545,7 @@ namespace G4UA{
h_rz_pion_spec->Write();
h_rz_prot_spec->Write();
h_rz_rest_spec->Write();
+ h_rz_tid_time->Write();
// normalize to volume element per bin
for(int i=0;i<h_full_rz_tid->GetNbinsX();i++) {
@@ -571,6 +602,13 @@ namespace G4UA{
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);
+ int vBinT = m_config.nBinsr*m_config.nBinslogT*i+j*m_config.nBinslogT+k;
+ val =maps.m_full_rz_tid_time[vBinT];
+ h_full_rz_tid_time->SetBinContent(kBin,val/vol);
+ }
if (!m_config.material.empty()) {
// need volume fraction only if particular material is selected
// VOL
@@ -595,6 +633,7 @@ namespace G4UA{
h_full_rz_pion_spec->Write();
h_full_rz_prot_spec->Write();
h_full_rz_rest_spec->Write();
+ h_full_rz_tid_time->Write();
// normalize to volume element per bin
for(int i=0;i<h_3d_tid->GetNbinsX();i++) { /* |z| */