GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CascadeReweight.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Julia Tena-Vidal <j.tena-vidal \at liverpool.ac.uk>
7  University of Liverpool
8 
9 */
10 //____________________________________________________________________________
11 
12 #include <cstdlib>
13 
27 
31 
32 using namespace genie;
33 using namespace genie::utils;
34 using namespace genie::constants;
35 
36 //___________________________________________________________________________
38  : EventRecordVisitorI("genie::CascadeReweight") {}
39 //___________________________________________________________________________
41  : EventRecordVisitorI("genie::CascadeReweight", config) {}
42 //___________________________________________________________________________
44 //___________________________________________________________________________
46  if (!evrec) {
47  LOG("CascadeReweight", pERROR) << "** Null input!";
48  return;
49  }
50  // Get Associated weight
51  double weight = GetEventWeight(*evrec);
52  // Set weight
53  evrec->SetWeight(weight);
54 
55  return;
56 }
57 //___________________________________________________________________________
58 double CascadeReweight::GetEventWeight(const GHepRecord &event) const {
59 
60  GHepParticle *p = 0;
61  TIter event_iter(&event);
62  double total_weight = 1.;
63  while ((p = dynamic_cast<GHepParticle *>(event_iter.Next()))) {
64  // Look only at stable particles in the nucleus:
65  if ( p->Status() != kIStHadronInTheNucleus ) continue ;
66  // Get particle fate
67  auto fate_rescatter = p->RescatterCode();
68  // Only look at particles that had FSI
69  if (fate_rescatter < 0 || fate_rescatter == kIHNFtUndefined)
70  continue;
71  INukeFateHN_t fate = (INukeFateHN_t)fate_rescatter;
72 
73  // Read map weight:
74  const auto map_it = fFateWeightsMap.find(fate);
75  // Get weight given a pdg code.
76  if (map_it != fFateWeightsMap.end()) {
77  int pdg_target = p->Pdg();
78  const auto weight_it = (map_it->second).find(pdg_target);
79  if (weight_it != (map_it->second).end()) {
80  total_weight *= weight_it->second;
81  continue;
82  }
83  }
84  // If fate is not in the pdg map, use default values:
85  const auto def_it = fDefaultMap.find(fate);
86  if (def_it != fDefaultMap.end()) {
87  total_weight *= def_it->second;
88  }
89  } // end loop over particles
90 
91  return total_weight;
92 }
93 //___________________________________________________________________________
94 void CascadeReweight::Configure(const Registry &config) {
95  Algorithm::Configure(config);
96  this->LoadConfig();
97 }
98 //___________________________________________________________________________
99 void CascadeReweight::Configure(string param_set) {
100  Algorithm::Configure(param_set);
101  this->LoadConfig();
102 }
103 //___________________________________________________________________________
105  bool good_config = true;
106 
107  // Clean maps
108  fDefaultMap.clear();
109  fFateWeightsMap.clear();
110 
111  // Create vector with list of possible keys (follows the order of the fates
112  // enumeration)
113  std::map<INukeFateHN_t, string> EINukeFate_map_keys = GetEINukeFateKeysMap();
114 
115  for (map<INukeFateHN_t, string>::iterator it_keys =
116  EINukeFate_map_keys.begin();
117  it_keys != EINukeFate_map_keys.end(); it_keys++) {
118  // Find fate specifications
119  std::string to_find_def =
120  "CascadeReweight-Default-Weight-" + (it_keys->second);
121 
122  auto kdef_list = GetConfig().FindKeys(to_find_def.c_str());
123  for (auto kiter = kdef_list.begin(); kiter != kdef_list.end(); ++kiter) {
124  const RgKey &key = *kiter;
125  double weight;
126  GetParam(key, weight);
127  // Add check weight > 0
128  if (weight < 0) {
129  LOG("CascadeReweight", pERROR)
130  << "The weight assigned to " << to_find_def << " is not positive";
131  good_config = false;
132  continue;
133  }
134  fDefaultMap[it_keys->first] = weight;
135  }
136 
137  // Find Pdg specifications
138  std::string to_find_pdg =
139  "CascadeReweight-Weight-" + (it_keys->second) + "@Pdg=";
140  auto kpdg_list = GetConfig().FindKeys(to_find_pdg.c_str());
141  std::map<int, double> WeightMap; // define map that stores <pdg, weight>
142  for (auto kiter = kpdg_list.begin(); kiter != kpdg_list.end(); ++kiter) {
143  const RgKey &key = *kiter;
144  vector<string> kv = genie::utils::str::Split(key, "=");
145  assert(kv.size() == 2);
146  int pdg_target = stoi(kv[1]);
147  if (!PDGLibrary::Instance()->Find(pdg_target)) {
148  LOG("CascadeReweight", pERROR)
149  << "The target Pdg code associated to " << to_find_pdg
150  << " is not valid : " << pdg_target;
151  good_config = false;
152  continue;
153  }
154  double weight;
155  GetParam(key, weight);
156  // Add check weight > 0
157  if (weight < 0) {
158  LOG("CascadeReweight", pERROR)
159  << "The weight assigned to " << to_find_pdg << " is not positive";
160  good_config = false;
161  continue;
162  }
163  // Add pdg and weight in map
164  WeightMap.insert(std::pair<int, double>(pdg_target, weight));
165  }
166  // store information in class member
167  fFateWeightsMap[it_keys->first] = std::move(WeightMap);
168  }
169 
170  if (!good_config) {
171  LOG("CascadeReweight", pFATAL) << "Configuration has failed.";
172  exit(78);
173  }
174 }
175 //___________________________________________________________________________
int RescatterCode(void) const
Definition: GHepParticle.h:65
virtual void SetWeight(double wght)
Definition: GHepRecord.h:130
void ProcessEventRecord(GHepRecord *event_rec) const
#define pERROR
Definition: Messenger.h:59
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
#define pFATAL
Definition: Messenger.h:56
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
int Pdg(void) const
Definition: GHepParticle.h:63
RgKeyList FindKeys(RgKey key_part) const
create list with all keys containing &#39;key_part&#39;
Definition: Registry.cxx:840
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
enum genie::EINukeFateHN_t INukeFateHN_t
std::map< INukeFateHN_t, map< int, double > > fFateWeightsMap
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void Configure(const Registry &config)
static const std::map< INukeFateHN_t, string > & GetEINukeFateKeysMap(void)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
string RgKey
double GetEventWeight(const GHepRecord &ev) const
get weight from fate and configuration
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
void LoadConfig(void)
read configuration from xml file
std::map< INukeFateHN_t, double > fDefaultMap
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39