Skip to content

Commit d5972d7

Browse files
amorschamorsch
committed
optional use of Geant4 fluence weighting added
1 parent 3ee5c9f commit d5972d7

6 files changed

Lines changed: 164 additions & 1 deletion

File tree

‎Common/SimConfig/CMakeLists.txt‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ o2_add_library(SimConfig
2020
src/MatMapParams.cxx
2121
src/InteractionDiamondParam.cxx
2222
src/GlobalProcessCutSimParam.cxx
23+
src/FluenceWeightCalculator.cxx
2324
PUBLIC_LINK_LIBRARIES O2::CommonUtils
2425
O2::DetectorsCommonDataFormats O2::SimulationDataFormat
2526
FairRoot::Base Boost::program_options)
@@ -35,7 +36,8 @@ o2_target_root_dictionary(SimConfig
3536
include/SimConfig/G4Params.h
3637
include/SimConfig/DetectorLists.h
3738
include/SimConfig/GlobalProcessCutSimParam.h
38-
include/SimConfig/MatMapParams.h)
39+
include/SimConfig/MatMapParams.h
40+
include/SimConfig/FluenceWeightCalculator.h)
3941

4042
o2_add_test(Config
4143
SOURCES test/TestConfig.cxx

‎Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
2+
#ifndef FluenceWeightCalculator_h
3+
#define FluenceWeightCalculator_h
4+
#include <vector>
5+
#include <string>
6+
#include <memory>
7+
#include "TGraph.h"
8+
//
9+
// Static container class for damage weight funnctions in form of TGraphs
10+
// The weights can be read from a csv file and stored in the graphs.
11+
//
12+
class FluenceWeightCalculator
13+
{
14+
public:
15+
FluenceWeightCalculator() = delete;
16+
static void InitWeights(const std::string& filename);
17+
static void InitWeightsFromCSV(const std::string& filename);
18+
static double GetWeight(const int pdg, const double ekin);
19+
20+
private:
21+
static std::unique_ptr<TGraph> neutronG;
22+
static std::unique_ptr<TGraph> protonG;
23+
static std::unique_ptr<TGraph> pionG;
24+
};
25+
#endif

‎Common/SimConfig/include/SimConfig/G4Params.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,11 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
4949

5050
EG4Nav navmode = EG4Nav::kTGeo; // geometry navigation mode (default TGeo)
5151

52+
std::string fluenceWeightFile = ""; // file containing the scoring weights (pdg, ekin, weight)
5253
std::string const& getPhysicsConfigString() const;
5354

55+
bool g4scoring = false;
56+
bool g4fluenceweight = false;
5457
O2ParamDef(G4Params, "G4");
5558
};
5659

‎Common/SimConfig/src/FluenceWeightCalculator.cxx‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
#include "SimConfig/FluenceWeightCalculator.h"
2+
#include <TFile.h>
3+
#include <fstream>
4+
#include <sstream>
5+
#include <iostream>
6+
7+
std::unique_ptr<TGraph> FluenceWeightCalculator::neutronG;
8+
std::unique_ptr<TGraph> FluenceWeightCalculator::protonG;
9+
std::unique_ptr<TGraph> FluenceWeightCalculator::pionG;
10+
11+
double FluenceWeightCalculator::GetWeight(const int pdg, const double kineticEnergy)
12+
{
13+
//
14+
// Obtain weight for given particle and kinetic energy
15+
if (!neutronG || !protonG || !pionG) {
16+
std::cerr << "FluenceWeightCalculator not initialized\n";
17+
return 0.;
18+
}
19+
switch (std::abs(pdg)) {
20+
case 2112: {
21+
return neutronG->Eval(kineticEnergy, 0, "S");
22+
}
23+
case 2212: {
24+
return ((kineticEnergy > 1e-3) ? protonG->Eval(kineticEnergy, 0, "S") : 0.);
25+
}
26+
case 211: {
27+
return ((kineticEnergy > 10.) ? pionG->Eval(kineticEnergy, 0, "S") : 0.);
28+
}
29+
default:
30+
return 0.0;
31+
}
32+
}
33+
34+
void FluenceWeightCalculator::InitWeights(const std::string& filename)
35+
{
36+
//
37+
// Read graphs from file
38+
TFile inFile(filename.c_str(), "READ");
39+
if (inFile.IsZombie()) {
40+
std::cerr << "Cannot open " << filename << "\n";
41+
return;
42+
}
43+
//
44+
TGraph* tmp = nullptr;
45+
inFile.GetObject("neutronDW", tmp);
46+
neutronG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
47+
if (!neutronG) {
48+
std::cerr << "Missing graph neutronDW\n";
49+
return;
50+
}
51+
neutronG->SetBit(TGraph::kIsSortedX);
52+
inFile.GetObject("protonDW", tmp);
53+
protonG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
54+
if (!protonG) {
55+
std::cerr << "Missing graph protonDW\n";
56+
return;
57+
}
58+
protonG->SetBit(TGraph::kIsSortedX);
59+
inFile.GetObject("pionDW", tmp);
60+
pionG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
61+
if (!pionG) {
62+
std::cerr << "Missing graph pionDW\n";
63+
return;
64+
}
65+
pionG->SetBit(TGraph::kIsSortedX);
66+
}
67+
68+
void FluenceWeightCalculator::InitWeightsFromCSV(const std::string& filename)
69+
{
70+
//
71+
// read the NIEL weights from input file and store them as TGraphs
72+
neutronG = std::make_unique<TGraph>();
73+
neutronG->SetName("neutronDW");
74+
auto neuN = 0;
75+
protonG = std::make_unique<TGraph>();
76+
protonG->SetName("protonDW");
77+
auto proN = 0;
78+
pionG = std::make_unique<TGraph>();
79+
pionG->SetName("pionDW");
80+
auto pioN = 0;
81+
82+
std::ifstream in(filename);
83+
if (!in.is_open()) {
84+
std::cerr << "Error: cannot open file with damage weights.\n";
85+
return;
86+
}
87+
std::string line;
88+
while (std::getline(in, line)) {
89+
if (line.empty() || line[0] == '#')
90+
continue;
91+
std::istringstream ss(line);
92+
std::string particle, e_str, w_str;
93+
if (!std::getline(ss, particle, ','))
94+
continue;
95+
if (!std::getline(ss, e_str, ','))
96+
continue;
97+
if (!std::getline(ss, w_str, ','))
98+
continue;
99+
auto e = std::stod(e_str);
100+
auto w = std::stod(w_str);
101+
auto pdg = std::stoi(particle);
102+
switch (pdg) {
103+
case 2112: {
104+
neutronG->SetPoint(neuN++, e, w);
105+
break;
106+
}
107+
case 2212: {
108+
protonG->SetPoint(proN++, e, w);
109+
break;
110+
}
111+
case 211: {
112+
pionG->SetPoint(pioN++, e, w);
113+
break;
114+
}
115+
default:;
116+
}
117+
}
118+
auto fout = new TFile("rd50_niel.root", "recreate");
119+
neutronG->Write();
120+
protonG->Write();
121+
pionG->Write();
122+
fout->Close();
123+
}

‎Detectors/gconfig/data/rd50_niel.root‎

23.2 KB
Binary file not shown.

‎Detectors/gconfig/g4Config.C‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,12 +50,14 @@ R__LOAD_LIBRARY(libgeant4vmc)
5050

5151
#if !defined(__CLING__) || defined(__ROOTCLING__)
5252
#include <iostream>
53+
#include <functional>
5354
#include "TGeant4.h"
5455
#include "TString.h"
5556
#include "FairRunSim.h"
5657
#include "TSystem.h"
5758
#include "TG4RunConfiguration.h"
5859
#include "SimConfig/G4Params.h"
60+
#include "SimConfig/FluenceWeightCalculator.h"
5961
#endif
6062
#include "commonConfig.C"
6163

@@ -111,8 +113,16 @@ void Config()
111113

112114
auto runConfiguration = new TG4RunConfiguration(geomNavStr, physicsSetup, "stepLimiter+specialCuts",
113115
specialStacking, mtMode);
116+
if (g4Params.g4scoring) {
117+
runConfiguration->SetUseOfG4Scoring();
118+
if (g4Params.g4fluenceweight) {
119+
FluenceWeightCalculator::InitWeightsFromCSV(g4Params.fluenceWeightFile);
120+
runConfiguration->SetScoreWeightCalculator(&FluenceWeightCalculator::GetWeight);
121+
}
122+
}
114123
/// avoid the use of G4BACKTRACE (it seems to inferfere with process logic in o2-sim)
115124
setenv("G4BACKTRACE", "none", 1);
125+
116126

117127
/// Create the G4 VMC
118128
TGeant4* geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration);

0 commit comments

Comments
 (0)