optional use of Geant4 fluence weighting added

Common/SimConfig/CMakeLists.txt

@@ -20,6 +20,7 @@ o2_add_library(SimConfig
                        src/MatMapParams.cxx
                        src/InteractionDiamondParam.cxx
                        src/GlobalProcessCutSimParam.cxx
+		       src/FluenceWeightCalculator.cxx
                PUBLIC_LINK_LIBRARIES O2::CommonUtils
                                      O2::DetectorsCommonDataFormats O2::SimulationDataFormat
                                      FairRoot::Base Boost::program_options)
@@ -35,7 +36,8 @@ o2_target_root_dictionary(SimConfig
                                   include/SimConfig/G4Params.h
                                   include/SimConfig/DetectorLists.h
                                   include/SimConfig/GlobalProcessCutSimParam.h
-                                  include/SimConfig/MatMapParams.h)
+                                  include/SimConfig/MatMapParams.h
+				  include/SimConfig/FluenceWeightCalculator.h)
 
 o2_add_test(Config
             SOURCES test/TestConfig.cxx

Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h

@@ -0,0 +1,35 @@
+// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef FluenceWeightCalculator_h
+#define FluenceWeightCalculator_h
+#include <vector>
+#include <string>
+#include <memory>
+#include "TGraph.h"
+//
+// Static container class for damage weight funnctions in form of TGraphs
+// The weights can be read from a csv file and stored in the graphs.
+//
+class FluenceWeightCalculator
+{
+ public:
+  FluenceWeightCalculator() = delete;
+  static void InitWeights(const std::string& filename);
+  static void InitWeightsFromCSV(const std::string& filename);
+  static double GetWeight(const int pdg, const double ekin);
+
+ private:
+  static std::unique_ptr<TGraph> neutronG;
+  static std::unique_ptr<TGraph> protonG;
+  static std::unique_ptr<TGraph> pionG;
+};
+#endif

Common/SimConfig/include/SimConfig/G4Params.h

@@ -49,8 +49,11 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
 
   EG4Nav navmode = EG4Nav::kTGeo; // geometry navigation mode (default TGeo)
 
+  std::string fluenceWeightFile = ""; // file containing the scoring weights (pdg, ekin, weight)
   std::string const& getPhysicsConfigString() const;
 
+  bool g4scoring = false;
+  bool g4fluenceweight = false;
   O2ParamDef(G4Params, "G4");
 };
 

Common/SimConfig/src/FluenceWeightCalculator.cxx

@@ -0,0 +1,138 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include "SimConfig/FluenceWeightCalculator.h"
+#include <TFile.h>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+
+std::unique_ptr<TGraph> FluenceWeightCalculator::neutronG;
+std::unique_ptr<TGraph> FluenceWeightCalculator::protonG;
+std::unique_ptr<TGraph> FluenceWeightCalculator::pionG;
+
+double FluenceWeightCalculator::GetWeight(const int pdg, const double kineticEnergy)
+{
+  //
+  // Obtain weight for given particle and kinetic energy
+  if (!neutronG || !protonG || !pionG) {
+    std::cerr << "FluenceWeightCalculator not initialized\n";
+    return 0.;
+  }
+  switch (std::abs(pdg)) {
+    case 2112: {
+      return neutronG->Eval(kineticEnergy, nullptr, "S");
+    }
+    case 2212: {
+      return ((kineticEnergy > 1e-3) ? protonG->Eval(kineticEnergy, nullptr, "S") : 0.);
+    }
+    case 211: {
+      return ((kineticEnergy > 10.) ? pionG->Eval(kineticEnergy, nullptr, "S") : 0.);
+    }
+    default:
+      return 0.0;
+  }
+}
+
+void FluenceWeightCalculator::InitWeights(const std::string& filename)
+{
+  //
+  // Read graphs from file
+  TFile inFile(filename.c_str(), "READ");
+  if (inFile.IsZombie()) {
+    std::cerr << "Cannot open " << filename << "\n";
+    return;
+  }
+  //
+  TGraph* tmp = nullptr;
+  inFile.GetObject("neutronDW", tmp);
+  neutronG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
+  if (!neutronG) {
+    std::cerr << "Missing graph neutronDW\n";
+    return;
+  }
+  neutronG->SetBit(TGraph::kIsSortedX);
+  inFile.GetObject("protonDW", tmp);
+  protonG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
+  if (!protonG) {
+    std::cerr << "Missing graph protonDW\n";
+    return;
+  }
+  protonG->SetBit(TGraph::kIsSortedX);
+  inFile.GetObject("pionDW", tmp);
+  pionG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
+  if (!pionG) {
+    std::cerr << "Missing graph pionDW\n";
+    return;
+  }
+  pionG->SetBit(TGraph::kIsSortedX);
+}
+
+void FluenceWeightCalculator::InitWeightsFromCSV(const std::string& filename)
+{
+  //
+  // read the NIEL weights from input file and store them as TGraphs
+  neutronG = std::make_unique<TGraph>();
+  neutronG->SetName("neutronDW");
+  auto neuN = 0;
+  protonG = std::make_unique<TGraph>();
+  protonG->SetName("protonDW");
+  auto proN = 0;
+  pionG = std::make_unique<TGraph>();
+  pionG->SetName("pionDW");
+  auto pioN = 0;
+
+  std::ifstream in(filename);
+  if (!in.is_open()) {
+    std::cerr << "Error: cannot open file with damage weights.\n";
+    return;
+  }
+  std::string line;
+  while (std::getline(in, line)) {
+    if (line.empty() || line[0] == '#') {
+      continue;
+    }
+    std::istringstream ss(line);
+    std::string particle, e_str, w_str;
+    if (!std::getline(ss, particle, ',')) {
+      continue;
+    }
+    if (!std::getline(ss, e_str, ',')) {
+      continue;
+    }
+    if (!std::getline(ss, w_str, ',')) {
+      continue;
+    }
+    auto e = std::stod(e_str);
+    auto w = std::stod(w_str);
+    auto pdg = std::stoi(particle);
+    switch (pdg) {
+      case 2112: {
+        neutronG->SetPoint(neuN++, e, w);
+        break;
+      }
+      case 2212: {
+        protonG->SetPoint(proN++, e, w);
+        break;
+      }
+      case 211: {
+        pionG->SetPoint(pioN++, e, w);
+        break;
+      }
+      default:;
+    }
+  }
+  auto fout = new TFile("rd50_niel.root", "recreate");
+  neutronG->Write();
+  protonG->Write();
+  pionG->Write();
+  fout->Close();
+}

Detectors/gconfig/g4Config.C

@@ -1,10 +1,13 @@
-/********************************************************************************
- *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
- *                                                                              *
- *              This software is distributed under the terms of the             *
- *         GNU Lesser General Public Licence version 3 (LGPL) version 3,        *
- *                  copied verbatim in the file "LICENSE"                       *
- ********************************************************************************/
+// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
 R__LOAD_LIBRARY(libG4ptl)
 R__LOAD_LIBRARY(libG4zlib)
 R__LOAD_LIBRARY(libG4expat)
@@ -50,12 +53,14 @@ R__LOAD_LIBRARY(libgeant4vmc)
 
 #if !defined(__CLING__) || defined(__ROOTCLING__)
 #include <iostream>
+#include <functional>
 #include "TGeant4.h"
 #include "TString.h"
 #include "FairRunSim.h"
 #include "TSystem.h"
 #include "TG4RunConfiguration.h"
 #include "SimConfig/G4Params.h"
+#include "SimConfig/FluenceWeightCalculator.h"
 #endif
 #include "commonConfig.C"
 
@@ -111,9 +116,15 @@ void Config()
 
   auto runConfiguration = new TG4RunConfiguration(geomNavStr, physicsSetup, "stepLimiter+specialCuts",
                                                   specialStacking, mtMode);
+  if (g4Params.g4scoring) {
+    runConfiguration->SetUseOfG4Scoring();
+    if (g4Params.g4fluenceweight) {
+      FluenceWeightCalculator::InitWeightsFromCSV(g4Params.fluenceWeightFile);
+      runConfiguration->SetScoreWeightCalculator(&FluenceWeightCalculator::GetWeight);
+    }
+  }
   /// avoid the use of G4BACKTRACE (it seems to inferfere with process logic in o2-sim)
   setenv("G4BACKTRACE", "none", 1);
-
   /// Create the G4 VMC
   TGeant4* geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration);
   std::cout << "Geant4 has been created." << std::endl;