gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
xsec.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
20 
21 #include <algorithm>
22 #include <cassert>
23 #include <cmath>
24 #include <omp.h>
27 
28 namespace Gambit
29 {
30  namespace ColliderBit
31  {
32 
33 
37 
40  _xsec(0),
41  _xsecerr(0),
42  _info_string(""),
43  _trust_level(1)
44  { }
45 
48  {
49  _xsec = 0;
50  _xsecerr = 0;
51  _info_string = "";
52  _trust_level = 1;
53  }
54 
56  double xsec_container::operator()() const { return _xsec; }
57  double xsec_container::xsec() const { return _xsec; }
58 
60  double xsec_container::xsec_err() const { return _xsecerr; }
61 
63  double xsec_container::xsec_relerr() const { return _xsec > 0 ? _xsecerr/_xsec : 0; }
64 
66  void xsec_container::set_xsec(double xs, double xserr) { _xsec = xs; _xsecerr = xserr; }
67 
69  void xsec_container::average_xsec(double other_xsec, double other_xsecerr)
70  {
71  if (other_xsec > 0)
72  {
73  if (_xsec <= 0)
74  {
75  set_xsec(other_xsec, other_xsecerr);
76  }
77  else
78  {
79  double w = 1./(_xsecerr*_xsecerr);
80  double other_w = 1./(other_xsecerr*other_xsecerr);
81  _xsec = (w * _xsec + other_w * other_xsec) / (w + other_w);
82  _xsecerr = 1./sqrt(w + other_w);
83  }
84  }
85  }
86 
88  {
89  average_xsec(other(), other.xsec_err());
90  }
91 
93  void xsec_container::sum_xsecs(double other_xsec, double other_xsecerr)
94  {
95  if (other_xsec > 0)
96  {
97  if (_xsec <= 0)
98  {
99  set_xsec(other_xsec, other_xsecerr);
100  }
101  else
102  {
103  _xsec += other_xsec;
104  _xsecerr = sqrt(_xsecerr * _xsecerr + other_xsecerr * other_xsecerr);
105  }
106  }
107  }
108 
110  {
111  sum_xsecs(other(), other.xsec_err());
112  }
113 
115  std::map<std::string, double> xsec_container::get_content_as_map() const
116  {
117  std::map<std::string, double> content_map;
118 
119  std::string key;
120  std::string base_key;
121  if (_info_string != "") base_key = _info_string + "__";
122  else base_key = "";
123 
124  key = base_key + "xsec_" + unit;
125  content_map[key] = this->xsec();
126 
127  key = base_key + "xsec_err_" + unit;
128  content_map[key] = this->xsec_err();
129 
130  key = base_key + "xsec_relerr";
131  content_map[key] = this->xsec_relerr();
132 
133  content_map["trust_level"] = static_cast<double>(this->trust_level());
134 
135  return content_map;
136  }
137 
139  void xsec_container::set_info_string(std::string info_string_in) { _info_string = info_string_in; }
140 
142  std::string xsec_container::info_string() const { return _info_string; }
143 
145  void xsec_container::set_trust_level(int trust_level_in) { _trust_level = trust_level_in; }
146 
149 
151  const std::string xsec_container::unit = "fb";
152 
153 
154 
158 
162  _ntot(0)
163  { }
164 
167  {
169  _ntot = 0;
170 
171  // Add this instance to the instances map if it's not there already.
172  int thread = omp_get_thread_num();
173  if (instances_map.count(thread) == 0)
174  {
175  #pragma omp critical
176  {
177  instances_map[thread] = this;
178  }
179  }
180  }
181 
184 
186  long long MC_xsec_container::num_events() const { return _ntot; }
187 
189  double MC_xsec_container::xsec_per_event() const { return (_xsec >= 0 && _ntot > 0) ? _xsec/_ntot : 0; }
190 
193 
195  void MC_xsec_container::average_xsec(double other_xsec, double other_xsecerr, long long other_ntot)
196  {
197  // Run base class function
198  xsec_container::average_xsec(other_xsec, other_xsecerr);
199 
200  // Update _ntot
201  if (other_xsec > 0)
202  {
203  if (_xsec <= 0) set_num_events(other_ntot);
204  else _ntot += other_ntot;
205  }
206  }
207 
209  {
210  MC_xsec_container::average_xsec(other(), other.xsec_err(), other.num_events());
211  }
212 
214  void MC_xsec_container::sum_xsecs(double other_xsec, double other_xsecerr, long long other_ntot)
215  {
216  // Run base class function
217  xsec_container::sum_xsecs(other_xsec, other_xsecerr);
218 
219  // Update _ntot
220  if (other_xsec > 0)
221  {
222  if (_xsec <= 0) set_num_events(other_ntot);
223  else _ntot += other_ntot;
224  }
225  }
226 
228  {
229  MC_xsec_container::sum_xsecs(other(), other.xsec_err(), other.num_events());
230  }
231 
232 
235  {
236  int this_thread = omp_get_thread_num();
237  for (auto& thread_xsec_pair : instances_map)
238  {
239  if (thread_xsec_pair.first == this_thread) continue;
240  const MC_xsec_container& other = (*thread_xsec_pair.second);
241  average_xsec(other(), other.xsec_err(), other.num_events());
242  }
243  }
244 
247  {
248  int this_thread = omp_get_thread_num();
249  for (auto& thread_xsec_pair : instances_map)
250  {
251  if (thread_xsec_pair.first == this_thread) continue;
252  const MC_xsec_container& other = (*thread_xsec_pair.second);
253  _ntot += other.num_events();
254  }
255  }
256 
258  std::map<std::string, double> MC_xsec_container::get_content_as_map() const
259  {
260  // Get content from base class
261  std::map<std::string, double> content_map = xsec_container::get_content_as_map();
262 
263  // Add content specific to this class
264  std::string key;
265  std::string base_key;
266  if (_info_string != "") base_key = _info_string + "__";
267  else base_key = "";
268 
269  key = base_key + "xsec_per_event_" + unit;
270  content_map[key] = this->xsec_per_event();
271 
272  key = base_key + "logged_events_" + unit;
273  content_map[key] = _ntot;
274 
275  return content_map;
276  }
277 
279  std::map<int, const MC_xsec_container*> MC_xsec_container::instances_map;
280 
281 
282 
286 
290  _process_code(-1),
291  _processes_sharing_xsec(std::vector<int>()),
292  _related_pid_pairs(std::vector<PID_pair>())
293  { }
294 
297  {
299  _process_code = -1;
300  _processes_sharing_xsec.clear();
301  _related_pid_pairs.clear();
302  }
303 
305  void process_xsec_container::average_xsec(double other_xsec, double other_xsecerr)
306  {
307  // Run base class function
308  xsec_container::average_xsec(other_xsec, other_xsecerr);
309  }
311  {
312  // Check that the process code of this instance is set
313  assert(_process_code > 0);
314  // Check that we are working with the same process code
315  assert(other.process_code() == _process_code);
316  // @todo Should we also check the content of the vectors
317  // _processes_sharing_xsec and _related_pid_pairs?
319  }
320 
322  void process_xsec_container::sum_xsecs(double other_xsec, double other_xsecerr)
323  {
324  // Run base class function
325  xsec_container::sum_xsecs(other_xsec, other_xsecerr);
326  }
328  {
329  // Check that the process code of this instance is set
330  assert(_process_code > 0);
331  // Check that we are working with the same process code
332  assert(other.process_code() == _process_code);
333  // @todo Should we also check the content of the vectors
334  // _processes_sharing_xsec and _related_pid_pairs?
336  }
337 
338 
341  { return _process_code; }
342 
344  void process_xsec_container::set_process_code(int process_code_in)
345  { _process_code = process_code_in; }
346 
350  const std::vector<int>& process_xsec_container::processes_sharing_xsec() const
351  { return _processes_sharing_xsec; }
352 
355  {
356  if(std::find(_processes_sharing_xsec.begin(), _processes_sharing_xsec.end(), process_code_in) == _processes_sharing_xsec.end())
357  {
358  _processes_sharing_xsec.push_back(process_code_in);
359  }
360  }
361 
363  const std::vector<PID_pair>& process_xsec_container::related_pid_pairs() const
364  { return _related_pid_pairs; }
365 
368  {
369  if(std::find(_related_pid_pairs.begin(), _related_pid_pairs.end(), pid_pair_in) == _related_pid_pairs.end())
370  {
371  _related_pid_pairs.push_back(pid_pair_in);
372  }
373  }
374 
375 
376 
380 
384  _pid_pair(PID_pair()),
385  _pid_pairs_sharing_xsec(std::vector<PID_pair>()),
386  _related_processes(std::vector<int>())
387  { }
388 
391  {
393  _pid_pair = PID_pair();
394  _pid_pairs_sharing_xsec.clear();
395  _related_processes.clear();
396  }
397 
399  void PID_pair_xsec_container::average_xsec(double other_xsec, double other_xsecerr)
400  {
401  // Run base class function
402  xsec_container::average_xsec(other_xsec, other_xsecerr);
403  }
405  {
406  // Check that the PID pair of this instance is set
407  assert((_pid_pair.pid1() != 0) && (_pid_pair.pid2() != 0));
408  // Check that we are working with the same PID pair
409  assert(other.pid_pair() == _pid_pair);
410  // @todo Should we also check the content of the vectors
411  // _pid_pairs_sharing_xsec and _related_processes?
413  }
414 
416  void PID_pair_xsec_container::sum_xsecs(double other_xsec, double other_xsecerr)
417  {
418  // Run base class function
419  xsec_container::sum_xsecs(other_xsec, other_xsecerr);
420  }
422  {
423  // Check that the PID pair of this instance is set
424  assert((_pid_pair.pid1() != 0) && (_pid_pair.pid2() != 0));
425  // Check that we are working with the same PID pair
426  assert(other.pid_pair() == _pid_pair);
427  // @todo Should we also check the content of the vectors
428  // _pid_pairs_sharing_xsec and _related_processes?
430  }
431 
434  { return _pid_pair; }
435 
438  { _pid_pair = pid_pair_in; }
439 
443  const std::vector<PID_pair>& PID_pair_xsec_container::pid_pairs_sharing_xsec() const
444  { return _pid_pairs_sharing_xsec; }
445 
448  {
449  if(std::find(_pid_pairs_sharing_xsec.begin(), _pid_pairs_sharing_xsec.end(), pid_pair_in) == _pid_pairs_sharing_xsec.end())
450  {
451  _pid_pairs_sharing_xsec.push_back(pid_pair_in);
452  }
453  }
454 
456  const std::vector<int>& PID_pair_xsec_container::related_processes() const
457  { return _related_processes; }
458 
461  {
462  if(std::find(_related_processes.begin(), _related_processes.end(), process_code_in) == _related_processes.end())
463  {
464  _related_processes.push_back(process_code_in);
465  }
466  }
467 
468 
469  }
470 }
A base class for holding cross-section info within ColliderBit.
Definition: xsec.hpp:37
std::string info_string() const
Get the info string.
Definition: xsec.cpp:142
double xsec_err() const
Return the cross-section error (in fb).
Definition: xsec.cpp:60
void average_xsec(double, double, long long)
Average cross-sections and combine errors.
Definition: xsec.cpp:195
int trust_level() const
Get the trust level.
Definition: xsec.cpp:148
const PID_pair & pid_pair() const
Return the PID pair.
Definition: xsec.cpp:433
const std::vector< PID_pair > & pid_pairs_sharing_xsec() const
Return the list of PID pairs that share this cross-section (This is due to the many-to-many mapping b...
Definition: xsec.cpp:443
void reset()
Reset this instance for reuse.
Definition: xsec.cpp:296
A class for holding the production cross-section for final state identified by the pair of PID codes...
Definition: xsec.hpp:195
void register_related_process(int)
Add a process code to the list of processes related to this cross-section.
Definition: xsec.cpp:460
void set_xsec(double, double)
Set the cross-section and its error (in fb).
Definition: xsec.cpp:66
Simple class for holding a sorted pair of particle ID (PID) codes.
Definition: PID_pair.hpp:31
void set_num_events(long long)
Set the total number of events seen so far.
Definition: xsec.cpp:192
void set_process_code(int)
Set the process code.
Definition: xsec.cpp:344
void average_xsec(double, double)
Average cross-sections and combine errors.
Definition: xsec.cpp:399
STL namespace.
int process_code() const
Return the process code.
Definition: xsec.cpp:340
process_xsec_container()
Definitions of process_xsec_container members.
Definition: xsec.cpp:288
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
int pid1() const
Definition: PID_pair.hpp:84
A class for holding a total cross-section calculated via MC across multiple threads.
Definition: xsec.hpp:96
std::vector< PID_pair > _pid_pairs_sharing_xsec
Definition: xsec.hpp:235
void register_process_sharing_xsec(int)
Add a process code to the list of processes sharing this cross-section.
Definition: xsec.cpp:354
void log_event()
Tell the xsec object that there has been a new event.
Definition: xsec.cpp:183
void average_xsec(double, double)
Average cross-sections and combine errors.
Definition: xsec.cpp:305
ColliderBit (production) cross-section class.
void register_related_pid_pair(PID_pair)
Add a PID pair to the list of processes related to this cross-section.
Definition: xsec.cpp:367
MC_xsec_container()
Definitions of MC_xsec_container members.
Definition: xsec.cpp:160
void average_xsec(double, double)
Average cross-sections and combine errors.
Definition: xsec.cpp:69
void gather_num_events()
Collect total events seen on all threads.
Definition: xsec.cpp:246
static std::map< int, const MC_xsec_container * > instances_map
A map with pointers to all instances of this class. The key is the OMP thread number.
Definition: xsec.hpp:142
void set_pid_pair(const PID_pair &)
Set the PID pair.
Definition: xsec.cpp:437
PID_pair_xsec_container()
Definitions of PID_pair_xsec_container members.
Definition: xsec.cpp:382
std::vector< PID_pair > _related_pid_pairs
Definition: xsec.hpp:189
std::map< std::string, double > get_content_as_map() const
Get content as map <string,double> map (for easy printing).
Definition: xsec.cpp:115
A class for holding the cross-section of a single Pythia process (identified by the Pythia process co...
Definition: xsec.hpp:148
int pid2() const
Definition: PID_pair.hpp:89
void register_pid_pair_sharing_xsec(PID_pair)
Add a PID pair to the list of PID pairs sharing this cross-section.
Definition: xsec.cpp:447
void sum_xsecs(double, double)
Sum cross-sections and add errors in quadrature.
Definition: xsec.cpp:93
void sum_xsecs(double, double, long long)
Sum cross-sections and add errors in quadrature.
Definition: xsec.cpp:214
void reset()
Reset this instance for reuse.
Definition: xsec.cpp:47
void reset()
Reset this instance for reuse.
Definition: xsec.cpp:390
DS5_MSPCTM DS_INTDOF int
void sum_xsecs(double, double)
Sum cross-sections and add errors in quadrature.
Definition: xsec.cpp:322
std::map< std::string, double > get_content_as_map() const
Get content as map <string,double> map (for easy printing).
Definition: xsec.cpp:258
xsec_container()
Definitions of xsec members.
Definition: xsec.cpp:39
const std::vector< int > & related_processes() const
Return the list of process codes related to this cross-section.
Definition: xsec.cpp:456
const std::vector< PID_pair > & related_pid_pairs() const
Return the list of PID pairs related to this cross-section.
Definition: xsec.cpp:363
Exception objects required for standalone compilation.
void gather_xsecs()
Collect xsec predictions from other threads and do a weighted combination.
Definition: xsec.cpp:234
long long num_events() const
Return the total number of events seen so far.
Definition: xsec.cpp:186
static const std::string unit
String Let&#39;s make it clear that we work with fb as unit.
Definition: xsec.hpp:84
void reset()
Reset this instance for reuse.
Definition: xsec.cpp:166
void set_info_string(std::string)
Set the info string.
Definition: xsec.cpp:139
double xsec_relerr() const
Return the cross-section relative error.
Definition: xsec.cpp:63
void set_trust_level(int)
Set the trust level.
Definition: xsec.cpp:145
double operator()() const
Return the full cross-section (in fb).
Definition: xsec.cpp:56
std::vector< int > _processes_sharing_xsec
Definition: xsec.hpp:188
double xsec_per_event() const
Return the cross-section per event seen (in fb).
Definition: xsec.cpp:189
TODO: see if we can use this one:
Definition: Analysis.hpp:33
void sum_xsecs(double, double)
Sum cross-sections and add errors in quadrature.
Definition: xsec.cpp:416