gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
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 
43 
44 #include "multimin/multimin.hpp"
45 
47 #include "Eigen/Eigenvalues"
49 
50 #include <gsl/gsl_sf_gamma.h>
51 
52 //#define COLLIDERBIT_DEBUG
53 #define DEBUG_PREFIX "DEBUG: OMP thread " << omp_get_thread_num() << ": "
54 
55 namespace Gambit
56 {
57 
58  namespace ColliderBit
59  {
60 
61 
64  {
65  using namespace Pipes::calc_LHC_signals;
66 
67  // Clear the result map
68  result.clear();
69 
70  std::stringstream summary_line;
71  summary_line << "LHC signals per SR: ";
72 
73  // Loop over analyses and collect the predicted events into the map
74  for (size_t analysis = 0; analysis < Dep::AllAnalysisNumbers->size(); ++analysis)
75  {
76  // AnalysisData for this analysis
77  const AnalysisData& adata = *(Dep::AllAnalysisNumbers->at(analysis));
78 
79  summary_line << adata.analysis_name << ": ";
80 
81  // Loop over the signal regions inside the analysis, and save the predicted number of events for each.
82  for (size_t SR = 0; SR < adata.size(); ++SR)
83  {
84  // Save SR numbers and absolute uncertainties
85  const SignalRegionData srData = adata[SR];
86  const str key = adata.analysis_name + "__" + srData.sr_label + "__i" + std::to_string(SR) + "__signal";
87  result[key] = srData.n_sig_scaled;
88  const double n_sig_scaled_err = srData.calc_n_sig_scaled_err();
89  result[key + "_uncert"] = n_sig_scaled_err;
90 
91  summary_line << srData.sr_label + "__i" + std::to_string(SR) << ":" << srData.n_sig_scaled << "+-" << n_sig_scaled_err << ", ";
92  }
93  }
94  logger() << LogTags::debug << summary_line.str() << EOM;
95  }
96 
97 
98 
99 
103  void _gsl_calc_Analysis_MinusLogLike(const size_t n, const double* unit_nuisances_dbl,
104  void* fixedparamspack, double* fval) {
105 
106  // Convert the array of doubles into an "Eigen view" of the nuisance params
107  Eigen::Map<const Eigen::ArrayXd> unit_nuisances(&unit_nuisances_dbl[0], n);
108 
109  // Convert the linearised array of doubles into "Eigen views" of the fixed params
110  double *fixedparamspack_dbl = (double*) fixedparamspack;
111  Eigen::Map<const Eigen::VectorXd> n_preds_nominal(&fixedparamspack_dbl[0], n);
112  Eigen::Map<const Eigen::ArrayXd> n_obss(&fixedparamspack_dbl[n], n);
113  Eigen::Map<const Eigen::ArrayXd> sqrtevals(&fixedparamspack_dbl[2*n], n);
114  Eigen::Map<const Eigen::MatrixXd> evecs(&fixedparamspack_dbl[3*n], n, n);
115 
116  // Rotate rate deltas into the SR basis and shift by SR mean rates
117  const Eigen::VectorXd n_preds = n_preds_nominal + evecs*(sqrtevals*unit_nuisances).matrix();
118 
119  // Calculate each SR's Poisson likelihood and add to composite likelihood calculation
120  double loglike_tot = n * log(1/sqrt(2*M_PI)); //< could also drop this, but it costs ~nothing
121  for (size_t j = 0; j < n; ++j) {
122 
123  // First the multivariate Gaussian bit (j = nuisance)
124  const double pnorm_j = -pow(unit_nuisances(j), 2)/2.;
125  loglike_tot += pnorm_j;
126 
127  // Then the Poisson bit (j = SR)
129  const double lambda_j = std::max(n_preds(j), 1e-3); //< manually avoid <= 0 rates
130  const double logfact_n_obs = 0; // gsl_sf_lngamma(n_obss(j) + 1); //< skipping log(n_obs!) computation
131  const double loglike_j = n_obss(j)*log(lambda_j) - lambda_j - logfact_n_obs;
132 
133  loglike_tot += loglike_j;
134  }
135 
136  // Output via argument (times -1 to return -LL for minimisation)
137  *fval = -loglike_tot;
138  }
139 
140 
142  void _gsl_calc_Analysis_MinusLogLikeGrad(const size_t n, const double* unit_nuisances_dbl,
143  void* fixedparamspack, double* fgrad) {
144 
145  // Convert the array of doubles into an "Eigen view" of the nuisance params
146  Eigen::Map<const Eigen::ArrayXd> unit_nuisances(&unit_nuisances_dbl[0], n);
147 
148  // Convert the linearised array of doubles into "Eigen views" of the fixed params
149  double *fixedparamspack_dbl = (double*) fixedparamspack;
150  Eigen::Map<const Eigen::VectorXd> n_preds_nominal(&fixedparamspack_dbl[0], n);
151  Eigen::Map<const Eigen::ArrayXd> n_obss(&fixedparamspack_dbl[n], n);
152  Eigen::Map<const Eigen::ArrayXd> sqrtevals(&fixedparamspack_dbl[2*n], n);
153  Eigen::Map<const Eigen::MatrixXd> evecs(&fixedparamspack_dbl[3*n], n, n);
154 
155  // Rotate rate deltas into the SR basis and shift by SR mean rates
156  const Eigen::VectorXd n_preds = n_preds_nominal + evecs*(sqrtevals*unit_nuisances).matrix();
157 
158  // Compute gradient elements
159  for (int j = 0; j < unit_nuisances.size(); ++j) {
160  double llgrad = 0;
161  for (int k = 0; k < unit_nuisances.size(); ++k) {
162  llgrad += (n_obss(k)/n_preds(k) - 1) * evecs(k,j);
163  }
164  llgrad = llgrad * sqrtevals(j) - unit_nuisances(j);
165  // Output via argument (times -1 to return -dLL for minimisation)
166  fgrad[j] = -llgrad;
167  }
168  }
169 
170 
171  void _gsl_calc_Analysis_MinusLogLikeAndGrad(const size_t n, const double* unit_nuisances_dbl,
172  void* fixedparamspack,
173  double* fval, double* fgrad) {
174  _gsl_calc_Analysis_MinusLogLike(n, unit_nuisances_dbl, fixedparamspack, fval);
175  _gsl_calc_Analysis_MinusLogLikeGrad(n, unit_nuisances_dbl, fixedparamspack, fgrad);
176  }
177 
178 
179  std::vector<double> _gsl_mkpackedarray(const Eigen::ArrayXd& n_preds,
180  const Eigen::ArrayXd& n_obss,
181  const Eigen::ArrayXd& sqrtevals,
182  const Eigen::MatrixXd& evecs) {
183  const size_t nSR = n_obss.size();
184  std::vector<double> fixeds(3*nSR + 2*nSR*nSR, 0.0);
185  for (size_t i = 0; i < nSR; ++i) {
186  fixeds[0+i] = n_preds(i);
187  fixeds[nSR+i] = n_obss(i);
188  fixeds[2*nSR+i] = sqrtevals(i);
189  for (size_t j = 0; j < nSR; ++j) {
190  fixeds[3*nSR+i*nSR+j] = evecs(j,i);
191  }
192  }
193 
194  return fixeds;
195  }
196 
197 
201  double profile_loglike_cov(const Eigen::ArrayXd& n_preds,
202  const Eigen::ArrayXd& n_obss,
203  const Eigen::ArrayXd& sqrtevals,
204  const Eigen::MatrixXd& evecs) {
205 
206  // Number of signal regions
207  const size_t nSR = n_obss.size();
208 
209  // Set initial guess for nuisances to zero
210  std::vector<double> nuisances(nSR, 0.0);
211 
212  // // Set nuisances to an informed starting position
213  // const Eigen::ArrayXd& err_n_preds = (evecs*sqrtevals.matrix()).array(); //< @todo CHECK
214  // std::vector<double> nuisances(nSR, 0.0);
215  // for (size_t j = 0; j < nSR; ++j) {
216  // // Calculate the max-L starting position, ignoring correlations
217  // const double obs = n_obss(j);
218  // const double rate = n_preds(j);
219  // const double delta = err_n_preds(j);
220  // const double a = delta;
221  // const double b = rate + delta*delta;
222  // const double c = delta * (rate - obs);
223  // const double d = b*b - 4*a*c;
224  // const double sqrtd = (d < 0) ? 0 : sqrt(d);
225  // if (sqrtd == 0) {
226  // nuisances[j] = -b / (2*a);
227  // } else {
228  // const double th0_a = (-b + sqrtd) / (2*a);
229  // const double th0_b = (-b - sqrtd) / (2*a);
230  // nuisances[j] = (fabs(th0_a) < fabs(th0_b)) ? th0_a : th0_b;
231  // }
232  // }
233 
234 
235  // Optimiser parameters
236  // Params: step1size, tol, maxiter, epsabs, simplex maxsize, method, verbosity
237  // Methods:
238  // 0: Fletcher-Reeves conjugate gradient
239  // 1: Polak-Ribiere conjugate gradient
240  // 2: Vector Broyden-Fletcher-Goldfarb-Shanno method
241  // 3: Steepest descent algorithm
242  // 4: Nelder-Mead simplex
243  // 5: Vector Broyden-Fletcher-Goldfarb-Shanno method ver. 2
244  // 6: Simplex algorithm of Nelder and Mead ver. 2
245  // 7: Simplex algorithm of Nelder and Mead: random initialization
246  using namespace Pipes::calc_LHC_LogLikes;
247  static const double INITIAL_STEP = runOptions->getValueOrDef<double>(0.1, "nuisance_prof_initstep");
248  static const double CONV_TOL = runOptions->getValueOrDef<double>(0.01, "nuisance_prof_convtol");
249  static const unsigned MAXSTEPS = runOptions->getValueOrDef<unsigned>(10000, "nuisance_prof_maxsteps");
250  static const double CONV_ACC = runOptions->getValueOrDef<double>(0.01, "nuisance_prof_convacc");
251  static const double SIMPLEX_SIZE = runOptions->getValueOrDef<double>(1e-5, "nuisance_prof_simplexsize");
252  static const unsigned METHOD = runOptions->getValueOrDef<unsigned>(6, "nuisance_prof_method");
253  static const unsigned VERBOSITY = runOptions->getValueOrDef<unsigned>(0, "nuisance_prof_verbosity");
254  static const struct multimin::multimin_params oparams = {INITIAL_STEP, CONV_TOL, MAXSTEPS, CONV_ACC, SIMPLEX_SIZE, METHOD, VERBOSITY};
255 
256  // Convert the linearised array of doubles into "Eigen views" of the fixed params
257  std::vector<double> fixeds = _gsl_mkpackedarray(n_preds, n_obss, sqrtevals, evecs);
258 
259  // Pass to the minimiser
260  double minusbestll = 999;
261  // _gsl_calc_Analysis_MinusLogLike(nSR, &nuisances[0], &fixeds[0], &minusbestll);
262  multimin::multimin(nSR, &nuisances[0], &minusbestll,
263  nullptr, nullptr, nullptr,
267  &fixeds[0], oparams);
268 
269  return -minusbestll;
270  }
271 
272 
273  double marg_loglike_nulike1sr(const Eigen::ArrayXd& n_preds,
274  const Eigen::ArrayXd& n_obss,
275  const Eigen::ArrayXd& sqrtevals) {
276  assert(n_preds.size() == 1);
277  assert(n_obss.size() == 1);
278  assert(sqrtevals.size() == 1);
279 
280  using namespace Pipes::calc_LHC_LogLikes;
281  auto marginaliser = (*BEgroup::lnlike_marg_poisson == "lnlike_marg_poisson_lognormal_error")
282  ? BEreq::lnlike_marg_poisson_lognormal_error : BEreq::lnlike_marg_poisson_gaussian_error;
283 
284  // Setting bkg above zero to avoid nulike special cases
285  const double sr_margll = marginaliser((int) n_obss(0), 0.001, n_preds(0), sqrtevals(0)/n_preds(0));
286  return sr_margll;
287  }
288 
289 
290  double marg_loglike_cov(const Eigen::ArrayXd& n_preds,
291  const Eigen::ArrayXd& n_obss,
292  const Eigen::ArrayXd& sqrtevals,
293  const Eigen::MatrixXd& evecs) {
294 
295  // Number of signal regions
296  const size_t nSR = n_obss.size();
297 
298  // Sample correlated SR rates from a rotated Gaussian defined by the covariance matrix and offset by the mean rates
299  using namespace Pipes::calc_LHC_LogLikes;
300  static const double CONVERGENCE_TOLERANCE_ABS = runOptions->getValueOrDef<double>(0.05, "nuisance_marg_convthres_abs");
301  static const double CONVERGENCE_TOLERANCE_REL = runOptions->getValueOrDef<double>(0.05, "nuisance_marg_convthres_rel");
302  static const size_t NSAMPLE_INPUT = runOptions->getValueOrDef<size_t>(100000, "nuisance_marg_nsamples_start");
303  static const bool NULIKE1SR = runOptions->getValueOrDef<bool>(false, "nuisance_marg_nulike1sr");
304 
305  // Optionally use nulike's more careful 1D marginalisation for one-SR cases
306  if (NULIKE1SR && nSR == 1) return marg_loglike_nulike1sr(n_preds, n_obss, sqrtevals);
307 
308  // Dynamic convergence control & test variables
309  size_t nsample = NSAMPLE_INPUT;
310  bool first_iteration = true;
311  double diff_abs = 9999;
312  double diff_rel = 1;
313 
314  // Likelihood variables (note use of long double to guard against blow-up of L as opposed to log(L1/L0))
315  long double ana_like_prev = 1;
316  long double ana_like = 1;
317  long double lsum_prev = 0;
318 
319  // Sampler for unit-normal nuisances
320  std::normal_distribution<double> unitnormdbn(0,1);
321 
322  // Log factorial of observed number of events.
323  // Currently use the ln(Gamma(x)) function gsl_sf_lngamma from GSL. (Need continuous function.)
324  // We may want to switch to using Stirling's approximation: ln(n!) ~ n*ln(n) - n
325  Eigen::ArrayXd logfact_n_obss(nSR);
326  for (size_t j = 0; j < nSR; ++j)
327  logfact_n_obss(j) = gsl_sf_lngamma(n_obss(j) + 1);
328 
329  // Check absolute difference between independent estimates
331  while ((diff_abs > CONVERGENCE_TOLERANCE_ABS && diff_rel > CONVERGENCE_TOLERANCE_REL) || 1.0/sqrt(nsample) > CONVERGENCE_TOLERANCE_ABS)
332  {
333  long double lsum = 0;
334 
341 
342  #pragma omp parallel
343  {
344  // Sample correlated SR rates from a rotated Gaussian defined by the covariance matrix and offset by the mean rates
345  double lsum_private = 0;
346  #pragma omp for nowait
347  for (size_t i = 0; i < nsample; ++i) {
348 
349  Eigen::VectorXd norm_samples(nSR);
350  for (size_t j = 0; j < nSR; ++j)
351  norm_samples(j) = sqrtevals(j) * unitnormdbn(Random::rng());
352 
353  // Rotate rate deltas into the SR basis and shift by SR mean rates
354  const Eigen::VectorXd n_pred_samples = n_preds + (evecs*norm_samples).array();
355 
356  // Calculate Poisson likelihood and add to composite likelihood calculation
357  double combined_loglike = 0;
358  for (size_t j = 0; j < nSR; ++j) {
359  const double lambda_j = std::max(n_pred_samples(j), 1e-3); //< manually avoid <= 0 rates
360  const double loglike_j = n_obss(j)*log(lambda_j) - lambda_j - logfact_n_obss(j);
361  combined_loglike += loglike_j;
362  }
363  // Add combined likelihood to running sums (to later calculate averages)
364  lsum_private += exp(combined_loglike);
365  }
366 
367  #pragma omp critical
368  {
369  lsum += lsum_private;
370  }
371  } // End omp parallel
372 
373  // Compare convergence to previous independent batch
374  if (first_iteration) // The first round must be generated twice
375  {
376  lsum_prev = lsum;
377  first_iteration = false;
378  }
379  else
380  {
381  ana_like_prev = lsum_prev / (double)nsample;
382  ana_like = lsum / (double)nsample;
383  diff_abs = fabs(ana_like_prev - ana_like);
384  diff_rel = diff_abs/ana_like;
385 
386  // Update variables
387  lsum_prev += lsum; // Aggregate result. This doubles the effective batch size for lsum_prev.
388  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.
389  }
390 
391  #ifdef COLLIDERBIT_DEBUG
392  cout << DEBUG_PREFIX
393  << "diff_rel: " << diff_rel << endl
394  << " diff_abs: " << diff_abs << endl
395  << " logl: " << log(ana_like) << endl;
396  cout << DEBUG_PREFIX << "nsample for the next iteration is: " << nsample << endl;
397  cout << DEBUG_PREFIX << endl;
398  #endif
399  }
400  // End convergence while-loop
401 
402  // Combine the independent estimates ana_like and ana_like_prev.
403  // Use equal weights since the estimates are based on equal batch sizes.
404  ana_like = 0.5*(ana_like + ana_like_prev);
405  const double ana_margll = log(ana_like);
406  #ifdef COLLIDERBIT_DEBUG
407  cout << DEBUG_PREFIX << "Combined estimate: ana_loglike: " << ana_margll << " (based on 2*nsample=" << 2*nsample << " samples)" << endl;
408  #endif
409 
410  return ana_margll;
411  }
412 
413 
416  {
417  // Read options
418  using namespace Pipes::calc_LHC_LogLikes;
419  // Use covariance matrices if available?
420  static const bool USE_COVAR = runOptions->getValueOrDef<bool>(true, "use_covariances");
421  // Use marginalisation rather than profiling (probably less stable)?
422  static const bool USE_MARG = runOptions->getValueOrDef<bool>(false, "use_marginalising");
423 
424  // Fix the profiling/marginalising function according to the option
425  auto marg_prof_fn = USE_MARG ? marg_loglike_cov : profile_loglike_cov;
426 
427 
428  // Clear the result map
429  result.clear();
430 
431  // Main loop over all analyses to compute DLL = LL_sb - LL_b
432  for (size_t analysis = 0; analysis < Dep::AllAnalysisNumbers->size(); ++analysis)
433  {
434 
435 
436  // AnalysisData for this analysis
437  const AnalysisData& adata = *(Dep::AllAnalysisNumbers->at(analysis));
438  const std::string ananame = adata.analysis_name;
439  const size_t nSR = adata.size();
440  const bool has_covar = adata.srcov.rows() > 0;
441 
442  #ifdef COLLIDERBIT_DEBUG
443  std::streamsize stream_precision = cout.precision(); // get current precision
444  cout.precision(2); // set precision
445  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: " << "Will print content of " << ananame << " signal regions:" << endl;
446  for (size_t SR = 0; SR < adata.size(); ++SR)
447  {
448  const SignalRegionData& srData = adata[SR];
449  cout << std::fixed << DEBUG_PREFIX
450  << "calc_LHC_LogLikes: " << ananame
451  << ", " << srData.sr_label
452  << ", n_b = " << srData.n_bkg << " +/- " << srData.n_bkg_err
453  << ", n_obs = " << srData.n_obs
454  << ", excess = " << srData.n_obs - srData.n_bkg << " +/- " << srData.n_bkg_err
455  << ", n_s = " << srData.n_sig_scaled
456  << ", (excess-n_s) = " << (srData.n_obs-srData.n_bkg) - srData.n_sig_scaled << " +/- " << srData.n_bkg_err
457  << ", n_s_MC = " << srData.n_sig_MC
458  << endl;
459  }
460  cout.precision(stream_precision); // restore previous precision
461  #endif
462 
463 
464  // Shortcut #1
465  //
466  // If no events have been generated (xsec veto) or too many events have
467  // failed, short-circut the loop and return delta log-likelihood = 0 for
468  // every SR in each analysis.
469  //
471  if (not Dep::RunMC->event_generation_began || Dep::RunMC->exceeded_maxFailedEvents)
472  {
473  // If this is an analysis with covariance info, only add a single 0-entry in the map
474  if (USE_COVAR && has_covar)
475  {
476  result[ananame].combination_sr_label = "none";
477  result[ananame].combination_sr_index = -1;
478  result[ananame].combination_loglike = 0.0;
479  }
480  // If this is an analysis without covariance info, add 0-entries for all SRs plus
481  // one for the combined LogLike
482  else
483  {
484  for (size_t SR = 0; SR < adata.size(); ++SR)
485  {
486  result[ananame].sr_indices[adata[SR].sr_label] = SR;
487  result[ananame].sr_loglikes[adata[SR].sr_label] = 0.0;
488  }
489  result[ananame].combination_sr_label = "none";
490  result[ananame].combination_sr_index = -1;
491  result[ananame].combination_loglike = 0.0;
492  }
493 
494  #ifdef COLLIDERBIT_DEBUG
495  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: " << ananame << "_LogLike : " << 0.0 << " (No events predicted / successfully generated. Skipped full calculation.)" << endl;
496  #endif
497 
498  // Continue to next analysis
499  continue;
500  }
501 
502 
503  // Shortcut #2
504  //
505  // If all SRs have 0 signal prediction, we know the delta log-likelihood is 0.
506  bool all_zero_signal = true;
507  for (size_t SR = 0; SR < nSR; ++SR)
508  {
509  if (adata[SR].n_sig_MC != 0)
510  {
511  all_zero_signal = false;
512  break;
513  }
514  }
515  if (all_zero_signal)
516  {
517  // Store result
518  if (!(USE_COVAR && has_covar))
519  {
520  for (size_t SR = 0; SR < adata.size(); ++SR)
521  {
522  result[ananame].sr_indices[adata[SR].sr_label] = SR;
523  result[ananame].sr_loglikes[adata[SR].sr_label] = 0.0;
524  }
525  }
526  result[ananame].combination_sr_label = "any";
527  result[ananame].combination_sr_index = -1;
528  result[ananame].combination_loglike = 0.0;
529 
530  #ifdef COLLIDERBIT_DEBUG
531  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: " << ananame << "_LogLike : " << 0.0 << " (No signal predicted. Skipped full calculation.)" << endl;
532  #endif
533 
534  // Continue to next analysis
535  continue;
536  }
537 
538 
539  // Work out the total (delta) log likelihood for this analysis, with correlations as available/instructed
540  double ana_dll = NAN;
541  if (USE_COVAR && has_covar)
542  {
543 
544 
557  #ifdef COLLIDERBIT_DEBUG
558  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: Analysis " << analysis << " has a covariance matrix: computing composite loglike." << endl;
559  #endif
560 
561 
562  // Construct vectors of SR numbers
564  Eigen::ArrayXd n_obs(nSR); // logfact_n_obs(nSR);
565  Eigen::ArrayXd n_pred_b(nSR), n_pred_sb(nSR), abs_unc_s(nSR);
566  for (size_t SR = 0; SR < nSR; ++SR)
567  {
568  const SignalRegionData& srData = adata[SR];
569 
570  // Actual observed number of events
571  n_obs(SR) = srData.n_obs;
572 
573  // Log factorial of observed number of events.
574  // Currently use the ln(Gamma(x)) function gsl_sf_lngamma from GSL. (Need continuous function.)
575  // We may want to switch to using Stirling's approximation: ln(n!) ~ n*ln(n) - n
576  //logfact_n_obs(SR) = gsl_sf_lngamma(n_obs(SR) + 1.);
577 
578  // A contribution to the predicted number of events that is not known exactly
579  n_pred_b(SR) = std::max(srData.n_bkg, 0.001); // <-- Avoid trouble with b==0
580  n_pred_sb(SR) = srData.n_sig_scaled + srData.n_bkg;
581 
582  // Absolute errors for n_predicted_uncertain_*
583  abs_unc_s(SR) = srData.calc_n_sig_scaled_err();
584  }
585 
586  // Diagonalise the background-only covariance matrix, extracting the correlation and rotation matrices
588  const Eigen::MatrixXd& srcov_b = adata.srcov;
589  Eigen::MatrixXd srcorr_b = srcov_b; // start with cov, then make corr
590  for (size_t SR = 0; SR < nSR; ++SR)
591  {
592  const double diagsd = sqrt(srcov_b(SR,SR));
593  srcorr_b.row(SR) /= diagsd;
594  srcorr_b.col(SR) /= diagsd;
595  }
596  const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig_b(adata.srcov);
597  const Eigen::ArrayXd Eb = eig_b.eigenvalues();
598  const Eigen::ArrayXd sqrtEb = Eb.sqrt();
599  const Eigen::MatrixXd Vb = eig_b.eigenvectors();
600 
601  // Construct and diagonalise the s+b covariance matrix, adding the diagonal signal uncertainties in quadrature
602  const Eigen::MatrixXd srcov_s = abs_unc_s.array().square().matrix().asDiagonal();
603  const Eigen::MatrixXd srcov_sb = srcov_b + srcov_s;
604  Eigen::MatrixXd srcorr_sb = srcov_sb;
605  for (size_t SR = 0; SR < nSR; ++SR)
606  {
607  const double diagsd = sqrt(srcov_sb(SR,SR));
608  srcorr_sb.row(SR) /= diagsd;
609  srcorr_sb.col(SR) /= diagsd;
610  }
611  const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig_sb(srcov_sb);
612  const Eigen::ArrayXd Esb = eig_sb.eigenvalues();
613  const Eigen::ArrayXd sqrtEsb = Esb.sqrt();
614  const Eigen::MatrixXd Vsb = eig_sb.eigenvectors();
615 
616  // cout << "B: " << srcorr_b << " " << srcov_b << endl;
617  // cout << "SB: " << srcorr_sb << " " << srcov_sb << endl;
618 
619  // Compute the single, correlated analysis-level DLL as the difference of s+b and b (partial) LLs
621  const double ll_b = marg_prof_fn(n_pred_b, n_obs, sqrtEb, Vb);
622  const double ll_sb = marg_prof_fn(n_pred_sb, n_obs, sqrtEsb, Vsb);
623  const double dll = ll_sb - ll_b;
624 
625  // Store result
626  ana_dll = dll;
627  result[ananame].combination_sr_label = "all";
628  result[ananame].combination_sr_index = -1;
629  result[ananame].combination_loglike = ana_dll;
630 
631  #ifdef COLLIDERBIT_DEBUG
632  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: " << ananame << "_LogLike : " << ana_dll << endl;
633  #endif
634 
635 
636  } else { // NO SR-CORRELATION INFO, OR USER CHOSE NOT TO USE IT:
637 
638 
639  // We either take the result from the SR *expected* to be most
640  // constraining under the s=0 assumption (default), or naively combine
641  // the loglikes for all SRs (if combine_SRs_without_covariances=true).
642  #ifdef COLLIDERBIT_DEBUG
643  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: Analysis " << analysis << " has no covariance matrix: computing single best-expected loglike." << endl;
644  #endif
645 
646  double bestexp_dll_exp = 0, bestexp_dll_obs = NAN;
647  str bestexp_sr_label;
648  int bestexp_sr_index;
649  double nocovar_srsum_dll_obs = 0;
650 
651  for (size_t SR = 0; SR < nSR; ++SR)
652  {
653  const SignalRegionData& srData = adata[SR];
654 
655  // Shortcut: If n_sig_MC == 0, we know the delta log-likelihood is 0.
656  if(srData.n_sig_MC == 0)
657  {
658  // Store (obs) result for this SR
659  result[ananame].sr_indices[srData.sr_label] = SR;
660  result[ananame].sr_loglikes[srData.sr_label] = 0.0;
661 
662  // Update the running best-expected-exclusion detail
663  if (0.0 < bestexp_dll_exp || SR == 0)
664  {
665  bestexp_dll_exp = 0.0;
666  bestexp_dll_obs = 0.0;
667  bestexp_sr_label = srData.sr_label;
668  bestexp_sr_index = SR;
669  }
670 
671  // Skip to next SR
672  continue;
673  }
674 
675  // A contribution to the predicted number of events that is not known exactly
676  const double n_pred_b = std::max(srData.n_bkg, 0.001); // <-- Avoid trouble with b==0
677  const double n_pred_sb = n_pred_b + srData.n_sig_scaled;
678 
679  // Actual observed number of events and predicted background, as integers cf. Poisson stats
680  const double n_obs = round(srData.n_obs);
681  const double n_pred_b_int = round(n_pred_b);
682 
683  // Absolute errors for n_predicted_uncertain_*
684  const double abs_uncertainty_b = std::max(srData.n_bkg_err, 0.001); // <-- Avoid trouble with b_err==0
685  const double abs_uncertainty_sb = std::max(srData.calc_n_sigbkg_err(), 0.001); // <-- Avoid trouble with sb_err==0
686 
687 
688  // Construct dummy 1-element Eigen objects for passing to the general likelihood calculator
690  Eigen::ArrayXd n_obss(1); n_obss(0) = n_obs;
691  Eigen::ArrayXd n_preds_b_int(1); n_preds_b_int(0) = n_pred_b_int;
692  Eigen::ArrayXd n_preds_b(1); n_preds_b(0) = n_pred_b;
693  Eigen::ArrayXd n_preds_sb(1); n_preds_sb(0) = n_pred_sb;
694  Eigen::ArrayXd sqrtevals_b(1); sqrtevals_b(0) = abs_uncertainty_b;
695  Eigen::ArrayXd sqrtevals_sb(1); sqrtevals_sb(0) = abs_uncertainty_sb;
696  Eigen::MatrixXd dummy(1,1); dummy(0,0) = 1.0;
697 
698 
699  // Compute this SR's DLLs as the differences of s+b and b (partial) LLs
702  const double ll_b_exp = marg_prof_fn(n_preds_b, n_preds_b_int, sqrtevals_b, dummy);
704  const double ll_b_obs = marg_prof_fn(n_preds_b, n_obss, sqrtevals_b, dummy);
705  const double ll_sb_exp = marg_prof_fn(n_preds_sb, n_preds_b_int, sqrtevals_sb, dummy);
706  const double ll_sb_obs = marg_prof_fn(n_preds_sb, n_obss, sqrtevals_sb, dummy);
707  const double dll_exp = ll_sb_exp - ll_b_exp;
708  const double dll_obs = ll_sb_obs - ll_b_obs;
709 
710  // Check for problems
711  if (Utils::isnan(ll_b_exp))
712  {
713  std::stringstream msg;
714  msg << "Computation of ll_b_exp for signal region " << srData.sr_label << " in analysis " << ananame << " returned NaN" << endl;
715  invalid_point().raise(msg.str());
716  }
717  if (Utils::isnan(ll_b_obs))
718  {
719  std::stringstream msg;
720  msg << "Computation of ll_b_obs for signal region " << srData.sr_label << " in analysis " << ananame << " returned NaN" << endl;
721  invalid_point().raise(msg.str());
722  }
723  if (Utils::isnan(ll_sb_exp))
724  {
725  std::stringstream msg;
726  msg << "Computation of ll_sb_exp for signal region " << srData.sr_label << " in analysis " << ananame << " returned NaN" << endl;
727  invalid_point().raise(msg.str());
728  }
729  if (Utils::isnan(ll_sb_obs))
730  {
731  std::stringstream msg;
732  msg << "Computation of ll_sb_obs for signal region " << srData.sr_label << " in analysis " << ananame << " returned NaN" << endl;
733  invalid_point().raise(msg.str());
734  }
735 
736  // Update the running best-expected-exclusion detail
737  if (dll_exp < bestexp_dll_exp || SR == 0)
738  {
739  bestexp_dll_exp = dll_exp;
740  bestexp_dll_obs = dll_obs;
741  bestexp_sr_label = srData.sr_label;
742  bestexp_sr_index = SR;
743  // #ifdef COLLIDERBIT_DEBUG
744  // cout << DEBUG_PREFIX << "Setting bestexp_sr_label to: " << bestexp_sr_label << ", LogL_exp = " << bestexp_dll_exp << ", LogL_obs = " << bestexp_dll_obs << endl;
745  // #endif
746  }
747 
748  // Store (obs) result for this SR
749  result[ananame].sr_indices[srData.sr_label] = SR;
750  result[ananame].sr_loglikes[srData.sr_label] = dll_obs;
751  // Also add the obs loglike to the no-correlations sum over SRs
752  nocovar_srsum_dll_obs += dll_obs;
753 
754  #ifdef COLLIDERBIT_DEBUG
755  cout << DEBUG_PREFIX << ananame << ", " << srData.sr_label << ", llsb_exp-llb_exp = " << dll_exp << ", llsb_obs-llb_obs= " << dll_obs << endl;
756  #endif
757 
758  }
759 
760  // Set this analysis' total obs DLL to that from the best-expected SR
761  ana_dll = bestexp_dll_obs;
762  result[ananame].combination_sr_label = bestexp_sr_label;
763  result[ananame].combination_sr_index = bestexp_sr_index;
764  result[ananame].combination_loglike = ana_dll;
765 
766  // Or should we use the naive sum of SR loglikes (without correlations) instead?
767  static const bool combine_nocovar_SRs = runOptions->getValueOrDef<bool>(false, "combine_SRs_without_covariances");
768  if (combine_nocovar_SRs)
769  {
770  result[ananame].combination_loglike = nocovar_srsum_dll_obs;
771  }
772 
773  #ifdef COLLIDERBIT_DEBUG
774  cout << DEBUG_PREFIX << "calc_LHC_LogLikes: " << ananame << "_" << bestexp_sr_label << "_LogLike : " << ana_dll << endl;
775  #endif
776 
777  } // end cov/no-cov
778 
779 
780  // Check for problems with the result
781  for(auto& s_d_pair : result[ananame].sr_loglikes)
782  {
783  if (Utils::isnan(s_d_pair.second))
784  {
785  std::stringstream msg;
786  msg << "Computation of loglike for signal region " << s_d_pair.first << " in analysis " << ananame << " returned NaN" << endl;
787  msg << "Will now print the signal region data for this analysis:" << endl;
788  for (size_t SR = 0; SR < nSR; ++SR)
789  {
790  const SignalRegionData& srData = adata[SR];
791  msg << srData.sr_label
792  << ", n_bkg = " << srData.n_bkg
793  << ", n_bkg_err = " << srData.n_bkg_err
794  << ", n_obs = " << srData.n_obs
795  << ", n_sig_scaled = " << srData.n_sig_scaled
796  << ", n_sig_MC = " << srData.n_sig_MC
797  << ", n_sig_MC_sys = " << srData.n_sig_MC_sys
798  << endl;
799  }
800  invalid_point().raise(msg.str());
801  }
802  }
803 
804  } // end analysis loop
805  }
806 
807 
810  {
811  using namespace Pipes::get_LHC_LogLike_per_analysis;
812 
813  std::stringstream summary_line;
814  summary_line << "LHC loglikes per analysis: ";
815 
816  for (const std::pair<str,AnalysisLogLikes>& pair : *Dep::LHC_LogLikes)
817  {
818  const str& analysis_name = pair.first;
819  const AnalysisLogLikes& analysis_loglikes = pair.second;
820 
821  result[analysis_name] = analysis_loglikes.combination_loglike;
822 
823  summary_line << analysis_name << ":" << analysis_loglikes.combination_loglike << ", ";
824  }
825  logger() << LogTags::debug << summary_line.str() << EOM;
826  }
827 
828 
831  {
832  using namespace Pipes::get_LHC_LogLike_per_SR;
833 
834  std::stringstream summary_line;
835  summary_line << "LHC loglikes per SR: ";
836 
837  for (const std::pair<str,AnalysisLogLikes>& pair_i : *Dep::LHC_LogLikes)
838  {
839  const str& analysis_name = pair_i.first;
840  const AnalysisLogLikes& analysis_loglikes = pair_i.second;
841 
842  summary_line << analysis_name << ": ";
843 
844  for (const std::pair<str,double>& pair_j : analysis_loglikes.sr_loglikes)
845  {
846  const str& sr_label = pair_j.first;
847  const double& sr_loglike = pair_j.second;
848  const int sr_index = analysis_loglikes.sr_indices.at(sr_label);
849 
850  const str key = analysis_name + "__" + sr_label + "__i" + std::to_string(sr_index) + "__LogLike";
851  result[key] = sr_loglike;
852 
853  summary_line << sr_label + "__i" + std::to_string(sr_index) << ":" << sr_loglike << ", ";
854  }
855 
856  result[analysis_name + "__combined_LogLike"] = analysis_loglikes.combination_loglike;
857 
858  summary_line << "combined_LogLike:" << analysis_loglikes.combination_loglike << ", ";
859  }
860  logger() << LogTags::debug << summary_line.str() << EOM;
861  }
862 
863 
866  {
867  using namespace Pipes::get_LHC_LogLike_per_SR;
868  for (const std::pair<str,AnalysisLogLikes>& pair_i : *Dep::LHC_LogLikes)
869  {
870  const str& analysis_name = pair_i.first;
871  const AnalysisLogLikes& analysis_loglikes = pair_i.second;
872 
873  result[analysis_name] = analysis_loglikes.combination_sr_label;
874  }
875  }
876 
877 
881  {
882  using namespace Pipes::get_LHC_LogLike_per_SR;
883 
884  std::stringstream summary_line;
885  summary_line << "LHC loglike SR indices: ";
886 
887  // Loop over analyses
888  for (const std::pair<str,AnalysisLogLikes>& pair_i : *Dep::LHC_LogLikes)
889  {
890  const str& analysis_name = pair_i.first;
891  const AnalysisLogLikes& analysis_loglikes = pair_i.second;
892 
893  result[analysis_name] = (double) analysis_loglikes.combination_sr_index;
894 
895  summary_line << analysis_name << ":" << analysis_loglikes.combination_sr_index << ", ";
896  }
897  logger() << LogTags::debug << summary_line.str() << EOM;
898  }
899 
900 
902  void calc_combined_LHC_LogLike(double& result)
903  {
904  using namespace Pipes::calc_combined_LHC_LogLike;
905  result = 0.0;
906 
907  static const bool write_summary_to_log = runOptions->getValueOrDef<bool>(false, "write_summary_to_log");
908 
909  std::stringstream summary_line_combined_loglike;
910  summary_line_combined_loglike << "calc_combined_LHC_LogLike: combined LogLike: ";
911  std::stringstream summary_line_skipped_analyses;
912  summary_line_skipped_analyses << "calc_combined_LHC_LogLike: skipped analyses: ";
913  std::stringstream summary_line_included_analyses;
914  summary_line_included_analyses << "calc_combined_LHC_LogLike: included analyses: ";
915 
916  // Read analysis names from the yaml file
917  std::vector<str> default_skip_analyses; // The default is empty lists of analyses to skip
918  static const std::vector<str> skip_analyses = runOptions->getValueOrDef<std::vector<str> >(default_skip_analyses, "skip_analyses");
919 
920  // If too many events have failed, do the conservative thing and return delta log-likelihood = 0
921  if (Dep::RunMC->exceeded_maxFailedEvents)
922  {
923  #ifdef COLLIDERBIT_DEBUG
924  cout << DEBUG_PREFIX << "calc_combined_LHC_LogLike: Too many failed events. Will be conservative and return a delta log-likelihood of 0." << endl;
925  #endif
926  return;
927  }
928 
929  // Loop over analyses and calculate the total observed dLL
930  for (auto const& analysis_loglike_pair : *Dep::LHC_LogLike_per_analysis)
931  {
932  const str& analysis_name = analysis_loglike_pair.first;
933  const double& analysis_loglike = analysis_loglike_pair.second;
934 
935  // If the analysis name is in skip_analyses, don't add its loglike to the total loglike.
936  if (std::find(skip_analyses.begin(), skip_analyses.end(), analysis_name) != skip_analyses.end())
937  {
938  #ifdef COLLIDERBIT_DEBUG
939  cout.precision(5);
940  cout << DEBUG_PREFIX << "calc_combined_LHC_LogLike: Leaving out analysis " << analysis_name << " with LogL = " << analysis_loglike << endl;
941  #endif
942 
943  // Add to log summary
944  if(write_summary_to_log)
945  {
946  summary_line_skipped_analyses << analysis_name << "__LogLike:" << analysis_loglike << ", ";
947  }
948 
949  continue;
950  }
951 
952  // Add analysis loglike.
953  // If using capped likelihood for each individual analysis, set analysis_loglike = min(analysis_loglike,0)
954  static const bool use_cap_loglike_individual = runOptions->getValueOrDef<bool>(false, "cap_loglike_individual_analyses");
955  if (use_cap_loglike_individual)
956  {
957  result += std::min(analysis_loglike, 0.0);
958  }
959  else
960  {
961  result += analysis_loglike;
962  }
963 
964  // Add to log summary
965  if(write_summary_to_log)
966  {
967  summary_line_included_analyses << analysis_name << "__LogLike:" << analysis_loglike << ", ";
968  }
969 
970  #ifdef COLLIDERBIT_DEBUG
971  cout.precision(5);
972  cout << DEBUG_PREFIX << "calc_combined_LHC_LogLike: Analysis " << analysis_name << " contributes with a LogL = " << analysis_loglike << endl;
973  #endif
974  }
975 
976  #ifdef COLLIDERBIT_DEBUG
977  cout << DEBUG_PREFIX << "calc_combined_LHC_LogLike: LHC_Combined_LogLike = " << result << endl;
978  #endif
979 
980  // If using a "global" capped likelihood, set result = min(result,0)
981  static const bool use_cap_loglike = runOptions->getValueOrDef<bool>(false, "cap_loglike");
982  if (use_cap_loglike)
983  {
984  result = std::min(result, 0.0);
985  }
986 
987  // Write log summary
988  if(write_summary_to_log)
989  {
990  summary_line_combined_loglike << result;
991 
992  logger() << summary_line_combined_loglike.str() << EOM;
993  logger() << summary_line_included_analyses.str() << EOM;
994  logger() << summary_line_skipped_analyses.str() << EOM;
995  }
996  }
997 
998 
1001  void calc_LHC_LogLike_scan_guide(double& result)
1002  {
1003  using namespace Pipes::calc_LHC_LogLike_scan_guide;
1004  result = 0.0;
1005 
1006  static const bool write_summary_to_log = runOptions->getValueOrDef<bool>(false, "write_summary_to_log");
1007  static const double target_LHC_loglike = runOptions->getValue<double>("target_LHC_loglike");
1008  static const double target_width = runOptions->getValue<double>("width_LHC_loglike");
1009 
1010  // Get the combined LHC loglike
1011  double LHC_loglike = *Dep::LHC_Combined_LogLike;
1012 
1013  // Calculate the dummy scan guide loglike using a gaussian centered on the target LHC loglike value
1014  result = Stats::gaussian_loglikelihood(LHC_loglike, target_LHC_loglike, 0.0, target_width, false);
1015 
1016  // Write log summary
1017  if(write_summary_to_log)
1018  {
1019  std::stringstream summary_line;
1020  summary_line << "LHC_LogLike_scan_guide: " << result;
1021  logger() << summary_line.str() << EOM;
1022  }
1023  }
1024 
1025  }
1026 }
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.
double n_sig_MC
The number of simulated model events passing selection for this signal region.
std::map< std::string, int > sr_indices
Eigen::MatrixXd srcov
Optional covariance matrix between SRs (0x0 null matrix = no correlation info)
double profile_loglike_cov(const Eigen::ArrayXd &n_preds, const Eigen::ArrayXd &n_obss, const Eigen::ArrayXd &sqrtevals, const Eigen::MatrixXd &evecs)
Return the best log likelihood.
Container for loglike information for an analysis.
void _gsl_calc_Analysis_MinusLogLikeGrad(const size_t n, const double *unit_nuisances_dbl, void *fixedparamspack, double *fgrad)
Loglike gradient-function wrapper to provide the signature for GSL multimin.
double n_bkg_err
The absolute error of n_bkg.
void _gsl_calc_Analysis_MinusLogLikeAndGrad(const size_t n, const double *unit_nuisances_dbl, void *fixedparamspack, double *fval, double *fgrad)
#define DEBUG_PREFIX
double marg_loglike_cov(const Eigen::ArrayXd &n_preds, const Eigen::ArrayXd &n_obss, const Eigen::ArrayXd &sqrtevals, const Eigen::MatrixXd &evecs)
void calc_LHC_LogLikes(map_str_AnalysisLogLikes &result)
Loop over all analyses and fill a map of AnalysisLogLikes objects.
double gaussian_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is Gaussian...
Definition: statistics.cpp:133
T * matrix(const int xN)
Definition: random_tools.hpp:9
double n_bkg
The number of standard model events expected to pass the selection for this signal region...
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::vector< double > _gsl_mkpackedarray(const Eigen::ArrayXd &n_preds, const Eigen::ArrayXd &n_obss, const Eigen::ArrayXd &sqrtevals, const Eigen::MatrixXd &evecs)
std::map< std::string, double > sr_loglikes
A container for the result of an analysis, potentially with many signal regions and correlations...
double n_obs
The number of events passing selection for this signal region as reported by the experiment.
Declarations of statistical utilities.
virtual void raise(const std::string &)
Raise the exception, i.e. throw it. Exact override of base method.
Definition: exceptions.cpp:422
std::string sr_label
A label for the particular signal region of the analysis.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
std::string analysis_name
Analysis name.
Header file that includes all GAMBIT headers required for a module source file.
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
double n_sig_MC_sys
The absolute systematic error of n_sig_MC.
Simple header file for turning compiler warnings back on after having included one of the begin_ignor...
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.
Pragma directives to suppress compiler warnings coming from including Eigen library headers...
std::map< std::string, std::string > map_str_str
Shorthand for a string-to-string map.
Definition: util_types.hpp:78
double pow(const double &a)
Outputs a^i.
double marg_loglike_nulike1sr(const Eigen::ArrayXd &n_preds, const Eigen::ArrayXd &n_obss, const Eigen::ArrayXd &sqrtevals)
invalid_point_exception & invalid_point()
Invalid point exceptions.
std::map< std::string, double > map_str_dbl
Shorthand for a string-to-double map.
double n_sig_scaled
n_sig_MC, scaled to luminosity * cross-section
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
void _gsl_calc_Analysis_MinusLogLike(const size_t n, const double *unit_nuisances_dbl, void *fixedparamspack, double *fval)
Loglike objective-function wrapper to provide the signature for GSL multimin.
size_t size() const
Number of analyses.
void calc_LHC_LogLike_scan_guide(double &result)
A dummy log-likelihood that helps the scanner track a given range of collider log-likelihood values...