RadiationMapsMaker: Added option to generate maps for just one specific

                    material. The volume fraction of the target material
                    is measured on the fly in each bin.

Former-commit-id: 415cbe2203bd8fa5fe114bbcf5d05785663e14d7

Simulation/G4Utilities/G4UserActions/G4UserActions/RadiationMapsMaker.h

@@ -6,6 +6,7 @@
 #define G4UserActions_RadiationMapsMaker_H
 
 #include <vector>
+#include <string>
 #include "G4UserRunAction.hh"
 #include "G4UserSteppingAction.hh"
 
@@ -27,6 +28,8 @@ namespace G4UA
 	/// Radiation Estimate Web tool for the default values given here.
 	/// They can be configured to other values/ranges for other purposes.
 
+	std::string material = std::string("");
+
 	int nBinsr = 120;
 	int nBinsz = 240;
 	
@@ -63,6 +66,13 @@ namespace G4UA
         /// vector of >20 MeV hadron flux seen by thread in zoomed area
 	std::vector<double> m_rz_h20;
 
+	///  next two vectors are used only in case maps are needed for a particular material instead of all
+
+        /// vector to measure volume fraction of target material in zoomed area
+	std::vector<double> m_rz_vol;
+        /// vector to normalize the volume fraction in zoomed area
+	std::vector<double> m_rz_norm;
+
         /// vector of tid seen by thread in full area
 	std::vector<double> m_full_rz_tid;
         /// vector of ionizing energy density seen by thread in full area
@@ -72,6 +82,13 @@ namespace G4UA
         /// vector of >20 MeV hadron flux seen by thread in full area
 	std::vector<double> m_full_rz_h20;
 
+	///  next two vectors are used only in case maps are needed for a particular material instead of all
+
+        /// vector to measure volume fraction of target material in full area
+	std::vector<double> m_full_rz_vol;
+        /// vector to normalize the volume fraction in full area
+	std::vector<double> m_full_rz_norm;
+
         /// vector of tid seen by thread in 3d
 	std::vector<double> m_3d_tid;
         /// vector of ionizing energy density seen by thread in 3d
@@ -81,6 +98,13 @@ namespace G4UA
         /// vector of >20 MeV hadron flux seen by thread in 3d
 	std::vector<double> m_3d_h20;
 
+	///  next two vectors are used only in case maps are needed for a particular material instead of all
+
+        /// vector to measure volume fraction of target material in 3d
+	std::vector<double> m_3d_vol;
+        /// vector to normalize the volume fraction in 3d
+	std::vector<double> m_3d_norm;
+
 	void merge(const Report& maps);
       };
 

Simulation/G4Utilities/G4UserActions/share/jobOptions.RadiationMapsMaker.py

@@ -5,6 +5,7 @@ from G4UserActions.G4UserActionsConf import G4UA__RadiationMapsMakerTool
 #
 # radmaptool = G4UA__RadiationMapsMakerTool("G4UA::RadiationMapsMakerTool")
 # radmaptool.RadMapsFileName = "RadMaps.root" 
+# radmaptool.Material        = "" # if left empty all materials are used (default)
 # radmaptool.NBinsR          =   120
 # radmaptool.NBinsZ          =   240
 # radmaptool.NBinsR3D        =    30

Simulation/G4Utilities/G4UserActions/src/RadiationMapsMaker.cxx

@@ -81,6 +81,14 @@ namespace G4UA{
       m_full_rz_h20 [i] += maps.m_full_rz_h20 [i];
     }
 
+    for(unsigned int i=0;i<maps.m_rz_vol.size();i++) {
+      // need volume fraction only if particular material is selected
+      m_rz_vol [i] += maps.m_rz_vol [i];
+      m_rz_norm[i] += maps.m_rz_norm[i];
+      m_full_rz_vol [i] += maps.m_full_rz_vol [i];
+      m_full_rz_norm[i] += maps.m_full_rz_norm[i];
+    }
+
     // all zoom 3d vectors have same size 
     for(unsigned int i=0;i<maps.m_3d_tid.size();i++) {
       // zoom 3d
@@ -89,6 +97,12 @@ namespace G4UA{
       m_3d_niel[i] += maps.m_3d_niel[i];
       m_3d_h20 [i] += maps.m_3d_h20 [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];
+      m_3d_norm[i] += maps.m_3d_norm[i];
+    }
   }
 
   void RadiationMapsMaker::BeginOfRunAction(const G4Run*){
@@ -110,6 +124,16 @@ namespace G4UA{
     m_maps.m_3d_niel.resize(0);
     m_maps.m_3d_h20 .resize(0);
 
+    if (!m_config.material.empty()) {
+      // need volume fraction only if particular material is selected
+      m_maps.m_rz_vol .resize(0);
+      m_maps.m_rz_norm.resize(0);
+      m_maps.m_full_rz_vol .resize(0);
+      m_maps.m_full_rz_norm.resize(0);
+      m_maps.m_3d_vol .resize(0);
+      m_maps.m_3d_norm.resize(0);
+    }
+
     // then resize to proper size and initialize with 0's 
 
     m_maps.m_rz_tid .resize(m_config.nBinsz*m_config.nBinsr,0.0);
@@ -127,6 +151,19 @@ namespace G4UA{
     m_maps.m_3d_niel.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
     m_maps.m_3d_h20 .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
 
+    if (!m_config.material.empty()) {
+      // need volume fraction only if particular material is selected
+      // 2d zoom
+      m_maps.m_rz_vol .resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      m_maps.m_rz_norm.resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      // 2d full
+      m_maps.m_full_rz_vol .resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      m_maps.m_full_rz_norm.resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      // 3d
+      m_maps.m_3d_vol .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
+      m_maps.m_3d_norm.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
+    }
+
     /// data files for NIEL damage factors in Si
     ///
     /// The data resides in the share directory of the package and should not be
@@ -227,6 +264,19 @@ namespace G4UA{
     
       double rho = aStep->GetTrack()->GetMaterial()->GetDensity()/CLHEP::g*CLHEP::cm3; 
 
+      bool goodMaterial(false);
+
+      if (m_config.material.empty()) {
+	// if no material is specified maps are made for all materials
+	goodMaterial = true;
+      } 
+      else {
+	// if a material is specified maps are made just for that one.
+	// Note that still all steps need to be treated to calculate the
+	// volume fraction of the target material in each bin!
+	goodMaterial = (m_config.material.compare(aStep->GetTrack()->GetMaterial()->GetName()) == 0);
+      } 
+
       // fluence dN/da = dl/dV 
       double x0 = aStep->GetPreStepPoint()->GetPosition().x()*0.1;
       double x1 = aStep->GetPostStepPoint()->GetPosition().x()*0.1;
@@ -342,20 +392,20 @@ namespace G4UA{
 	    }
 	  }
 	  
-	  if ( vBinZoom >=0 ) {
+	  if ( goodMaterial && vBinZoom >=0 ) {
 	    m_maps.m_rz_tid [vBinZoom] += dE_ION/rho;
 	    m_maps.m_rz_eion[vBinZoom] += dE_ION;
 	  }
-	  if ( vBinFull >=0 ) {
+	  if ( goodMaterial && vBinFull >=0 ) {
 	    m_maps.m_full_rz_tid [vBinFull] += dE_ION/rho;
 	    m_maps.m_full_rz_eion[vBinFull] += dE_ION;
 	  }
-	  if ( vBin3d >=0 ) {
+	  if ( goodMaterial && vBin3d >=0 ) {
 	    m_maps.m_3d_tid [vBin3d] += dE_ION/rho;
 	    m_maps.m_3d_eion[vBin3d] += dE_ION;
 	  }
 	    
-	  if ( pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 ) {
+	  if ( goodMaterial && (pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 )) {
 	    
 	    if ( weight > 0 ) {
 	      if ( vBinZoom >=0 ) {
@@ -380,6 +430,40 @@ namespace G4UA{
 	      }
 	    }
 	  }
+	  if (!m_config.material.empty()) {
+	    // need volume fraction only if particular material is selected
+	    if ( eKin > 1 && (pdgid == 6 || pdgid == 7)) {
+	      // count all neutron > 1 MeV track length 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
+	      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_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;
+		}
+	      }
+	    }
+	  }
 	}
       }
     }

Simulation/G4Utilities/G4UserActions/src/RadiationMapsMakerTool.cxx

@@ -19,6 +19,8 @@ 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);
     /// map granularities 
     /// number of bins in r and z for all 2D maps
     declareProperty("NBinsR"    , m_config.nBinsr);
@@ -50,6 +52,7 @@ namespace G4UA{
   {
     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   << ", " << 
@@ -90,6 +93,19 @@ namespace G4UA{
     m_report.m_3d_niel.resize(0);
     m_report.m_3d_h20 .resize(0);
 
+    if (!m_config.material.empty()) {
+      // need volume fraction only if particular material is selected
+      // 2d zoom
+      m_report.m_rz_vol .resize(0);
+      m_report.m_rz_norm.resize(0);
+      // 2d full
+      m_report.m_full_rz_vol.resize(0);
+      m_report.m_full_rz_norm.resize(0);
+      // 3d
+      m_report.m_3d_vol .resize(0);
+      m_report.m_3d_norm.resize(0);
+    }
+
     // then resize to proper size and initialize with 0's 
 
     m_report.m_rz_tid .resize(m_config.nBinsz*m_config.nBinsr,0.0);
@@ -107,6 +123,19 @@ namespace G4UA{
     m_report.m_3d_niel.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
     m_report.m_3d_h20 .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
 
+    if (!m_config.material.empty()) {
+      // need volume fraction only if particular material is selected
+      // 2d zoom
+      m_report.m_rz_vol .resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      m_report.m_rz_norm.resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      // 2d full
+      m_report.m_full_rz_vol .resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      m_report.m_full_rz_norm.resize(m_config.nBinsz*m_config.nBinsr,0.0);
+      // 3d
+      m_report.m_3d_vol .resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
+      m_report.m_3d_norm.resize(m_config.nBinsz3d*m_config.nBinsr3d*m_config.nBinsphi3d,0.0);
+    }
+
     // merge radiation map vectors from threads
     mergeReports();
 
@@ -177,6 +206,53 @@ namespace G4UA{
     h_3d_niel ->SetTitle("NIEL [n_{eq}/cm^{2}]");
     h_3d_h20  ->SetTitle("SEE [h_{>20 MeV}/cm^{2}]");
 
+    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("|z| [cm]");
+      h_rz_norm ->SetXTitle("|z| [cm]");
+      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("|z| [cm]");
+      h_full_rz_norm ->SetXTitle("|z| [cm]");
+      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("|z| [cm]");
+      h_3d_norm ->SetXTitle("|z| [cm]");
+      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++) { 
@@ -201,6 +277,15 @@ namespace G4UA{
 	// SEE
 	val =m_report.m_rz_h20[vBin];
 	h_rz_h20->SetBinContent(iBin,val/vol);
+	if (!m_config.material.empty()) {
+	  // need volume fraction only if particular material is selected
+	  // VOL
+	  val =m_report.m_rz_vol[vBin];
+	  h_rz_vol->SetBinContent(iBin,val/vol);
+	  // NORM
+	  val =m_report.m_rz_norm[vBin];
+	  h_rz_norm->SetBinContent(iBin,val/vol);
+	}
       }
     }
     h_rz_tid->Write();
@@ -231,6 +316,15 @@ namespace G4UA{
 	// SEE
 	val =m_report.m_full_rz_h20[vBin];
 	h_full_rz_h20->SetBinContent(iBin,val/vol);
+	if (!m_config.material.empty()) {
+	  // need volume fraction only if particular material is selected
+	  // VOL
+	  val =m_report.m_full_rz_vol[vBin];
+	  h_full_rz_vol->SetBinContent(iBin,val/vol);
+	  // NORM
+	  val =m_report.m_full_rz_norm[vBin];
+	  h_full_rz_norm->SetBinContent(iBin,val/vol);
+	}
       }
     }
     h_full_rz_tid->Write();
@@ -269,6 +363,15 @@ namespace G4UA{
 	  // SEE
 	  val =m_report.m_3d_h20[vBin];
 	  h_3d_h20->SetBinContent(iBin,val/vol);
+	  if (!m_config.material.empty()) {
+	    // need volume fraction only if particular material is selected
+	    // VOL
+	    val =m_report.m_3d_vol[vBin];
+	    h_3d_vol->SetBinContent(iBin,val/vol);
+	    // NORM
+	    val =m_report.m_3d_norm[vBin];
+	    h_3d_norm->SetBinContent(iBin,val/vol);
+	  }
 	}
       }
     }
@@ -277,6 +380,16 @@ namespace G4UA{
     h_3d_niel->Write();
     h_3d_h20->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();
+    }
+
     f->Close();
 
     return StatusCode::SUCCESS;