Skip to content
Snippets Groups Projects
Commit aa717740 authored by Frank Winklmeier's avatar Frank Winklmeier
Browse files

Merge branch 'master-RadMaps-Direction' into 'master' 

RadiationMapsMaker: Direction info for spectra and memory footprint improvement in finalize

See merge request atlas/athena!34080
parents 131666c6 621bea71
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 566 additions and 536 deletions
Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py
+ 4
− 0
View file @ aa717740
......@@ -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
#
......
Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx
+ 85
− 1
View file @ aa717740
......@@ -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()) {
......
Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.h
+ 25
− 1
View file @ aa717740
......@@ -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
......
Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.h
+ 1
− 1
View file @ aa717740
/*
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
......
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