Skip to content
Snippets Groups Projects
Commit 0363b74c authored by Adam Edward Barton's avatar Adam Edward Barton
Browse files

Merge branch 'master-RadMaps' into 'master' 

RadiationMapsMaker: manual port of latest 21.0 additions to master (Time-dependent TID-maps)

See merge request atlas/athena!24601
parents e83ae3c3 6428da66
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h
+ 13
− 1
View file @ 0363b74c
......@@ -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);
};
......
Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
+ 3
− 0
View file @ 0363b74c
......@@ -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')
Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
+ 48
− 0
View file @ 0363b74c
......@@ -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;
......
Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx
+ 81
− 42
View file @ 0363b74c
......@@ -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| */
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment