gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
LHC_likelihoods.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
36 
37 #include <string>
38 #include <sstream>
39 
42 
43 #include "Eigen/Eigenvalues"
44 #include <gsl/gsl_sf_gamma.h>
45 
46 // #define COLLIDERBIT_DEBUG
47 
48 namespace Gambit
49 {
50 
51  namespace ColliderBit
52  {
53 
56  {
57  using namespace Pipes::calc_LHC_signals;
58 
59  // Clear the result map
60  result.clear();
61 
62  std::stringstream summary_line;
63  summary_line << "LHC signals per SR: ";
64 
65  // Loop over analyses and collect the predicted events into the map
66  for (size_t analysis = 0; analysis < Dep::AllAnalysisNumbers->size(); ++analysis)
67  {
68  // AnalysisData for this analysis
69  const AnalysisData& adata = *(Dep::AllAnalysisNumbers->at(analysis));
70 
71  summary_line << adata.analysis_name << ": ";
72 
73  // Loop over the signal regions inside the analysis, and save the predicted number of events for each.
74  for (size_t SR = 0; SR < adata.size(); ++SR)
75  {
76  // Save SR numbers and absolute uncertainties
77  const SignalRegionData srData = adata[SR];
78  const str key = adata.analysis_name + "__" + srData.sr_label + "__i" + std::to_string(SR) + "__signal";
79  result[key] = srData.n_signal_at_lumi;
80  const double abs_uncertainty_s_stat = (srData.n_signal == 0 ? 0 : sqrt(srData.n_signal) * (srData.n_signal_at_lumi/srData.n_signal));
81  const double abs_uncertainty_s_sys = srData.signal_sys;
82  const double combined_uncertainty = HEPUtils::add_quad(abs_uncertainty_s_stat, abs_uncertainty_s_sys);
83  result[key + "_uncert"] = combined_uncertainty;
84 
85  summary_line << srData.sr_label + "__i" + std::to_string(SR) << ":" << srData.n_signal_at_lumi << "+-" << combined_uncertainty << ", ";
86  }
87  }
88  logger() << LogTags::debug << summary_line.str() << EOM;
89  }
90 
91 
94  {
95  using namespace Pipes::calc_LHC_LogLikes;
96 
97  // Use covariance matrix when available?
98  static const bool use_covar = runOptions->getValueOrDef<bool>(true, "use_covariances");
99 
100  // Clear the result map
101  result.clear();
102 
103  // Loop over analyses and calculate the observed dLL for each
104  for (size_t analysis = 0; analysis < Dep::AllAnalysisNumbers->size(); ++analysis)
105  {
106  // AnalysisData for this analysis
107  const AnalysisData& adata = *(Dep::AllAnalysisNumbers->at(analysis));
108 
114  if (not Dep::RunMC->event_generation_began || Dep::RunMC->exceeded_maxFailedEvents)
115  {
116  // If this is an anlysis with covariance info, only add a single 0-entry in the map
117  if (use_covar && adata.srcov.rows() > 0)
118  {
119  result[adata.analysis_name].combination_sr_label = "none";
120  result[adata.analysis_name].combination_loglike = 0.0;
121  continue;
122  }
123  // If this is an anlysis without covariance info, add 0-entries for all SRs plus
124  // one for the combined LogLike
125  else
126  {
127  for (size_t SR = 0; SR < adata.size(); ++SR)
128  {
129  result[adata.analysis_name].sr_indices[adata[SR].sr_label] = SR;
130  result[adata.analysis_name].sr_loglikes[adata[SR].sr_label] = 0.0;
131  continue;
132  }
133  result[adata.analysis_name].combination_sr_label = "none";
134  result[adata.analysis_name].combination_loglike = 0.0;
135  continue;
136  }
137 
138  }
139 
140 
141  #ifdef COLLIDERBIT_DEBUG
142  std::streamsize stream_precision = cout.precision(); // get current precision
143  cout.precision(2); // set precision
144  cout << debug_prefix() << "calc_LHC_LogLikes: " << "Will print content of " << adata.analysis_name << " signal regions:" << endl;
145  for (size_t SR = 0; SR < adata.size(); ++SR)
146  {
147  const SignalRegionData& srData = adata[SR];
148  cout << std::fixed << debug_prefix()
149  << "calc_LHC_LogLikes: " << adata.analysis_name
150  << ", " << srData.sr_label
151  << ", n_b = " << srData.n_background << " +/- " << srData.background_sys
152  << ", n_obs = " << srData.n_observed
153  << ", excess = " << srData.n_observed - srData.n_background << " +/- " << srData.background_sys
154  << ", n_s = " << srData.n_signal_at_lumi
155  << ", (excess-n_s) = " << (srData.n_observed-srData.n_background) - srData.n_signal_at_lumi << " +/- " << srData.background_sys
156  << ", n_s_MC = " << srData.n_signal
157  << endl;
158  }
159  cout.precision(stream_precision); // restore previous precision
160  #endif
161 
162 
163  // Loop over the signal regions inside the analysis, and work out the total (delta) log likelihood for this analysis
166  if (use_covar && adata.srcov.rows() > 0)
167  {
184 
185  #ifdef COLLIDERBIT_DEBUG
186  cout << debug_prefix() << "calc_LHC_LogLikes: Analysis " << analysis << " has a covariance matrix: computing composite loglike." << endl;
187  #endif
188 
189 
190  // Shortcut: if all SRs have 0 signal prediction, we know the Delta LogLike is 0.
191  bool all_zero_signal = true;
192  for (size_t SR = 0; SR < adata.size(); ++SR)
193  {
194  if (adata[SR].n_signal != 0)
195  {
196  all_zero_signal = false;
197  break;
198  }
199  }
200  if (all_zero_signal)
201  {
202  // Store result
203  result[adata.analysis_name].combination_sr_label = "all";
204  result[adata.analysis_name].combination_sr_index = -1;
205  result[adata.analysis_name].combination_loglike = 0.0;
206 
207  #ifdef COLLIDERBIT_DEBUG
208  cout << debug_prefix() << "calc_LHC_LogLikes: " << adata.analysis_name << "_LogLike : " << 0.0 << " (No signal predicted. Skipped covariance calculation.)" <<endl;
209  #endif
210 
211  // Continue to next analysis
212  continue;
213  }
214 
215  // Construct vectors of SR numbers
216  Eigen::ArrayXd n_obs(adata.size()), logfact_n_obs(adata.size()), n_pred_b(adata.size()), n_pred_sb(adata.size()), abs_unc_s(adata.size());
217  for (size_t SR = 0; SR < adata.size(); ++SR)
218  {
219  const SignalRegionData srData = adata[SR];
220 
221  // Actual observed number of events
222  n_obs(SR) = srData.n_observed;
223 
224  // Log factorial of observed number of events.
225  // Currently use the ln(Gamma(x)) function gsl_sf_lngamma from GSL. (Need continuous function.)
226  // We may want to switch to using Sterlings approximation: ln(n!) ~ n*ln(n) - n
227  logfact_n_obs(SR) = gsl_sf_lngamma(n_obs(SR) + 1.);
228 
229  // A contribution to the predicted number of events that is not known exactly
230  n_pred_b(SR) = srData.n_background;
231  n_pred_sb(SR) = srData.n_signal_at_lumi + srData.n_background;
232 
233  // Absolute errors for n_predicted_uncertain_*
234  const double abs_uncertainty_s_stat = (srData.n_signal == 0 ? 0 : sqrt(srData.n_signal) * (srData.n_signal_at_lumi/srData.n_signal));
235  const double abs_uncertainty_s_sys = srData.signal_sys;
236  abs_unc_s(SR) = HEPUtils::add_quad(abs_uncertainty_s_stat, abs_uncertainty_s_sys);
237  }
238 
239  // Diagonalise the background-only covariance matrix, extracting the rotation matrix
243  const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig_b(adata.srcov);
244  const Eigen::ArrayXd Eb = eig_b.eigenvalues();
245  const Eigen::ArrayXd sqrtEb = Eb.sqrt();
246  const Eigen::MatrixXd Vb = eig_b.eigenvectors();
247  //const Eigen::MatrixXd Vbinv = Vb.inverse();
248 
249  // Construct and diagonalise the s+b covariance matrix, adding the diagonal signal uncertainties in quadrature
251  const Eigen::MatrixXd srcov_s = abs_unc_s.array().square().matrix().asDiagonal();
252  const Eigen::MatrixXd srcov_sb = adata.srcov + srcov_s;
253  const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig_sb(srcov_sb);
254  const Eigen::ArrayXd Esb = eig_sb.eigenvalues();
255  const Eigen::ArrayXd sqrtEsb = Esb.sqrt();
256  const Eigen::MatrixXd Vsb = eig_sb.eigenvectors();
257  //const Eigen::MatrixXd Vsbinv = Vsb.inverse();
258 
261 
262  // Sample correlated SR rates from a rotated Gaussian defined by the covariance matrix and offset by the mean rates
263  static const double CONVERGENCE_TOLERANCE_ABS = runOptions->getValueOrDef<double>(0.05, "covariance_marg_convthres_abs");
264  static const double CONVERGENCE_TOLERANCE_REL = runOptions->getValueOrDef<double>(0.05, "covariance_marg_convthres_rel");
265  static const size_t nsample_input = runOptions->getValueOrDef<size_t>(100000, "covariance_nsamples_start");
266  size_t NSAMPLE = nsample_input;
267 
268  // Dynamic convergence control & test variables
269  bool first_iteration = true;
270  double diff_abs = 9999;
271  double diff_rel = 1;
272 
273  // Likelihood variables (note use of long double to guard against blow-up of L as opposed to log(L1/L0))
274  long double ana_like_b_prev = 1;
275  long double ana_like_sb_prev = 1;
276  long double ana_like_b = 1;
277  long double ana_like_sb = 1;
278  long double lsum_b_prev = 0;
279  long double lsum_sb_prev = 0;
280 
281  std::normal_distribution<double> unitnormdbn(0,1);
282 
283  // Check absolute difference between independent estimates
285  while ((diff_abs > CONVERGENCE_TOLERANCE_ABS && diff_rel > CONVERGENCE_TOLERANCE_REL) || 1.0/sqrt(NSAMPLE) > CONVERGENCE_TOLERANCE_ABS)
286  {
287  long double lsum_b = 0;
288  long double lsum_sb = 0;
289 
290  // typedef Eigen::Array<long double, Eigen::Dynamic, 1> ArrayXld;
291 
299 
300  const bool COVLOGNORMAL = false;
301  if (!COVLOGNORMAL)
302  {
303 
304  #pragma omp parallel
305  {
306  double lsum_b_private = 0;
307  double lsum_sb_private = 0;
308 
309  // Sample correlated SR rates from a rotated Gaussian defined by the covariance matrix and offset by the mean rates
310  #pragma omp for nowait
311  for (size_t i = 0; i < NSAMPLE; ++i) {
312 
313  Eigen::VectorXd norm_sample_b(adata.size()), norm_sample_sb(adata.size());
314  for (size_t j = 0; j < adata.size(); ++j) {
315  norm_sample_b(j) = sqrtEb(j) * unitnormdbn(Random::rng());
316  norm_sample_sb(j) = sqrtEsb(j) * unitnormdbn(Random::rng());
317  }
318 
319  // Rotate rate deltas into the SR basis and shift by SR mean rates
320  const Eigen::VectorXd n_pred_b_sample = n_pred_b + (Vb*norm_sample_b).array();
321  const Eigen::VectorXd n_pred_sb_sample = n_pred_sb + (Vsb*norm_sample_sb).array();
322 
323  // Calculate Poisson likelihood and add to composite likelihood calculation
324  double combined_loglike_b = 0;
325  double combined_loglike_sb = 0;
326  for (size_t j = 0; j < adata.size(); ++j) {
327  const double lambda_b_j = std::max(n_pred_b_sample(j), 1e-3); //< manually avoid <= 0 rates
328  const double lambda_sb_j = std::max(n_pred_sb_sample(j), 1e-3); //< manually avoid <= 0 rates
329  const double loglike_b_j = n_obs(j)*log(lambda_b_j) - lambda_b_j - logfact_n_obs(j);
330  const double loglike_sb_j = n_obs(j)*log(lambda_sb_j) - lambda_sb_j - logfact_n_obs(j);
331  combined_loglike_b += loglike_b_j;
332  combined_loglike_sb += loglike_sb_j;
333  }
334  // Add combined likelihood to running sums (to later calculate averages)
335  lsum_b_private += exp(combined_loglike_b);
336  lsum_sb_private += exp(combined_loglike_sb);
337  }
338  #pragma omp critical
339  {
340  lsum_b += lsum_b_private;
341  lsum_sb += lsum_sb_private;
342  }
343  } // End omp parallel
344  } // End if !COVLOGNORMAL
345 
346  // /// @todo Check that this log-normal sampling works as expected.
347  // else // COVLOGNORMAL
348  // {
349 
350  // const Eigen::ArrayXd ln_n_pred_b = n_pred_b.log();
351  // const Eigen::ArrayXd ln_n_pred_sb = n_pred_sb.log();
352  // const Eigen::ArrayXd ln_sqrtEb = (n_pred_b + sqrtEb).log() - ln_n_pred_b;
353  // const Eigen::ArrayXd ln_sqrtEsb = (n_pred_sb + sqrtEsb).log() - ln_n_pred_sb;
354 
355  // #pragma omp parallel
356  // {
357  // std::normal_distribution<> unitnormdbn{0,1};
358  // Eigen::ArrayXd llrsums_private = Eigen::ArrayXd::Zero(adata.size());
359 
360  // #pragma omp for nowait
361 
362  // // Sample correlated SR rates from a rotated Gaussian defined by the covariance matrix and offset by the mean rates
363  // for (size_t i = 0; i < NSAMPLE; ++i) {
364  // Eigen::VectorXd ln_norm_sample_b(adata.size()), ln_norm_sample_sb(adata.size());
365  // for (size_t j = 0; j < adata.size(); ++j) {
366  // ln_norm_sample_b(j) = ln_sqrtEb(j) * unitnormdbn(Random::rng());
367  // ln_norm_sample_sb(j) = ln_sqrtEsb(j) * unitnormdbn(Random::rng());
368  // }
369 
370  // // Rotate rate deltas into the SR basis and shift by SR mean rates
371  // const Eigen::ArrayXd delta_ln_n_pred_b_sample = Vb*ln_norm_sample_b;
372  // const Eigen::ArrayXd delta_ln_n_pred_sb_sample = Vsb*ln_norm_sample_sb;
373  // const Eigen::ArrayXd n_pred_b_sample = (ln_n_pred_b + delta_ln_n_pred_b_sample).exp();
374  // const Eigen::ArrayXd n_pred_sb_sample = (ln_n_pred_sb + delta_ln_n_pred_sb_sample).exp();
375 
376  // // Calculate Poisson LLR and add to aggregated LL calculation
377  // for (size_t j = 0; j < adata.size(); ++j) {
378  // const double lambda_b_j = std::max(n_pred_b_sample(j), 1e-3); //< shouldn't be needed in log-space sampling
379  // const double lambda_sb_j = std::max(n_pred_sb_sample(j), 1e-3); //< shouldn't be needed in log-space sampling
380  // const double llr_j = n_obs(j)*log(lambda_sb_j/lambda_b_j) - (lambda_sb_j - lambda_b_j);
381  // llrsums_private(j) += llr_j;
382  // }
383  // }
384 
385  // #pragma omp critical
386  // {
387  // for (size_t j = 0; j < adata.size(); ++j) { llrsums(j) += llrsums_private(j); }
388  // }
389  // } // End omp parallel
390  // }
391 
392  // Compare convergence to previous independent batch
393  if (first_iteration) // The first round must be generated twice
394  {
395  lsum_b_prev = lsum_b;
396  lsum_sb_prev = lsum_sb;
397  first_iteration = false;
398  }
399  else
400  {
401  ana_like_b_prev = lsum_b_prev / (double)NSAMPLE;
402  ana_like_sb_prev = lsum_sb_prev / (double)NSAMPLE;
403  ana_like_b = lsum_b / (double)NSAMPLE;
404  ana_like_sb = lsum_sb / (double)NSAMPLE;
405  //
406  const double diff_abs_b = std::abs(ana_like_b_prev - ana_like_b);
407  const double diff_abs_sb = std::abs(ana_like_sb_prev - ana_like_sb);
408  const double diff_rel_b = diff_abs_b/ana_like_b;
409  const double diff_rel_sb = diff_abs_sb/ana_like_sb;
410  //
411  diff_rel = std::max(diff_rel_b, diff_rel_sb); // Relative convergence check
412  diff_abs = std::max(diff_abs_b, diff_abs_sb); // Absolute convergence check
413 
414  // Update variables
415  lsum_b_prev += lsum_b; // Aggregate result. This doubles the effective batch size for lsum_prev.
416  lsum_sb_prev += lsum_sb; // Aggregate result. This doubles the effective batch size for lsum_prev.
417  NSAMPLE *=2; // This ensures that the next batch for lsum is as big as the current batch size for lsum_prev, so they can be compared directly.
418  }
419 
420  #ifdef COLLIDERBIT_DEBUG
421  cout << debug_prefix()
422  << "diff_rel: " << diff_rel << endl
423  << " diff_abs: " << diff_abs << endl
424  << " ana_llr_prev: " << log(ana_like_sb_prev/ana_like_b_prev) << endl
425  << " ana_dll: " << log(ana_like_sb/ana_like_b) << endl
426  << " logl_sb: " << log(ana_like_sb) << endl
427  << " logl_b: " << log(ana_like_b) << endl;
428  cout << debug_prefix() << "NSAMPLE for the next iteration is: " << NSAMPLE << endl;
429  cout << debug_prefix() << endl;
430  #endif
431  } // End while loop
432 
433  // Combine the independent estimates ana_like and ana_like_prev.
434  // Use equal weights since the estimates are based on equal batch sizes.
435  ana_like_b = 0.5*(ana_like_b + ana_like_b_prev);
436  ana_like_sb = 0.5*(ana_like_sb + ana_like_sb_prev);
437 
438  // Compute LLR from mean s+b and b likelihoods
439  const double ana_dll = log(ana_like_sb) - log(ana_like_b);
440  #ifdef COLLIDERBIT_DEBUG
441  cout << debug_prefix() << "Combined estimate: ana_dll: " << ana_dll << " (based on 2*NSAMPLE=" << 2*NSAMPLE << " samples)" << endl;
442  #endif
443 
444  // Check for problem
445  if (Utils::isnan(ana_dll))
446  {
447  std::stringstream msg;
448  msg << "Computation of composite loglike for analysis " << adata.analysis_name << " returned NaN. Will now print the signal region data for this analysis:" << endl;
449  for (size_t SR = 0; SR < adata.size(); ++SR)
450  {
451  const SignalRegionData& srData = adata[SR];
452  msg << srData.sr_label
453  << ", n_background = " << srData.n_background
454  << ", background_sys = " << srData.background_sys
455  << ", n_observed = " << srData.n_observed
456  << ", n_signal_at_lumi = " << srData.n_signal_at_lumi
457  << ", n_signal = " << srData.n_signal
458  << ", signal_sys = " << srData.signal_sys
459  << endl;
460  }
461  invalid_point().raise(msg.str());
462  }
463 
464  // Store result
465  result[adata.analysis_name].combination_sr_label = "all";
466  result[adata.analysis_name].combination_sr_index = -1;
467  result[adata.analysis_name].combination_loglike = ana_dll;
468 
469  #ifdef COLLIDERBIT_DEBUG
470  cout << debug_prefix() << "calc_LHC_LogLikes: " << adata.analysis_name << "_LogLike : " << ana_dll << endl;
471  #endif
472 
473  }
474 
475  else
476  {
477  // No SR-correlation info, or user chose not to use it.
478  // Then we either take the result from the SR *expected* to be most constraining
479  // under the s=0 assumption (default), or naively combine the loglikes for
480  // all SRs (if combine_SRs_without_covariances=true).
481  #ifdef COLLIDERBIT_DEBUG
482  cout << debug_prefix() << "calc_LHC_LogLikes: Analysis " << analysis << " has no covariance matrix: computing single best-expected loglike." << endl;
483  #endif
484 
485  double bestexp_dll_exp = 0, bestexp_dll_obs = 0;
486  str bestexp_sr_label;
487  int bestexp_sr_index;
488  double nocovar_srsum_dll_obs = 0;
489 
490  for (size_t SR = 0; SR < adata.size(); ++SR)
491  {
492  const SignalRegionData &srData = adata[SR];
493 
494  // Actual observed number of events
495  const int n_obs = (int) round(srData.n_observed);
496 
497  // A contribution to the predicted number of events that is known exactly
498  // (e.g. from data-driven background estimate)
499  const double n_predicted_exact = 0;
500 
501  // A contribution to the predicted number of events that is not known exactly
502  const double n_predicted_uncertain_s = srData.n_signal_at_lumi;
503  const double n_predicted_uncertain_b = srData.n_background;
504  const double n_predicted_uncertain_sb = n_predicted_uncertain_s + n_predicted_uncertain_b;
505 
506  // Absolute errors for n_predicted_uncertain_*
507  const double abs_uncertainty_s_stat = (srData.n_signal == 0 ? 0 : sqrt(srData.n_signal) * (srData.n_signal_at_lumi/srData.n_signal));
508  const double abs_uncertainty_s_sys = srData.signal_sys;
509  const double abs_uncertainty_b = srData.background_sys;
510  const double abs_uncertainty_sb = HEPUtils::add_quad(abs_uncertainty_s_stat, abs_uncertainty_s_sys, abs_uncertainty_b);
511 
512  // Relative errors for n_predicted_uncertain_*
513  const double frac_uncertainty_b = abs_uncertainty_b / n_predicted_uncertain_b;
514  const double frac_uncertainty_sb = abs_uncertainty_sb / n_predicted_uncertain_sb;
515 
516  // Predicted total background, as an integer for use in Poisson functions
517  const int n_predicted_total_b_int = (int) round(n_predicted_exact + n_predicted_uncertain_b);
518 
519  // Marginalise over systematic uncertainties on mean rates
520  // Use a log-normal/Gaussia distribution for the nuisance parameter, as requested
521  auto marginaliser = (*BEgroup::lnlike_marg_poisson == "lnlike_marg_poisson_lognormal_error")
522  ? BEreq::lnlike_marg_poisson_lognormal_error : BEreq::lnlike_marg_poisson_gaussian_error;
523  const double llb_exp = marginaliser(n_predicted_total_b_int, n_predicted_exact, n_predicted_uncertain_b, frac_uncertainty_b);
524  const double llsb_exp = marginaliser(n_predicted_total_b_int, n_predicted_exact, n_predicted_uncertain_sb, frac_uncertainty_sb);
525  const double llb_obs = marginaliser(n_obs, n_predicted_exact, n_predicted_uncertain_b, frac_uncertainty_b);
526  const double llsb_obs = marginaliser(n_obs, n_predicted_exact, n_predicted_uncertain_sb, frac_uncertainty_sb);
527 
528  // Calculate the expected dll and set the bestexp values for exp and obs dll if this one is the best so far
529  const double dll_exp = llb_exp - llsb_exp; //< note positive dll convention -> more exclusion here
530  #ifdef COLLIDERBIT_DEBUG
531  cout << debug_prefix() << adata.analysis_name << ", " << srData.sr_label << ", llsb_exp-llb_exp = " << llsb_exp-llb_exp << ", llsb_obs-llb_obs= " << llsb_obs - llb_obs << endl;
532  #endif
533  if (dll_exp > bestexp_dll_exp || SR == 0)
534  {
535  bestexp_dll_exp = dll_exp;
536  bestexp_dll_obs = llb_obs - llsb_obs;
537  bestexp_sr_label = srData.sr_label;
538  bestexp_sr_index = SR;
539  // #ifdef COLLIDERBIT_DEBUG
540  // cout << debug_prefix() << "Setting bestexp_sr_label to: " << bestexp_sr_label << ", LogL_exp = " << -bestexp_dll_exp << ", LogL_obs = " << -bestexp_dll_obs << endl;
541  // #endif
542  }
543 
544  // Store "observed LogLike" result for this SR
545  result[adata.analysis_name].sr_indices[srData.sr_label] = SR;
546  result[adata.analysis_name].sr_loglikes[srData.sr_label] = llsb_obs - llb_obs;
547 
548  // Add loglike to the no-correlations loglike sum over SRs
549  nocovar_srsum_dll_obs += llsb_obs - llb_obs;
550  }
551 
552  // Check for problem
553  if (Utils::isnan(bestexp_dll_obs))
554  {
555  std::stringstream msg;
556  msg << "Computation of single-SR loglike for analysis " << adata.analysis_name << " returned NaN, from signal region: " << bestexp_sr_label << endl;
557  msg << "Will now print the signal region data for this analysis:" << endl;
558  for (size_t SR = 0; SR < adata.size(); ++SR)
559  {
560  const SignalRegionData& srData = adata[SR];
561  msg << srData.sr_label
562  << ", n_background = " << srData.n_background
563  << ", background_sys = " << srData.background_sys
564  << ", n_observed = " << srData.n_observed
565  << ", n_signal_at_lumi = " << srData.n_signal_at_lumi
566  << ", n_signal = " << srData.n_signal
567  << ", signal_sys = " << srData.signal_sys
568  << endl;
569  }
570  invalid_point().raise(msg.str());
571  }
572 
573  // Set this analysis' total obs dLL to that from the best-expected SR (with conversion to more negative dll = more exclusion convention)
574  // result[adata.analysis_name] = -bestexp_dll_obs;
575  result[adata.analysis_name].combination_sr_label = bestexp_sr_label;
576  result[adata.analysis_name].combination_sr_index = bestexp_sr_index;
577  result[adata.analysis_name].combination_loglike = -bestexp_dll_obs;
578 
579  // Should we use the naive sum of SR loglikes (without correlations), instead of the best-expected SR?
580  static const bool combine_nocovar_SRs = runOptions->getValueOrDef<bool>(false, "combine_SRs_without_covariances");
581  if (combine_nocovar_SRs)
582  {
583  result[adata.analysis_name].combination_loglike = nocovar_srsum_dll_obs;
584  }
585 
586  #ifdef COLLIDERBIT_DEBUG
587  cout << debug_prefix() << "calc_LHC_LogLikes: " << adata.analysis_name << "_" << bestexp_sr_label << "_LogLike : " << -bestexp_dll_obs << endl;
588  #endif
589  }
590 
591  }
592 
593  }
594 
595 
598  {
599  using namespace Pipes::get_LHC_LogLike_per_analysis;
600 
601  std::stringstream summary_line;
602  summary_line << "LHC loglikes per analysis: ";
603 
604  for (const std::pair<str,AnalysisLogLikes>& pair : *Dep::LHC_LogLikes)
605  {
606  const str& analysis_name = pair.first;
607  const AnalysisLogLikes& analysis_loglikes = pair.second;
608 
609  result[analysis_name] = analysis_loglikes.combination_loglike;
610 
611  summary_line << analysis_name << ":" << analysis_loglikes.combination_loglike << ", ";
612  }
613  logger() << LogTags::debug << summary_line.str() << EOM;
614  }
615 
616 
619  {
620  using namespace Pipes::get_LHC_LogLike_per_SR;
621 
622  std::stringstream summary_line;
623  summary_line << "LHC loglikes per SR: ";
624 
625  for (const std::pair<str,AnalysisLogLikes>& pair_i : *Dep::LHC_LogLikes)
626  {
627  const str& analysis_name = pair_i.first;
628  const AnalysisLogLikes& analysis_loglikes = pair_i.second;
629 
630  summary_line << analysis_name << ": ";
631 
632  for (const std::pair<str,double>& pair_j : analysis_loglikes.sr_loglikes)
633  {
634  const str& sr_label = pair_j.first;
635  const double& sr_loglike = pair_j.second;
636  const int sr_index = analysis_loglikes.sr_indices.at(sr_label);
637 
638  const str key = analysis_name + "__" + sr_label + "__i" + std::to_string(sr_index) + "__LogLike";
639  result[key] = sr_loglike;
640 
641  summary_line << sr_label + "__i" + std::to_string(sr_index) << ":" << sr_loglike << ", ";
642  }
643 
644  result[analysis_name + "__combined_LogLike"] = analysis_loglikes.combination_loglike;
645 
646  summary_line << "combined_LogLike:" << analysis_loglikes.combination_loglike << ", ";
647  }
648  logger() << LogTags::debug << summary_line.str() << EOM;
649  }
650 
651 
654  {
655  using namespace Pipes::get_LHC_LogLike_per_SR;
656  for (const std::pair<str,AnalysisLogLikes>& pair_i : *Dep::LHC_LogLikes)
657  {
658  const str& analysis_name = pair_i.first;
659  const AnalysisLogLikes& analysis_loglikes = pair_i.second;
660 
661  result[analysis_name] = analysis_loglikes.combination_sr_label;
662  }
663  }
664 
665 
669  {
670  using namespace Pipes::get_LHC_LogLike_per_SR;
671 
672  std::stringstream summary_line;
673  summary_line << "LHC loglike SR indices: ";
674 
675  // Loop over analyses
676  for (const std::pair<str,AnalysisLogLikes>& pair_i : *Dep::LHC_LogLikes)
677  {
678  const str& analysis_name = pair_i.first;
679  const AnalysisLogLikes& analysis_loglikes = pair_i.second;
680 
681  result[analysis_name] = (double) analysis_loglikes.combination_sr_index;
682 
683  summary_line << analysis_name << ":" << analysis_loglikes.combination_sr_index << ", ";
684  }
685  logger() << LogTags::debug << summary_line.str() << EOM;
686  }
687 
688 
690  void calc_combined_LHC_LogLike(double& result)
691  {
692  using namespace Pipes::calc_combined_LHC_LogLike;
693  result = 0.0;
694 
695  // Read analysis names from the yaml file
696  std::vector<str> default_skip_analyses; // The default is empty lists of analyses to skip
697  static const std::vector<str> skip_analyses = runOptions->getValueOrDef<std::vector<str> >(default_skip_analyses, "skip_analyses");
698 
699  // If too many events have failed, do the conservative thing and return delta log-likelihood = 0
700  if (Dep::RunMC->exceeded_maxFailedEvents)
701  {
702  #ifdef COLLIDERBIT_DEBUG
703  cout << debug_prefix() << "calc_combined_LHC_LogLike: Too many failed events. Will be conservative and return a delta log-likelihood of 0." << endl;
704  #endif
705  return;
706  }
707 
708  // Loop over analyses and calculate the total observed dLL
709  for (auto const& analysis_loglike_pair : *Dep::LHC_LogLike_per_analysis)
710  {
711  const str& analysis_name = analysis_loglike_pair.first;
712  const double& analysis_loglike = analysis_loglike_pair.second;
713 
714  // If the analysis name is in skip_analyses, don't add its loglike to the total loglike.
715  if (std::find(skip_analyses.begin(), skip_analyses.end(), analysis_name) != skip_analyses.end())
716  {
717  #ifdef COLLIDERBIT_DEBUG
718  cout.precision(5);
719  cout << debug_prefix() << "calc_combined_LHC_LogLike: Leaving out analysis " << analysis_name << " with LogL = " << analysis_loglike << endl;
720  #endif
721  continue;
722  }
723 
724  // Add analysis loglike
725  result += analysis_loglike;
726 
727  #ifdef COLLIDERBIT_DEBUG
728  cout.precision(5);
729  cout << debug_prefix() << "calc_combined_LHC_LogLike: Analysis " << analysis_name << " contributes with a LogL = " << analysis_loglike << endl;
730  #endif
731  }
732 
733  #ifdef COLLIDERBIT_DEBUG
734  cout << debug_prefix() << "calc_combined_LHC_LogLike: LHC_Combined_LogLike = " << result << endl;
735  #endif
736 
737  // If using capped likelihood, set result = min(result,0)
738  static const bool use_cap_loglike = runOptions->getValueOrDef<bool>(false, "cap_loglike");
739  if (use_cap_loglike)
740  {
741  result = std::min(result, 0.0);
742 
743  #ifdef COLLIDERBIT_DEBUG
744  cout << debug_prefix() << "calc_combined_LHC_LogLike: LHC_Combined_LogLike (capped) = " << result << endl;
745  #endif
746  }
747 
748  std::stringstream summary_line;
749  summary_line << "LHC combined loglike:" << result;
750  logger() << LogTags::debug << summary_line.str() << EOM;
751  }
752 
753  }
754 
755 }
void get_LHC_LogLike_per_SR(map_str_dbl &result)
Extract the log likelihood for each SR.
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
Rollcall header for ColliderBit module.
void calc_combined_LHC_LogLike(double &result)
Compute the total likelihood combining all analyses.
void get_LHC_LogLike_SR_indices(map_str_dbl &result)
Extract the indices for the SRs used in the analysis loglikes.
std::map< std::string, AnalysisLogLikes > map_str_AnalysisLogLikes
Typedef for a string-to-AnalysisLogLikes map.
std::map< std::string, int > sr_indices
Eigen::MatrixXd srcov
Optional covariance matrix between SRs (0x0 null matrix = no correlation info)
double n_signal_at_lumi
n_signal, scaled to the experimental luminosity
Container for loglike information for an analysis.
double n_signal
The number of simulated model events passing selection for this signal region.
double n_observed
The number of events passing selection for this signal region as reported by the experiment.
void calc_LHC_LogLikes(map_str_AnalysisLogLikes &result)
Loop over all analyses and fill a map of AnalysisLogLikes objects.
void get_LHC_LogLike_per_analysis(map_str_dbl &result)
Extract the combined log likelihood for each analysis.
void get_LHC_LogLike_SR_labels(map_str_str &result)
Extract the labels for the SRs used in the analysis loglikes.
std::map< std::string, double > sr_loglikes
A container for the result of an analysis, potentially with many signal regions and correlations...
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
std::string sr_label
A label for the particular signal region of the analysis.
double n_background
The number of standard model events expected to pass the selection for this signal region...
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
std::string analysis_name
Analysis name.
Header file that includes all GAMBIT headers required for a module source file.
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
A simple container for the result of one signal region from one analysis.
void calc_LHC_signals(map_str_dbl &result)
Loop over all analyses and fill a map of predicted counts.
DS5_MSPCTM DS_INTDOF int
std::map< std::string, std::string > map_str_str
Shorthand for a string-to-string map.
Definition: util_types.hpp:74
double background_sys
The absolute systematic error of n_background.
invalid_point_exception & invalid_point()
Invalid point exceptions.
str debug_prefix()
Debug prefix string giving thread number.
std::map< std::string, double > map_str_dbl
Shorthand for a string-to-double map.
AnalysisDataPointers lnlike_marg_poisson_lognormal_error
static Utils::threadsafe_rng & rng()
Return a threadsafe wrapper for the chosen RNG engine (to be passed to e.g.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
size_t size() const
Number of analyses.
double signal_sys
The absolute systematic error of n_signal.