diff --git a/Common/SimConfig/CMakeLists.txt b/Common/SimConfig/CMakeLists.txt index f8e007209eacc..65d30935904ad 100644 --- a/Common/SimConfig/CMakeLists.txt +++ b/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 diff --git a/Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h b/Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h new file mode 100644 index 0000000000000..15d74ba27ab1b --- /dev/null +++ b/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 +#include +#include +#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 neutronG; + static std::unique_ptr protonG; + static std::unique_ptr pionG; +}; +#endif diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index aa8aa05263c0a..2a333a39e4242 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -49,8 +49,11 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { 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"); }; diff --git a/Common/SimConfig/src/FluenceWeightCalculator.cxx b/Common/SimConfig/src/FluenceWeightCalculator.cxx new file mode 100644 index 0000000000000..63828f71286e2 --- /dev/null +++ b/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 +#include +#include +#include + +std::unique_ptr FluenceWeightCalculator::neutronG; +std::unique_ptr FluenceWeightCalculator::protonG; +std::unique_ptr 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(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(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(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(); + neutronG->SetName("neutronDW"); + auto neuN = 0; + protonG = std::make_unique(); + protonG->SetName("protonDW"); + auto proN = 0; + pionG = std::make_unique(); + 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(); +} diff --git a/Detectors/gconfig/data/rd50_niel.root b/Detectors/gconfig/data/rd50_niel.root new file mode 100644 index 0000000000000..97174a714112a Binary files /dev/null and b/Detectors/gconfig/data/rd50_niel.root differ diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index c2b1fbd433e4b..16374cf9fd4a3 100644 --- a/Detectors/gconfig/g4Config.C +++ b/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 +#include #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;