gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
collider_event_weights.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
21 
23 
24 #define COLLIDERBIT_DEBUG
25 #define DEBUG_PREFIX "DEBUG: OMP thread " << omp_get_thread_num() << ": "
26 
27 namespace Gambit
28 {
29 
30  namespace ColliderBit
31  {
32 
34  void _setEventWeight_unity(HEPUtils::Event& event, const BaseCollider*) // <-- Ignoring second argument
35  {
36  event.set_weight(1.0);
37  event.set_weight_err(0.0);
38  }
39 
43  {
44  using namespace Pipes::setEventWeight_unity;
45  result = std::bind(_setEventWeight_unity, std::placeholders::_1, std::placeholders::_2);
46  }
47 
48 
49 
51  void _setEventWeight_fromCrossSection(HEPUtils::Event& event, const BaseCollider* HardScatteringSim_ptr, const map_int_process_xsec& ProcessCrossSectionsMap, const int use_trust_level)
52  {
53  // Initialize weight
54  double weight = 1.0;
55  double weight_err = 0.0;
56 
57  // Get process code from the generator
58  int process_code = HardScatteringSim_ptr->process_code();
59 
60  #ifdef COLLIDERBIT_DEBUG
61  cerr << DEBUG_PREFIX << "Current process_code: " << process_code << endl;
62  #endif
63 
64  // Get the process_xsec_container instance that holds the externally provided cross-section for this process
65  process_xsec_container xs = ProcessCrossSectionsMap.at(process_code);
66 
67  // Get the generator cross-section for this process
68  double process_xsec_generator = HardScatteringSim_ptr->xsec_fb(process_code);
69  double process_xsec_err_generator_sq = pow(HardScatteringSim_ptr->xsecErr_fb(process_code), 2);
70 
71  #ifdef COLLIDERBIT_DEBUG
72  cerr << DEBUG_PREFIX << "- process_code: " << process_code << ", xsec_fb: " << HardScatteringSim_ptr->xsec_fb(process_code)
73  << ", xsec_err_fb: " << HardScatteringSim_ptr->xsecErr_fb(process_code) << endl;
74  #endif
75 
76  // Determine what to do based on the trust_level of the externally provided cross-section:
77  if (xs.trust_level() >= use_trust_level)
78  {
79  // Add the generator cross-sections for other process codes which also
80  // contribute to the externaly provided cross-section
81  for (int other_process_code : xs.processes_sharing_xsec())
82  {
83  process_xsec_generator += HardScatteringSim_ptr->xsec_fb(other_process_code);
84  process_xsec_err_generator_sq += pow(HardScatteringSim_ptr->xsecErr_fb(other_process_code), 2);
85  #ifdef COLLIDERBIT_DEBUG
86  cerr << DEBUG_PREFIX << "- process_code: " << other_process_code << ", xsec_fb: " << HardScatteringSim_ptr->xsec_fb(other_process_code)
87  << ", xsec_err_fb: " << HardScatteringSim_ptr->xsecErr_fb(other_process_code) << endl;
88  #endif
89  }
90  double process_xsec_err_generator = sqrt(process_xsec_err_generator_sq);
91 
92  // Event weight = [external cross-section] / [sum of contributing generator cross-sections]
93  if (process_xsec_generator > 0.0)
94  {
95  weight = xs.xsec() / process_xsec_generator;
96  weight_err = sqrt( pow(xs.xsec_err() / process_xsec_generator, 2)
97  + pow(xs.xsec() * process_xsec_err_generator / pow(process_xsec_generator, 2), 2) );
98  }
99  else
100  {
101  std::stringstream errmsg_ss;
102  errmsg_ss << "Generated an event of process " << process_code << " for which the generator itself predicts a cross-section <= 0. Not sure what to do with that...";
103  ColliderBit_error().raise(LOCAL_INFO, errmsg_ss.str());
104  }
105  }
106  else
107  {
108  // Too low trust_level. Will fall back to use the generator cross-section
109  #ifdef COLLIDERBIT_DEBUG
110  cerr << DEBUG_PREFIX << "process_xsec trust_level too low (" << xs.trust_level() << "). Setting event weight to 1.0." << endl;
111  #endif
112  weight = 1.0;
113  weight_err = 0.0;
114  }
115 
116  #ifdef COLLIDERBIT_DEBUG
117  cerr << DEBUG_PREFIX << "total process_xsec: " << xs.xsec() << ", process_xsec_MC: " << process_xsec_generator << ", weight: " << weight << ", weight_err: " << weight_err << ", trust_level: " << xs.trust_level() << endl;
118  #endif
119 
120  event.set_weight(weight);
121  event.set_weight_err(weight_err);
122  }
123 
127  {
129 
130  const static int use_trust_level = runOptions->getValueOrDef<int>(1, "use_cross_section_trust_level");
131 
132  if(*Loop::iteration < 0) return;
133 
134  result = std::bind(_setEventWeight_fromCrossSection,
135  std::placeholders::_1,
136  std::placeholders::_2,
137  *Dep::ProcessCrossSectionsMap,
138  use_trust_level);
139  }
140 
141 
142  }
143 }
144 
145 
#define DEBUG_PREFIX
double xsec_err() const
Return the cross-section error (in fb).
Definition: xsec.cpp:60
int trust_level() const
Get the trust level.
Definition: xsec.cpp:148
virtual int process_code() const =0
Report an integer process code for the last generated event.
#define LOCAL_INFO
Definition: local_info.hpp:34
const std::vector< int > & processes_sharing_xsec() const
Return the list of process codes that share this cross-section (This is due to the many-to-many mappi...
Definition: xsec.cpp:350
Declarations common to all ColliderBit event loop functions.
void setEventWeight_fromCrossSection(EventWeighterFunctionType &result)
Module function providing an instance of EventWeighterFunctionType pointing to _setEventWeight_fromCr...
A class for holding the cross-section of a single Pythia process (identified by the Pythia process co...
Definition: xsec.hpp:148
void _setEventWeight_fromCrossSection(HEPUtils::Event &event, const BaseCollider *HardScatteringSim_ptr, const map_int_process_xsec &ProcessCrossSectionsMap, const int use_trust_level)
A function that sets the event weight based on the process cross-sections.
std::map< int, process_xsec_container > map_int_process_xsec
virtual double xsecErr_fb() const =0
Report the uncertainty in the total or process-specific cross section (in fb or pb).
double pow(const double &a)
Outputs a^i.
void setEventWeight_unity(EventWeighterFunctionType &result)
Module function providing an instance of EventWeighterFunctionType pointing to _setEventWeight_unity...
An abstract base class for collider simulators within ColliderBit.
std::function< void(HEPUtils_Event &, const BaseCollider *)> EventWeighterFunctionType
TODO: see if we can use this one:
Definition: Analysis.hpp:33
virtual double xsec_fb() const =0
void _setEventWeight_unity(HEPUtils::Event &event, const BaseCollider *)
A function that sets the event weight to unity, with zero uncertainty.