gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
Gambit::ColliderBit Namespace Reference

Namespaces

 ATLAS
 ATLAS-specific efficiency and smearing functions for super fast detector simulation.
 
 CMS
 CMS-specific efficiency and smearing functions for super fast detector simulation.
 

Classes

class  ALEPHSelectronLimitAt208GeV
 A class to contain the limit data from ALEPH_PLB526_2002_206, figure 3a. More...
 
class  ALEPHSmuonLimitAt208GeV
 A class to contain the limit data from ALEPH_PLB526_2002_206, figure 3b. More...
 
class  ALEPHStauLimitAt208GeV
 A class to contain the limit data from ALEPH_PLB526_2002_206, figure 3c. More...
 
class  Analysis
 A class for collider analyses within ColliderBit. More...
 
class  AnalysisContainer
 A class for managing collections of Analysis instances. More...
 
struct  AnalysisData
 A container for the result of an analysis, potentially with many signal regions and correlations. More...
 
struct  AnalysisLogLikes
 Container for loglike information for an analysis. More...
 
class  BaseCollider
 An abstract base class for collider simulators within ColliderBit. More...
 
class  BaseDetector
 An abstract base class for detector simulators within ColliderBit. More...
 
class  BaseLimitContainer
 Base class for experimental limit curve interpolation. More...
 
class  BuckFast
 A base class for BuckFast simple smearing simulations within ColliderBit. More...
 
class  ColliderPythia
 A specializable, recyclable class interfacing ColliderBit and Pythia. More...
 
struct  convergence_settings
 Type for holding Monte Carlo convergence settings. More...
 
struct  Cutflow
 A tracker of numbers & fractions of events passing sequential cuts. More...
 
struct  Cutflows
 A container for several Cutflow objects, with some convenient batch access. More...
 
class  ImageLimit
 
class  L3ChargedGauginoSmallDeltaMAnySneutrinoLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_482_2000_31, figure 5b. More...
 
class  L3ChargedGauginoSmallDeltaMWithHeavySneutrinoLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_482_2000_31, figure 5a. More...
 
class  L3ChargedHiggsinoSmallDeltaMLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_482_2000_31, figure 5c. More...
 
class  L3CharginoAllChannelsLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_472_2000_420, figure 2a. More...
 
class  L3CharginoLeptonicLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_472_2000_420, figure 2b. More...
 
class  L3NeutralinoAllChannelsLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_472_2000_420, figure 3a. More...
 
class  L3NeutralinoLeptonicLimitAt188pt6GeV
 A class to contain the limit data from L3PLB_472_2000_420, figure 3b. More...
 
class  L3SelectronLimitAt205GeV
 A class to contain the limit data from L3_PLB580_2004_37, figure 2a. More...
 
class  L3SmuonLimitAt205GeV
 A class to contain the limit data from L3_PLB580_2004_37, figure 2b. More...
 
class  L3StauLimitAt205GeV
 A class to contain the limit data from L3_PLB580_2004_37, figure 2c. More...
 
class  LineSegment
 A simple container for a line segment on an xy plane. More...
 
class  MC_convergence_checker
 Helper class for testing for convergence of analyses. More...
 
struct  MCLoopInfo
 Container for event loop status data and settings. More...
 
class  OPALCharginoAllChannelsLimitAt208GeV
 A class to contain the limit data from OPAL_EPJC35_2004_1, figure 8b. More...
 
class  OPALCharginoHadronicLimitAt208GeV
 A class to contain the limit data from OPAL_EPJC35_2004_1, figure 5b. More...
 
class  OPALCharginoLeptonicLimitAt208GeV
 A class to contain the limit data from OPAL_EPJC35_2004_1, figure 7b. More...
 
class  OPALCharginoSemiLeptonicLimitAt208GeV
 A class to contain the limit data from OPAL_EPJC35_2004_1, figure 6b. More...
 
class  OPALDegenerateCharginoLimitAt208GeV
 A class to contain the limit data from OPAL, hep-ex/0210043, figure 5a (in colour) More...
 
class  OPALNeutralinoHadronicLimitAt208GeV
 A class to contain the limit data from OPAL_EPJC35_2004_1, figure 9b. More...
 
class  OPALNeutralinoHadronicViaZLimitAt208GeV
 A class to contain the limit data from OPAL_EPJC35_2004_1, figure 9b. More...
 
class  P2
 A simple container for a point on an xy plane. More...
 
class  Perf_Plot
 
struct  SignalRegionData
 A simple container for the result of one signal region from one analysis. More...
 
class  xsec
 A class for holding cross-section info within ColliderBit. More...
 

Typedefs

typedef std::vector< AnalysisDataAnalysisNumbers
 Container for data from multiple analyses and SRs. More...
 
typedef std::vector< const AnalysisData * > AnalysisDataPointers
 
typedef std::map< std::string, AnalysisLogLikesmap_str_AnalysisLogLikes
 Typedef for a string-to-AnalysisLogLikes map. More...
 
typedef std::vector< AnalysisContainerAnalysisContainers
 Container for multiple analysis containers. More...
 
typedef std::chrono::milliseconds ms
 
typedef std::chrono::steady_clock steady_clock
 
typedef std::chrono::steady_clock::time_point tp
 
typedef std::map< std::string, doubletimer_map_type
 
typedef std::vector< std::vector< double > > MixMatrix
 
typedef std::vector< std::vector< double > > data_type
 
typedef std::vector< HEPUtils::Particle * > ParticlePtrs
 Typedef for a vector of Particle pointers. More...
 
typedef std::vector< const HEPUtils::Particle * > ConstParticlePtrs
 Typedef for a vector of const Particle pointers. More...
 
typedef std::vector< HEPUtils::Jet * > JetPtrs
 Typedef for a vector of Jet pointers. More...
 
typedef std::vector< const HEPUtils::Jet * > ConstJetPtrs
 Typedef for a vector of const Jet pointers. More...
 

Enumerations

enum  specialIterations {
  BASE_INIT = -1, COLLIDER_INIT = -2, START_SUBPROCESS = -3, COLLECT_CONVERGENCE_DATA = -4,
  CHECK_CONVERGENCE = -5, END_SUBPROCESS = -6, COLLIDER_FINALIZE = -7, BASE_FINALIZE = -8
}
 Special iteration labels for the loop controlled by operateLHCLoop. More...
 

Functions

AnalysismkAnalysis (const str &name)
 Create a new analysis based on a name string. More...
 
str getDetector (const str &name)
 Return the detector to be used for a given analysis name, checking that the analysis exists. More...
 
ostream & operator<< (ostream &os, const Cutflow &cf)
 Print a Cutflow to a stream. More...
 
ostream & operator<< (ostream &os, const Cutflows &cfs)
 Print a Cutflows to a stream. More...
 
str debug_prefix ()
 Debug prefix string giving thread number. More...
 
template<typename PythiaT , typename EventT >
void generateEventColliderPythia (EventT &result, const MCLoopInfo &RunMC, const ColliderPythia< PythiaT, EventT > &HardScatteringSim, const int iteration, void(*wrapup)())
 Generate a hard scattering event with Pythia. More...
 
template<typename EventT >
BaseDetector< EventT > * getBuckFast (const str &detname, const MCLoopInfo &RunMC, int iteration, const Options &runOptions)
 Get a BuckFast detector simulation. More...
 
template<typename PythiaT , typename EventT >
void getColliderPythia (ColliderPythia< PythiaT, EventT > &result, const MCLoopInfo &RunMC, const Spectrum &spectrum, const DecayTable &decay_rates, const str model_suffix, const int iteration, void(*wrapup)(), const Options &runOptions, bool is_SUSY)
 Retrieve a Pythia hard-scattering Monte Carlo simulation. More...
 
LineSegment makeLine (const P2 &pt1, const P2 &pt2)
 Factory function for lines. More...
 
double addInQuad (const double &a, const double &b)
 Add two numbers in quadrature. More...
 
template<typename EventT >
void smearEvent (HEPUtils::Event &result, const EventT &HardScatteringEvent, const BaseDetector< EventT > &Smearer, const MCLoopInfo &RunMC, const int iteration, const str &ID, const str &detname, void(*wrapup)())
 Smear an event. More...
 
bool sortByPT13 (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByPT13_sharedptr (std::shared_ptr< HEPUtils::Jet > jet1, std::shared_ptr< HEPUtils::Jet > jet2)
 
bool sortByMass (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByMass_sharedptr (std::shared_ptr< HEPUtils::Jet > jet1, std::shared_ptr< HEPUtils::Jet > jet2)
 
double calcMT (HEPUtils::P4 jetMom, HEPUtils::P4 metMom)
 
bool sortByPT_1l (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByPT_1l_sharedptr (std::shared_ptr< HEPUtils::Jet > jet1, std::shared_ptr< HEPUtils::Jet > jet2)
 
bool sortByMass_1l (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByMass_1l_sharedptr (std::shared_ptr< HEPUtils::Jet > jet1, std::shared_ptr< HEPUtils::Jet > jet2)
 
double calcMT_1l (HEPUtils::P4 jetMom, HEPUtils::P4 metMom)
 
bool sortByPT_j (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByPT_l (HEPUtils::Particle *lep1, HEPUtils::Particle *lep2)
 
bool sortByPT_jet (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByPT_lep (HEPUtils::Particle *lep1, HEPUtils::Particle *lep2)
 
bool sortByPT_RJ3L (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortLepByPT_RJ3L (HEPUtils::Particle *lep1, HEPUtils::Particle *lep2)
 
bool SortLeptons (const pair< TLorentzVector, int > lv1, const pair< TLorentzVector, int > lv2)
 
bool SortJets (const TLorentzVector jv1, const TLorentzVector jv2)
 
 DEFINE_ANALYSIS_FACTORY (ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb)
 
bool sortByPT (HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
 
bool sortByPT_2lep (HEPUtils::Particle *lep1, HEPUtils::Particle *lep2)
 
double Phi_mpi_pi (double x)
 
double _Phi_mpi_pi (double x)
 
 MAP_ANALYSES_WITH_ROOT_RESTFRAMES (DECLARE_ANALYSIS_FACTORY) MAP_ANALYSES_WITH_ROOT(DECLARE_ANALYSIS_FACTORY) MAP_ANALYSES(DECLARE_ANALYSIS_FACTORY) Analysis *mkAnalysis(const str &name)
 Forward declarations using DECLARE_ANALYSIS_FACTORY(ANAME) More...
 
void getAnalysisContainer (AnalysisContainer &result, const str &detname, const MCLoopInfo &RunMC, const xsec &CrossSection, int iteration)
 Retrieve an analysis container for a specific detector. More...
 
 GET_ANALYSIS_CONTAINER (getATLASAnalysisContainer, ATLAS) GET_ANALYSIS_CONTAINER(getCMSAnalysisContainer
 
void runAnalyses (AnalysisDataPointers &result, const str &, const MCLoopInfo &RunMC, const AnalysisContainer &Container, const HEPUtils::Event &SmearedEvent, int iteration, void(*wrapup)())
 Run all the analyses in a given container. More...
 
 RUN_ANALYSES (runATLASAnalyses, ATLAS, ATLASSmearedEvent) RUN_ANALYSES(runCMSAnalyses
 
void getDummyColliderObservable (double &result)
 Dummy observable that creates a dependency on TestModel1D. More...
 
void operateLHCLoop (MCLoopInfo &result)
 LHC Loop Manager. More...
 
void getLHCEventLoopInfo (map_str_dbl &result)
 Store some information about the event generation. More...
 
void CollectAnalyses (AnalysisDataPointers &result)
 Loop over all analyses and collect them in one place. More...
 
void set_CS (hb_ModelParameters &result, const HiggsCouplingsTable &couplings, int n_neutral_higgses)
 Helper function to set HiggsBounds/Signals parameters cross-section ratios from a GAMBIT HiggsCouplingsTable. More...
 
void set_SMLikeHiggs_ModelParameters (const SubSpectrum &spec, const HiggsCouplingsTable &couplings, hb_ModelParameters &result)
 Helper function for populating a HiggsBounds/Signals ModelParameters object for SM-like Higgs. More...
 
void SMHiggs_ModelParameters (hb_ModelParameters &result)
 SM Higgs model parameters for HiggsBounds/Signals. More...
 
void SMLikeHiggs_ModelParameters (hb_ModelParameters &result)
 SM-like (SM + possible invisibles) Higgs model parameters for HiggsBounds/Signals. More...
 
void MSSMHiggs_ModelParameters (hb_ModelParameters &result)
 MSSM Higgs model parameters. More...
 
void calc_HB_LEP_LogLike (double &result)
 Get a LEP chisq from HiggsBounds. More...
 
void calc_HS_LHC_LogLike (double &result)
 Get an LHC chisq from HiggsSignals. More...
 
void FH_HiggsProd (fh_HiggsProd &result)
 Higgs production cross-sections from FeynHiggs. More...
 
double limit_LLike (double x, double x95, double sigma)
 LEP limit likelihood function. More...
 
bool is_xsec_sane (const triplet< double > &xsecWithError)
 LEP limit debugging function. More...
 
void LEP207_SLHA1_convention_xsec_chi00_11 (triplet< double > &result)
 
void L3_Gravitino_LLike (double &result)
 
void getMCxsec (xsec &result)
 Compute a cross-section from Monte Carlo. More...
 
void getNLLFastxsec (xsec &result)
 Get a cross-section from NLL-FAST. More...
 
void calc_LHC_signals (map_str_dbl &result)
 Loop over all analyses and fill a map of predicted counts. More...
 
void calc_LHC_LogLikes (map_str_AnalysisLogLikes &result)
 Loop over all analyses and fill a map of AnalysisLogLikes objects. More...
 
void get_LHC_LogLike_per_analysis (map_str_dbl &result)
 Extract the combined log likelihood for each analysis. More...
 
void get_LHC_LogLike_per_SR (map_str_dbl &result)
 Extract the log likelihood for each SR. More...
 
void get_LHC_LogLike_SR_labels (map_str_str &result)
 Extract the labels for the SRs used in the analysis loglikes. More...
 
void get_LHC_LogLike_SR_indices (map_str_dbl &result)
 Extract the indices for the SRs used in the analysis loglikes. More...
 
void calc_combined_LHC_LogLike (double &result)
 Compute the total likelihood combining all analyses. More...
 
 GET_SPECIFIC_PYTHIA (getPythia_EM, Pythia_EM_default, EM_spectrum, _EM, NOT_SUSY) GET_PYTHIA_AS_BASE_COLLIDER(getPythia_EMAsBase) GET_PYTHIA_EVENT(generateEventPythia_EM
 
Pythia_EM_default::Pythia8::Event GET_BUCKFAST_AS_BASE_DETECTOR (getBuckFastATLASPythia_EM, Pythia_EM_default::Pythia8::Event, ATLAS) GET_BUCKFAST_AS_BASE_DETECTOR(getBuckFastCMSPythia_EM
 
Pythia_EM_default::Pythia8::Event CMS GET_BUCKFAST_AS_BASE_DETECTOR (getBuckFastIdentityPythia_EM, Pythia_EM_default::Pythia8::Event, Identity) SMEAR_EVENT(smearEventATLAS_EM
 
Pythia_EM_default::Pythia8::Event CMS ATLAS SMEAR_EVENT (smearEventCMS_EM, CMS) SMEAR_EVENT(copyEvent_EM
 
 GET_SPECIFIC_PYTHIA (getPythia, Pythia_default, MSSM_spectrum,, IS_SUSY) GET_PYTHIA_AS_BASE_COLLIDER(getPythiaAsBase) GET_PYTHIA_EVENT(generateEventPythia
 
Pythia_default::Pythia8::Event GET_BUCKFAST_AS_BASE_DETECTOR (getBuckFastATLASPythia, Pythia_default::Pythia8::Event, ATLAS) GET_BUCKFAST_AS_BASE_DETECTOR(getBuckFastCMSPythia
 
Pythia_default::Pythia8::Event CMS GET_BUCKFAST_AS_BASE_DETECTOR (getBuckFastIdentityPythia, Pythia_default::Pythia8::Event, Identity) SMEAR_EVENT(smearEventATLAS
 
Pythia_default::Pythia8::Event CMS ATLAS SMEAR_EVENT (smearEventCMS, CMS) SMEAR_EVENT(copyEvent
 
Converters to/from Pythia8's native 4-vector
template<typename Vec4T >
FJNS::PseudoJet mk_pseudojet (const Vec4T &p)
 
template<typename Vec4T >
HEPUtils::P4 mk_p4 (const Vec4T &p)
 
Detailed Pythia8 event record walking/mangling functions
template<typename EventT >
bool fromBottom (int n, const EventT &evt)
 
template<typename EventT >
bool fromTau (int n, const EventT &evt)
 
template<typename EventT >
bool fromHadron (int n, const EventT &evt)
 
template<typename EventT >
bool isFinalB (int n, const EventT &evt)
 
template<typename EventT >
bool isFinalTau (int n, const EventT &evt)
 
template<typename EventT >
bool isParton (int n, const EventT &evt)
 
template<typename EventT >
bool isFinalParton (int n, const EventT &evt)
 
template<typename EventT >
bool isFinalPhoton (int n, const EventT &evt)
 
template<typename EventT >
bool isFinalLepton (int n, const EventT &evt)
 
void get_sigma_ee_ll (triplet< double > &result, const double sqrts, const int generation, const int l_chirality, const int lbar_chirality, const double gtol, const double ftol, const bool gpt_error, const bool fpt_error, const Spectrum &spec, const double gammaZ, const bool l_are_gauge_es)
 High-level cross section routines. More...
 
void get_sigma_ee_chi00 (triplet< double > &result, const double sqrts, const int chi_first, const int chi_second, const double tol, const bool pt_error, const Spectrum &spec, const double gammaZ)
 Retrieve the production cross-section at an e+e- collider for neutralino pairs. More...
 
void get_sigma_ee_chipm (triplet< double > &result, const double sqrts, const int chi_plus, const int chi_minus, const double tol, const bool pt_error, const Spectrum &spec, const double gammaZ)
 Retrieve the production cross-section at an e+e- collider for chargino pairs. More...
 
double xsec_sleislej (int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix F, MixMatrix N, const double mN[4], double alpha, double mZ, double gZ, double sin2thetaW, bool warn_on_CP_violating_masses=true)
 Low-level cross section routines. More...
 
double xsec_neuineuj (int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix N, const double mS[2], double tanb, double alpha, double mZ, double gZ, double sin2thetaW)
 Cross section [pb] for $ e^+e^- -> \tilde\chi^0_i \tilde\chi^0_j $ Masses mi and mj for the neutralinos are signed. More...
 
double xsec_chaichaj (int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix V, MixMatrix U, double msn, double alpha, double mZ, double gZ, double sin2thetaW)
 Cross section [pb] for $ e^+e^- -> \tilde\chi^+_i \tilde\chi^-_j $ Masses mi and mj for the charginos are signed. More...
 
void SLHA2BFM_NN (MixMatrix &NN, double tanb, double sin2thetaW)
 Conversion between SLHA and BFM conventions. More...
 
void SLHA2BFM_VV (MixMatrix &VV)
 Converts the chargino mixing matrix V in SLHA conventions to BFM conventions. More...
 
void BFM2SLHA_NN (MixMatrix &NN, double tanb, double sin2thetaW)
 Converts a neutralino mixing matrix in BFM conventions to SLHA conventions, $\tan\beta$ is as defined in SLHA. More...
 
void BFM2SLHA_VV (MixMatrix &VV)
 Converts the chargino mixing matrix V in BFM conventions to SLHA conventions. More...
 
MixMatrix multiply (MixMatrix A, MixMatrix B)
 Helper function to multiply matrices. More...
 
MixMatrix transpose (MixMatrix A)
 Helper function to find matrix transpose. More...
 
void print (MixMatrix A)
 Helper function to print a matrix. More...
 
External operators and string representation for P2
P2 operator+ (const P2 &a, const P2 &b)
 
P2 operator- (const P2 &a, const P2 &b)
 
P2 operator* (const P2 &a, double f)
 
P2 operator* (double f, const P2 &a)
 
P2 operator/ (const P2 &a, double f)
 
std::string to_str (const P2 &p2)
 Make a string representation of the vector. More...
 
std::ostream & operator<< (std::ostream &ostr, const P2 &p2)
 Write a string representation of the vector to the provided stream. More...
 
String representation for LineSegment
std::string to_str (const LineSegment &lineseg)
 Make a string representation of the LineSegment. More...
 
std::ostream & operator<< (std::ostream &ostr, const LineSegment &lineseg)
 Write a string representation of the vector to the provided stream. More...
 
Particle-list filtering by cuts
template<typename CONTAINER , typename RMFN >
void iremoveerase (CONTAINER &c, const RMFN &fn)
 Convenience combination of remove_if and erase. More...
 
void ifilter_reject (ParticlePtrs &particles, std::function< bool(Particle *)> rejfn, bool do_delete=true)
 In-place filter a supplied particle vector by rejecting those which fail a supplied cut. More...
 
void ifilter_select (ParticlePtrs &particles, std::function< bool(Particle *)> selfn, bool do_delete=true)
 In-place filter a supplied particle vector by keeping those which pass a supplied cut. More...
 
ParticlePtrs filter_reject (const ParticlePtrs &particles, std::function< bool(Particle *)> rejfn, bool do_delete=true)
 Filter a supplied particle vector by rejecting those which fail a supplied cut. More...
 
ParticlePtrs filter_select (const ParticlePtrs &particles, std::function< bool(Particle *)> selfn, bool do_delete=true)
 Filter a supplied particle vector by keeping those which pass a supplied cut. More...
 
Jet-list filtering by cuts
void ifilter_reject (JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
 In-place filter a supplied jet vector by rejecting those which fail a supplied cut. More...
 
void ifilter_select (JetPtrs &jets, std::function< bool(Jet *)> selfn, bool do_delete=true)
 In-place filter a supplied jet vector by keeping those which pass a supplied cut. More...
 
JetPtrs filter_reject (const JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
 Filter a supplied particle vector by rejecting those which fail a supplied cut. More...
 
JetPtrs filter_select (const JetPtrs &jets, std::function< bool(Jet *)> selfn, bool do_delete=true)
 Filter a supplied particle vector by keeping those which pass a supplied cut. More...
 
Random booleans sampled from efficiency maps
bool random_bool (double eff)
 Return a random true/false at a success rate given by a number. More...
 
bool random_bool (const HEPUtils::BinnedFn1D< double > &effmap, double x)
 Return a random true/false at a success rate given by a 1D efficiency map. More...
 
bool random_bool (const HEPUtils::BinnedFn2D< double > &effmap, double x, double y)
 Return a random true/false at a success rate given by a 2D efficiency map. More...
 
Random filtering by efficiency
void filtereff (std::vector< HEPUtils::Particle *> &particles, double eff, bool do_delete=false)
 Utility function for filtering a supplied particle vector by sampling wrt an efficiency scalar. More...
 
void filtereff (std::vector< HEPUtils::Particle *> &particles, std::function< double(HEPUtils::Particle *)> eff_fn, bool do_delete=false)
 Utility function for filtering a supplied particle vector by sampling an efficiency returned by a provided function object. More...
 
void filtereff_pt (std::vector< HEPUtils::Particle *> &particles, const HEPUtils::BinnedFn1D< double > &eff_pt, bool do_delete=false)
 Utility function for filtering a supplied particle vector by sampling wrt a binned 1D efficiency map in pT. More...
 
void filtereff_etapt (std::vector< HEPUtils::Particle *> &particles, const HEPUtils::BinnedFn2D< double > &eff_etapt, bool do_delete=false)
 Utility function for filtering a supplied particle vector by sampling wrt a binned 2D efficiency map in |eta| and pT. More...
 
Tagging
bool has_tag (const HEPUtils::BinnedFn2D< double > &effmap, double eta, double pt)
 Randomly get a tag result (can be anything) from a 2D |eta|-pT efficiency map. More...
 
template<typename NUM1 , typename NUM2 >
size_t binIndex (NUM1 val, const std::vector< NUM2 > &binedges, bool allow_overflow=false)
 
std::vector< doublemk_bin_values (const std::vector< double > &binEdgeValues)
 Make a vector of central bin values from a vector of bin edge values using linear interpolation. More...
 
std::vector< doublemakeBinValues (const std::vector< double > &binEdgeValues)
 Alias. More...
 
template<typename MOM >
std::vector< std::shared_ptr< HEPUtils::Jet > > get_jets (const std::vector< MOM *> &moms, double R, double ptmin=0 *GeV, FJNS::JetAlgorithm alg=FJNS::antikt_algorithm)
 Run jet clustering from any P4-compatible momentum type. More...
 
bool object_in_cone (const HEPUtils::Event &e, const HEPUtils::P4 &p4, double ptmin, double rmax, double rmin=0.05)
 Check if there's a physics object above ptmin in an annulus rmin..rmax around the given four-momentum p4. More...
 
template<typename MOMPTRS1 , typename MOMPTRS2 >
void removeOverlap (MOMPTRS1 &momstofilter, const MOMPTRS2 &momstocmp, double deltaRMax, bool use_rapidity=false)
 Overlap removal – discard from first list if within deltaRMax of any from the second list. More...
 
template<typename CONTAINER , typename FN >
bool all_of (const CONTAINER &c, const FN &f)
 Non-iterator version of std::all_of. More...
 
template<typename CONTAINER , typename FN >
bool any_of (const CONTAINER &c, const FN &f)
 Non-iterator version of std::any_of. More...
 
template<typename CONTAINER , typename FN >
bool none_of (const CONTAINER &c, const FN &f)
 Non-iterator version of std::none_of. More...
 
std::vector< std::vector< HEPUtils::Particle * > > getSFOSpairs (std::vector< HEPUtils::Particle *> particles)
 Utility function for returning a collection of same-flavour, oppsosite-sign particle pairs. More...
 
std::vector< std::vector< HEPUtils::Particle * > > getOSpairs (std::vector< HEPUtils::Particle *> particles)
 Utility function for returning a collection of oppsosite-sign particle pairs. More...
 
std::vector< std::vector< HEPUtils::Particle * > > getSSpairs (std::vector< HEPUtils::Particle *> particles)
 Utility function for returning a collection of same-sign particle pairs. More...
 
Sorting
void sortBy (ParticlePtrs &particles, std::function< bool(const Particle *, const Particle *)> cmpfn)
 Particle-sorting function. More...
 
bool cmpParticlesByPt (const HEPUtils::Particle *lep1, const HEPUtils::Particle *lep2)
 Comparison function to give a particles sorting order decreasing by pT. More...
 
void sortByPt (ParticlePtrs &particles)
 
void sortBy (JetPtrs &jets, std::function< bool(const Jet *, const Jet *)> cmpfn)
 Jet-sorting function. More...
 
bool cmpJetsByPt (const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
 Comparison function to give a jets sorting order decreasing by pT. More...
 
void LEP208_SLHA1_convention_xsec_selselbar (triplet< double > &result)
 ee –> selectron pair production cross-sections at 208 GeV More...
 
void LEP208_SLHA1_convention_xsec_selserbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_serserbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_serselbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_se1se1bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_se1se2bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_se2se2bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_se2se1bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smulsmulbar (triplet< double > &result)
 ee –> smuon pair production cross-sections at 208 GeV More...
 
void LEP208_SLHA1_convention_xsec_smulsmurbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smursmurbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smursmulbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smu1smu1bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smu1smu2bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smu2smu2bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_smu2smu1bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_staulstaulbar (triplet< double > &result)
 ee –> stau pair production cross-sections at 208 GeV More...
 
void LEP208_SLHA1_convention_xsec_staulstaurbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_staurstaurbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_staurstaulbar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_stau1stau1bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_stau1stau2bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_stau2stau2bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_stau2stau1bar (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_11 (triplet< double > &result)
 ee –> neutralino pair production cross-sections at 208 GeV More...
 
void LEP208_SLHA1_convention_xsec_chi00_12 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_13 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_14 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_22 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_23 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_24 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_33 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_34 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chi00_44 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chipm_11 (triplet< double > &result)
 ee –> chargino pair production cross-sections at 208 GeV More...
 
void LEP208_SLHA1_convention_xsec_chipm_12 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chipm_22 (triplet< double > &result)
 
void LEP208_SLHA1_convention_xsec_chipm_21 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_selselbar (triplet< double > &result)
 ee –> selectron pair production cross-sections at 205 GeV More...
 
void LEP205_SLHA1_convention_xsec_selserbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_serserbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_serselbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_se1se1bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_se1se2bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_se2se2bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_se2se1bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smulsmulbar (triplet< double > &result)
 ee –> smuon pair production cross-sections at 205 GeV More...
 
void LEP205_SLHA1_convention_xsec_smulsmurbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smursmurbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smursmulbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smu1smu1bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smu1smu2bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smu2smu2bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_smu2smu1bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_staulstaulbar (triplet< double > &result)
 ee –> stau pair production cross-sections at 205 GeV More...
 
void LEP205_SLHA1_convention_xsec_staulstaurbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_staurstaurbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_staurstaulbar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_stau1stau1bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_stau1stau2bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_stau2stau2bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_stau2stau1bar (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_11 (triplet< double > &result)
 ee –> neutralino pair production cross-sections at 205 GeV More...
 
void LEP205_SLHA1_convention_xsec_chi00_12 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_13 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_14 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_22 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_23 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_24 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_33 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_34 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chi00_44 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chipm_11 (triplet< double > &result)
 ee –> chargino pair production cross-sections at 205 GeV More...
 
void LEP205_SLHA1_convention_xsec_chipm_12 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chipm_22 (triplet< double > &result)
 
void LEP205_SLHA1_convention_xsec_chipm_21 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_selselbar (triplet< double > &result)
 ee –> selectron pair production cross-sections at 188.6 GeV More...
 
void LEP188_SLHA1_convention_xsec_selserbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_serserbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_serselbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_se1se1bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_se1se2bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_se2se2bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_se2se1bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smulsmulbar (triplet< double > &result)
 ee –> smuon pair production cross-sections at 188.6 GeV More...
 
void LEP188_SLHA1_convention_xsec_smulsmurbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smursmurbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smursmulbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smu1smu1bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smu1smu2bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smu2smu2bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_smu2smu1bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_staulstaulbar (triplet< double > &result)
 ee –> stau pair production cross-sections at 188.6 GeV More...
 
void LEP188_SLHA1_convention_xsec_staulstaurbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_staurstaurbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_staurstaulbar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_stau1stau1bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_stau1stau2bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_stau2stau2bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_stau2stau1bar (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_11 (triplet< double > &result)
 ee –> neutralino pair production cross-sections at 188.6 GeV More...
 
void LEP188_SLHA1_convention_xsec_chi00_12 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_13 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_14 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_22 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_23 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_24 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_33 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_34 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chi00_44 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chipm_11 (triplet< double > &result)
 ee –> chargino pair production cross-sections at 188.6 GeV More...
 
void LEP188_SLHA1_convention_xsec_chipm_12 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chipm_22 (triplet< double > &result)
 
void LEP188_SLHA1_convention_xsec_chipm_21 (triplet< double > &result)
 
void ALEPH_Selectron_Conservative_LLike (double &result)
 LEP Slepton Log-Likelihoods. More...
 
void ALEPH_Smuon_Conservative_LLike (double &result)
 
void ALEPH_Stau_Conservative_LLike (double &result)
 
void L3_Selectron_Conservative_LLike (double &result)
 
void L3_Smuon_Conservative_LLike (double &result)
 
void L3_Stau_Conservative_LLike (double &result)
 
void L3_Neutralino_All_Channels_Conservative_LLike (double &result)
 LEP Gaugino Log-Likelihoods. More...
 
void L3_Neutralino_Leptonic_Conservative_LLike (double &result)
 
void L3_Chargino_All_Channels_Conservative_LLike (double &result)
 
void L3_Chargino_Leptonic_Conservative_LLike (double &result)
 
void OPAL_Chargino_Hadronic_Conservative_LLike (double &result)
 
void OPAL_Chargino_SemiLeptonic_Conservative_LLike (double &result)
 
void OPAL_Chargino_Leptonic_Conservative_LLike (double &result)
 
void OPAL_Degenerate_Chargino_Conservative_LLike (double &result)
 
void OPAL_Chargino_All_Channels_Conservative_LLike (double &result)
 
void OPAL_Neutralino_Hadronic_Conservative_LLike (double &result)
 
double I1 (double s, double m1, double m2, double mk, double ml)
 Integrals for t-channel neutralino diagrams m1 and m2 are masses of final state sleptons mk and ml are neutralino masses. More...
 
double I2 (double s, double m1, double m2, double mk, double ml)
 
double I3 (double s, double m1, double m2, double mk)
 

Variables

 CMS
 

Typedef Documentation

◆ AnalysisContainers

Container for multiple analysis containers.

Definition at line 92 of file ColliderBit_types.hpp.

◆ AnalysisDataPointers

Definition at line 69 of file ColliderBit_types.hpp.

◆ AnalysisNumbers

Container for data from multiple analyses and SRs.

Definition at line 68 of file ColliderBit_types.hpp.

◆ ConstJetPtrs

typedef std::vector<const HEPUtils::Jet*> Gambit::ColliderBit::ConstJetPtrs

Typedef for a vector of const Jet pointers.

Definition at line 39 of file Utils.hpp.

◆ ConstParticlePtrs

typedef std::vector<const HEPUtils::Particle*> Gambit::ColliderBit::ConstParticlePtrs

Typedef for a vector of const Particle pointers.

Definition at line 33 of file Utils.hpp.

◆ data_type

typedef std::vector<std::vector<double> > Gambit::ColliderBit::data_type

Definition at line 33 of file ImageLimit.hpp.

◆ JetPtrs

typedef std::vector<HEPUtils::Jet*> Gambit::ColliderBit::JetPtrs

Typedef for a vector of Jet pointers.

Definition at line 37 of file Utils.hpp.

◆ map_str_AnalysisLogLikes

Typedef for a string-to-AnalysisLogLikes map.

Definition at line 89 of file ColliderBit_types.hpp.

◆ MixMatrix

typedef std::vector< std::vector<double> > Gambit::ColliderBit::MixMatrix

Definition at line 56 of file lep_mssm_xsecs.hpp.

◆ ms

typedef std::chrono::milliseconds Gambit::ColliderBit::ms

Definition at line 94 of file ColliderBit_types.hpp.

◆ ParticlePtrs

typedef std::vector<HEPUtils::Particle*> Gambit::ColliderBit::ParticlePtrs

Typedef for a vector of Particle pointers.

Definition at line 31 of file Utils.hpp.

◆ steady_clock

typedef std::chrono::steady_clock Gambit::ColliderBit::steady_clock

Definition at line 95 of file ColliderBit_types.hpp.

◆ timer_map_type

typedef std::map<std::string,double> Gambit::ColliderBit::timer_map_type

Definition at line 97 of file ColliderBit_types.hpp.

◆ tp

typedef std::chrono::steady_clock::time_point Gambit::ColliderBit::tp

Definition at line 96 of file ColliderBit_types.hpp.

Enumeration Type Documentation

◆ specialIterations

Function Documentation

◆ _Phi_mpi_pi()

◆ addInQuad()

double Gambit::ColliderBit::addInQuad ( const double a,
const double b 
)
inline

Add two numbers in quadrature.

Definition at line 13 of file PointsAndLines.hpp.

Referenced by Gambit::ColliderBit::P2::abs(), and Gambit::ColliderBit::P2::r().

13 { return sqrt(a*a + b*b); }
START_MODEL b
Definition: demo.hpp:235
Here is the caller graph for this function:

◆ ALEPH_Selectron_Conservative_LLike()

void Gambit::ColliderBit::ALEPH_Selectron_Conservative_LLike ( double result)

LEP Slepton Log-Likelihoods.

Definition at line 1210 of file ColliderBit_LEP.cpp.

References ALEPH_Selectron_Conservative_LLike, Gambit::triplet< TYPE >::central, Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), LEP208_xsec_selselbar, LEP208_xsec_serserbar, limit_LLike(), Gambit::ColliderBit::BaseLimitContainer::limitAverage(), Gambit::triplet< TYPE >::lower, Gambit::slhahelp::mass_es_from_gauge_es(), MSSM_spectrum, Gambit::Par::Pole_Mass, Gambit::Scanner::pow(), selectron_l_decay_rates, selectron_r_decay_rates, and Gambit::triplet< TYPE >::upper.

1211  {
1213  using std::pow;
1214  using std::log;
1215 
1216  const Spectrum& spec = *Dep::MSSM_spectrum;
1217 
1218  double max_mixing;
1219  const SubSpectrum& mssm = spec.get_HE();
1220  str sel_string = slhahelp::mass_es_from_gauge_es("~e_L", max_mixing, mssm);
1221  str ser_string = slhahelp::mass_es_from_gauge_es("~e_R", max_mixing, mssm);
1222  const double mass_seL=spec.get(Par::Pole_Mass,sel_string);
1223  const double mass_neut1 = spec.get(Par::Pole_Mass,1000022, 0);
1224  const double mass_seR = spec.get(Par::Pole_Mass,ser_string);
1225  const double mZ = spec.get(Par::Pole_Mass,23, 0);
1226  triplet<double> xsecWithError;
1227  double xsecLimit;
1228 
1229  static const ALEPHSelectronLimitAt208GeV limitContainer;
1230  // #ifdef COLLIDERBIT_DEBUG
1231  //static bool dumped=false;
1232  //if(!dumped)
1233  //{
1234  // limitContainer.dumpPlotData(45., 115., 0., 100., mZ, "lepLimitPlanev2/ALEPHSelectronLimitAt208GeV.dump");
1235  // dumped=true;
1236  //}
1237  // #endif
1238 
1239  result = 0;
1240  // Due to the nature of the analysis details of the model independent limit in
1241  // the paper, the best we can do is to try these two processes individually:
1242 
1243  // se_L, se_L
1244  xsecLimit = limitContainer.limitAverage(mass_seL, mass_neut1, mZ);
1245  xsecWithError = *Dep::LEP208_xsec_selselbar;
1246  xsecWithError.upper *= pow(Dep::selectron_l_decay_rates->BF("~chi0_1", "e-"), 2);
1247  xsecWithError.central *= pow(Dep::selectron_l_decay_rates->BF("~chi0_1", "e-"), 2);
1248  xsecWithError.lower *= pow(Dep::selectron_l_decay_rates->BF("~chi0_1", "e-"), 2);
1249 
1250  if (xsecWithError.central < xsecLimit)
1251  {
1252  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1253  }
1254  else
1255  {
1256  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1257  }
1258 
1259  // se_R, se_R
1260  xsecLimit = limitContainer.limitAverage(mass_seR, mass_neut1, mZ);
1261 
1262  xsecWithError = *Dep::LEP208_xsec_serserbar;
1263  xsecWithError.upper *= pow(Dep::selectron_r_decay_rates->BF("~chi0_1", "e-"), 2);
1264  xsecWithError.central *= pow(Dep::selectron_r_decay_rates->BF("~chi0_1", "e-"), 2);
1265  xsecWithError.lower *= pow(Dep::selectron_r_decay_rates->BF("~chi0_1", "e-"), 2);
1266 
1267  if (xsecWithError.central < xsecLimit)
1268  {
1269  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1270  }
1271  else
1272  {
1273  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1274  }
1275 
1276  }
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 ALEPH_Selectron_Conservative_LLike
double limit_LLike(double x, double x95, double sigma)
LEP limit likelihood function.
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 selectron_r_decay_rates
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 selectron_l_decay_rates
Spectrum Spectrum
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
double pow(const double &a)
Outputs a^i.
DecayTable::Entry LEP208_xsec_serserbar
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
Here is the call graph for this function:

◆ ALEPH_Smuon_Conservative_LLike()

void Gambit::ColliderBit::ALEPH_Smuon_Conservative_LLike ( double result)

Definition at line 1278 of file ColliderBit_LEP.cpp.

References Gambit::triplet< TYPE >::central, Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), LEP208_xsec_smulsmulbar, LEP208_xsec_smursmurbar, limit_LLike(), Gambit::ColliderBit::BaseLimitContainer::limitAverage(), Gambit::triplet< TYPE >::lower, Gambit::slhahelp::mass_es_from_gauge_es(), MSSM_spectrum, Gambit::Par::Pole_Mass, Gambit::Scanner::pow(), smuon_l_decay_rates, smuon_r_decay_rates, and Gambit::triplet< TYPE >::upper.

1279  {
1280  using namespace Pipes::ALEPH_Smuon_Conservative_LLike;
1281  using std::pow;
1282  using std::log;
1283 
1284  const Spectrum& spec = *Dep::MSSM_spectrum;
1285 
1286  double max_mixing;
1287  const SubSpectrum& mssm = spec.get_HE();
1288  str smul_string = slhahelp::mass_es_from_gauge_es("~mu_L", max_mixing, mssm);
1289  str smur_string = slhahelp::mass_es_from_gauge_es("~mu_R", max_mixing, mssm);
1290  const double mass_smuL=spec.get(Par::Pole_Mass,smul_string);
1291  const double mass_neut1 = spec.get(Par::Pole_Mass,1000022, 0);
1292  const double mass_smuR = spec.get(Par::Pole_Mass,smur_string);
1293  const double mZ = spec.get(Par::Pole_Mass,23, 0);
1294  triplet<double> xsecWithError;
1295  double xsecLimit;
1296 
1297  static const ALEPHSmuonLimitAt208GeV limitContainer;
1298  // #ifdef COLLIDERBIT_DEBUG
1299  // static bool dumped=false;
1300  // if(!dumped)
1301  // {
1302  // limitContainer.dumpPlotData(45., 115., 0., 100., mZ, "lepLimitPlanev2/ALEPHSmuonLimitAt208GeV.dump");
1303  // dumped=true;
1304  // }
1305  // #endif
1306 
1307  result = 0;
1308  // Due to the nature of the analysis details of the model independent limit in
1309  // the paper, the best we can do is to try these two processes individually:
1310 
1311  // smu_L, smu_L
1312  xsecLimit = limitContainer.limitAverage(mass_smuL, mass_neut1, mZ);
1313  xsecWithError = *Dep::LEP208_xsec_smulsmulbar;
1314  xsecWithError.upper *= pow(Dep::smuon_l_decay_rates->BF("~chi0_1", "mu-"), 2);
1315  xsecWithError.central *= pow(Dep::smuon_l_decay_rates->BF("~chi0_1", "mu-"), 2);
1316  xsecWithError.lower *= pow(Dep::smuon_l_decay_rates->BF("~chi0_1", "mu-"), 2);
1317 
1318  if (xsecWithError.central < xsecLimit)
1319  {
1320  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1321  }
1322  else
1323  {
1324  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1325  }
1326 
1327  // smu_R, smu_R
1328  xsecLimit = limitContainer.limitAverage(mass_smuR, mass_neut1, mZ);
1329 
1330  xsecWithError = *Dep::LEP208_xsec_smursmurbar;
1331  xsecWithError.upper *= pow(Dep::smuon_r_decay_rates->BF("~chi0_1", "mu-"), 2);
1332  xsecWithError.central *= pow(Dep::smuon_r_decay_rates->BF("~chi0_1", "mu-"), 2);
1333  xsecWithError.lower *= pow(Dep::smuon_r_decay_rates->BF("~chi0_1", "mu-"), 2);
1334 
1335  if (xsecWithError.central < xsecLimit)
1336  {
1337  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1338  }
1339  else
1340  {
1341  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1342  }
1343 
1344  }
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry LEP208_xsec_smulsmulbar
double limit_LLike(double x, double x95, double sigma)
LEP limit likelihood function.
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry LEP208_xsec_smursmurbar
Spectrum Spectrum
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
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 DecayTable::Entry DecayTable::Entry smuon_l_decay_rates
double pow(const double &a)
Outputs a^i.
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
void ALEPH_Smuon_Conservative_LLike(double &result)
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 DecayTable::Entry DecayTable::Entry smuon_r_decay_rates
Here is the call graph for this function:

◆ ALEPH_Stau_Conservative_LLike()

void Gambit::ColliderBit::ALEPH_Stau_Conservative_LLike ( double result)

Definition at line 1346 of file ColliderBit_LEP.cpp.

References ALEPH_Stau_Conservative_LLike, Gambit::triplet< TYPE >::central, Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), LEP208_xsec_stau1stau1bar, LEP208_xsec_stau2stau2bar, limit_LLike(), Gambit::ColliderBit::BaseLimitContainer::limitAverage(), LOCAL_INFO, Gambit::triplet< TYPE >::lower, Gambit::slhahelp::mass_es_closest_to_family(), MSSM_spectrum, Gambit::Par::Pole_Mass, Gambit::Scanner::pow(), stau_1_decay_rates, stau_2_decay_rates, and Gambit::triplet< TYPE >::upper.

1347  {
1348  using namespace Pipes::ALEPH_Stau_Conservative_LLike;
1349  using std::pow;
1350  using std::log;
1351 
1352  const Spectrum& spec = *Dep::MSSM_spectrum;
1353 
1354  const SubSpectrum& mssm = spec.get_HE();
1355  static const double tol = runOptions->getValueOrDef<double>(1e-5, "family_mixing_tolerance");
1356  static const bool pterror = runOptions->getValueOrDef<bool>(false, "family_mixing_tolerance_invalidates_point_only");
1357  str stau1_string = slhahelp::mass_es_closest_to_family("~tau_1", mssm,tol,LOCAL_INFO,pterror);
1358  str stau2_string = slhahelp::mass_es_closest_to_family("~tau_2", mssm,tol,LOCAL_INFO,pterror);
1359  const double mass_stau1=spec.get(Par::Pole_Mass,stau1_string);
1360  const double mass_neut1 = spec.get(Par::Pole_Mass,1000022, 0);
1361  const double mass_stau2 = spec.get(Par::Pole_Mass,stau2_string);
1362  const double mZ = spec.get(Par::Pole_Mass,23, 0);
1363  triplet<double> xsecWithError;
1364  double xsecLimit;
1365 
1366  static const ALEPHStauLimitAt208GeV limitContainer;
1367  // #ifdef COLLIDERBIT_DEBUG
1368  // static bool dumped=false;
1369  // if(!dumped)
1370  // {
1371  // limitContainer.dumpPlotData(45., 115., 0., 100., mZ, "lepLimitPlanev2/ALEPHStauLimitAt208GeV.dump");
1372  // dumped=true;
1373  // }
1374  // #endif
1375 
1376  result = 0;
1377  // Due to the nature of the analysis details of the model independent limit in
1378  // the paper, the best we can do is to try these two processes individually:
1379 
1380  // stau_1, stau_1
1381  xsecLimit = limitContainer.limitAverage(mass_stau1, mass_neut1, mZ);
1382 
1383  xsecWithError = *Dep::LEP208_xsec_stau1stau1bar;
1384  xsecWithError.upper *= pow(Dep::stau_1_decay_rates->BF("~chi0_1", "tau-"), 2);
1385  xsecWithError.central *= pow(Dep::stau_1_decay_rates->BF("~chi0_1", "tau-"), 2);
1386  xsecWithError.lower *= pow(Dep::stau_1_decay_rates->BF("~chi0_1", "tau-"), 2);
1387 
1388  if (xsecWithError.central < xsecLimit)
1389  {
1390  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1391  }
1392  else
1393  {
1394  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1395  }
1396 
1397  // stau_2, stau_2
1398  xsecLimit = limitContainer.limitAverage(mass_stau2, mass_neut1, mZ);
1399 
1400  xsecWithError = *Dep::LEP208_xsec_stau2stau2bar;
1401  xsecWithError.upper *= pow(Dep::stau_2_decay_rates->BF("~chi0_1", "tau-"), 2);
1402  xsecWithError.central *= pow(Dep::stau_2_decay_rates->BF("~chi0_1", "tau-"), 2);
1403  xsecWithError.lower *= pow(Dep::stau_2_decay_rates->BF("~chi0_1", "tau-"), 2);
1404 
1405  if (xsecWithError.central < xsecLimit)
1406  {
1407  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1408  }
1409  else
1410  {
1411  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1412  }
1413 
1414  }
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 DecayTable::Entry stau_2_decay_rates
#define LOCAL_INFO
Definition: local_info.hpp:34
double limit_LLike(double x, double x95, double sigma)
LEP limit likelihood function.
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry LEP208_xsec_stau1stau1bar
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 DecayTable::Entry stau_1_decay_rates
Spectrum Spectrum
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
str mass_es_closest_to_family(str familystate, const SubSpectrum &mssm)
identify the mass eigenstate corresponding to family state takes string and returns only requested st...
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry LEP208_xsec_stau2stau2bar
double pow(const double &a)
Outputs a^i.
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 DecayTable::Entry ALEPH_Stau_Conservative_LLike
Here is the call graph for this function:

◆ all_of()

template<typename CONTAINER , typename FN >
bool Gambit::ColliderBit::all_of ( const CONTAINER &  c,
const FN &  f 
)
inline

Non-iterator version of std::all_of.

Definition at line 244 of file Utils.hpp.

244  {
245  return std::all_of(std::begin(c), std::end(c), f);
246  }
bool all_of(const CONTAINER &c, const FN &f)
Non-iterator version of std::all_of.
Definition: Utils.hpp:244

◆ any_of()

template<typename CONTAINER , typename FN >
bool Gambit::ColliderBit::any_of ( const CONTAINER &  c,
const FN &  f 
)
inline

Non-iterator version of std::any_of.

Definition at line 250 of file Utils.hpp.

250  {
251  return std::any_of(std::begin(c), std::end(c), f);
252  }
bool any_of(const CONTAINER &c, const FN &f)
Non-iterator version of std::any_of.
Definition: Utils.hpp:250

◆ BFM2SLHA_NN()

void Gambit::ColliderBit::BFM2SLHA_NN ( MixMatrix NN,
double  tanb,
double  sin2thetaW 
)

Converts a neutralino mixing matrix in BFM conventions to SLHA conventions, $\tan\beta$ is as defined in SLHA.

Definition at line 746 of file lep_mssm_xsecs.cpp.

References multiply(), and transpose().

747  {
748  // Define conversion matrix
749  double sinthetaW = sqrt(sin2thetaW);
750  double costhetaW = sqrt(1.-sin2thetaW);
751  double tanv = 1./tanb; // Needed because of convention difference
752  double sinv = sin(atan(tanv));
753  double cosv = cos(atan(tanv));
754  MixMatrix T(4,std::vector<double>(4));
755  T[0][0] = costhetaW; T[0][1] = -sinthetaW;
756  T[1][0] = sinthetaW; T[1][1] = costhetaW;
757  T[2][2] = sinv; T[2][3] = cosv;
758  T[3][2] = -cosv; T[3][3] = sinv;
759  // Multiply N_{SLHA} = N_{BFM} T^T
760  NN = multiply(NN,transpose(T));
761  }
MixMatrix transpose(MixMatrix A)
Helper function to find matrix transpose.
MixMatrix multiply(MixMatrix A, MixMatrix B)
Helper function to multiply matrices.
std::vector< std::vector< double > > MixMatrix
Here is the call graph for this function:

◆ BFM2SLHA_VV()

void Gambit::ColliderBit::BFM2SLHA_VV ( MixMatrix VV)

Converts the chargino mixing matrix V in BFM conventions to SLHA conventions.

Definition at line 764 of file lep_mssm_xsecs.cpp.

References SLHA2BFM_VV().

765  {
766  SLHA2BFM_VV(VV);
767  }
void SLHA2BFM_VV(MixMatrix &VV)
Converts the chargino mixing matrix V in SLHA conventions to BFM conventions.
Here is the call graph for this function:

◆ binIndex()

template<typename NUM1 , typename NUM2 >
size_t Gambit::ColliderBit::binIndex ( NUM1  val,
const std::vector< NUM2 > &  binedges,
bool  allow_overflow = false 
)
inline

< Below/out of histo range

< Above/out of histo range

Definition at line 182 of file Utils.hpp.

182  {
183  if (val < binedges.front()) return -1;
184  if (val >= binedges.back()) return allow_overflow ? int(binedges.size())-1 : -1;
185  return std::distance(binedges.begin(), --std::upper_bound(binedges.begin(), binedges.end(), val));
186  }

◆ calc_combined_LHC_LogLike()

void Gambit::ColliderBit::calc_combined_LHC_LogLike ( double result)

Compute the total likelihood combining all analyses.

Definition at line 690 of file LHC_likelihoods.cpp.

References Gambit::LogTags::debug, debug_prefix(), Gambit::EOM, and Gambit::logger().

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  }
void calc_combined_LHC_LogLike(double &result)
Compute the total likelihood combining all analyses.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
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
str debug_prefix()
Debug prefix string giving thread number.
Here is the call graph for this function:

◆ calc_HB_LEP_LogLike()

void Gambit::ColliderBit::calc_HB_LEP_LogLike ( double result)

Get a LEP chisq from HiggsBounds.

Definition at line 301 of file ColliderBit_Higgs.cpp.

References LOCAL_INFO.

302  {
303  using namespace Pipes::calc_HB_LEP_LogLike;
304 
305  hb_ModelParameters ModelParam = *Dep::HB_ModelParameters;
306 
307  Farray<double, 1,3, 1,3> CS_lep_hjhi_ratio;
308  Farray<double, 1,3, 1,3> BR_hjhihi;
309  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
310  {
311  CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
312  BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
313  }
314 
315  BEreq::HiggsBounds_neutral_input_part(&ModelParam.Mh[0], &ModelParam.hGammaTot[0], &ModelParam.CP[0],
316  &ModelParam.CS_lep_hjZ_ratio[0], &ModelParam.CS_lep_bbhj_ratio[0],
317  &ModelParam.CS_lep_tautauhj_ratio[0], CS_lep_hjhi_ratio,
318  &ModelParam.CS_gg_hj_ratio[0], &ModelParam.CS_bb_hj_ratio[0],
319  &ModelParam.CS_bg_hjb_ratio[0], &ModelParam.CS_ud_hjWp_ratio[0],
320  &ModelParam.CS_cs_hjWp_ratio[0], &ModelParam.CS_ud_hjWm_ratio[0],
321  &ModelParam.CS_cs_hjWm_ratio[0], &ModelParam.CS_gg_hjZ_ratio[0],
322  &ModelParam.CS_dd_hjZ_ratio[0], &ModelParam.CS_uu_hjZ_ratio[0],
323  &ModelParam.CS_ss_hjZ_ratio[0], &ModelParam.CS_cc_hjZ_ratio[0],
324  &ModelParam.CS_bb_hjZ_ratio[0], &ModelParam.CS_tev_vbf_ratio[0],
325  &ModelParam.CS_tev_tthj_ratio[0], &ModelParam.CS_lhc7_vbf_ratio[0],
326  &ModelParam.CS_lhc7_tthj_ratio[0], &ModelParam.CS_lhc8_vbf_ratio[0],
327  &ModelParam.CS_lhc8_tthj_ratio[0], &ModelParam.BR_hjss[0],
328  &ModelParam.BR_hjcc[0], &ModelParam.BR_hjbb[0],
329  &ModelParam.BR_hjmumu[0], &ModelParam.BR_hjtautau[0],
330  &ModelParam.BR_hjWW[0], &ModelParam.BR_hjZZ[0],
331  &ModelParam.BR_hjZga[0], &ModelParam.BR_hjgaga[0],
332  &ModelParam.BR_hjgg[0], &ModelParam.BR_hjinvisible[0], BR_hjhihi);
333 
334  BEreq::HiggsBounds_charged_input(&ModelParam.MHplus[0], &ModelParam.HpGammaTot[0], &ModelParam.CS_lep_HpjHmi_ratio[0],
335  &ModelParam.BR_tWpb, &ModelParam.BR_tHpjb[0], &ModelParam.BR_Hpjcs[0],
336  &ModelParam.BR_Hpjcb[0], &ModelParam.BR_Hptaunu[0]);
337 
338  BEreq::HiggsBounds_set_mass_uncertainties(&ModelParam.deltaMh[0],&ModelParam.deltaMHplus[0]);
339 
340  // run Higgs bounds 'classic'
341  double obsratio;
342  int HBresult, chan, ncombined;
343  BEreq::run_HiggsBounds_classic(HBresult,chan,obsratio,ncombined);
344 
345  // extract the LEP chisq
346  double chisq_withouttheory,chisq_withtheory;
347  int chan2;
348  double theor_unc = 1.5; // theory uncertainty
349  BEreq::HB_calc_stats(theor_unc,chisq_withouttheory,chisq_withtheory,chan2);
350 
351  // Catch HiggsBound's error value, chisq = -999
352  if( fabs(chisq_withouttheory - (-999.)) < 1e-6)
353  {
354  ColliderBit_warning().raise(LOCAL_INFO, "Got chisq=-999 from HB_calc_stats in HiggsBounds, indicating a cross-section outside tabulated range. Will use chisq=0.");
355  chisq_withouttheory = 0.0;
356  }
357 
358  result = -0.5*chisq_withouttheory;
359  }
#define LOCAL_INFO
Definition: local_info.hpp:34
void calc_HB_LEP_LogLike(double &result)
Get a LEP chisq from HiggsBounds.

◆ calc_HS_LHC_LogLike()

void Gambit::ColliderBit::calc_HS_LHC_LogLike ( double result)

Get an LHC chisq from HiggsSignals.

Definition at line 362 of file ColliderBit_Higgs.cpp.

References combine_hdf5::f.

363  {
364  using namespace Pipes::calc_HS_LHC_LogLike;
365 
366  hb_ModelParameters ModelParam = *Dep::HB_ModelParameters;
367 
368  Farray<double, 1,3, 1,3> CS_lep_hjhi_ratio;
369  Farray<double, 1,3, 1,3> BR_hjhihi;
370  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
371  {
372  CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
373  BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
374  }
375 
376  BEreq::HiggsBounds_neutral_input_part_HS(&ModelParam.Mh[0], &ModelParam.hGammaTot[0], &ModelParam.CP[0],
377  &ModelParam.CS_lep_hjZ_ratio[0], &ModelParam.CS_lep_bbhj_ratio[0],
378  &ModelParam.CS_lep_tautauhj_ratio[0], CS_lep_hjhi_ratio,
379  &ModelParam.CS_gg_hj_ratio[0], &ModelParam.CS_bb_hj_ratio[0],
380  &ModelParam.CS_bg_hjb_ratio[0], &ModelParam.CS_ud_hjWp_ratio[0],
381  &ModelParam.CS_cs_hjWp_ratio[0], &ModelParam.CS_ud_hjWm_ratio[0],
382  &ModelParam.CS_cs_hjWm_ratio[0], &ModelParam.CS_gg_hjZ_ratio[0],
383  &ModelParam.CS_dd_hjZ_ratio[0], &ModelParam.CS_uu_hjZ_ratio[0],
384  &ModelParam.CS_ss_hjZ_ratio[0], &ModelParam.CS_cc_hjZ_ratio[0],
385  &ModelParam.CS_bb_hjZ_ratio[0], &ModelParam.CS_tev_vbf_ratio[0],
386  &ModelParam.CS_tev_tthj_ratio[0], &ModelParam.CS_lhc7_vbf_ratio[0],
387  &ModelParam.CS_lhc7_tthj_ratio[0], &ModelParam.CS_lhc8_vbf_ratio[0],
388  &ModelParam.CS_lhc8_tthj_ratio[0], &ModelParam.BR_hjss[0],
389  &ModelParam.BR_hjcc[0], &ModelParam.BR_hjbb[0],
390  &ModelParam.BR_hjmumu[0], &ModelParam.BR_hjtautau[0],
391  &ModelParam.BR_hjWW[0], &ModelParam.BR_hjZZ[0],
392  &ModelParam.BR_hjZga[0], &ModelParam.BR_hjgaga[0],
393  &ModelParam.BR_hjgg[0], &ModelParam.BR_hjinvisible[0], BR_hjhihi);
394 
395  BEreq::HiggsBounds_charged_input_HS(&ModelParam.MHplus[0], &ModelParam.HpGammaTot[0], &ModelParam.CS_lep_HpjHmi_ratio[0],
396  &ModelParam.BR_tWpb, &ModelParam.BR_tHpjb[0], &ModelParam.BR_Hpjcs[0],
397  &ModelParam.BR_Hpjcb[0], &ModelParam.BR_Hptaunu[0]);
398 
399  BEreq::HiggsSignals_neutral_input_MassUncertainty(&ModelParam.deltaMh[0]);
400 
401  // add uncertainties to cross-sections and branching ratios
402  // double dCS[5] = {0.,0.,0.,0.,0.};
403  // double dBR[5] = {0.,0.,0.,0.,0.};
404  // BEreq::setup_rate_uncertainties(dCS,dBR);
405 
406  // run HiggsSignals
407  int mode = 1; // 1- peak-centered chi2 method (recommended)
408  double csqmu, csqmh, csqtot, Pvalue;
409  int nobs;
410  BEreq::run_HiggsSignals(mode, csqmu, csqmh, csqtot, nobs, Pvalue);
411 
412  result = -0.5*csqtot;
413 
414  #ifdef COLLIDERBIT_DEBUG
415  std::ofstream f;
416  f.open ("HB_ModelParameters_contents.dat");
417  f<<"LHC log-likleihood";
418  for (int i = 0; i < 3; i++) f<<
419  " higgs index" <<
420  " "<<i<<":CP"<<
421  " "<<i<<":Mh"<<
422  " "<<i<<":hGammaTot"<<
423  " "<<i<<":CS_lep_hjZ_ratio"<<
424  " "<<i<<":CS_tev_vbf_ratio"<<
425  " "<<i<<":CS_lep_bbhj_ratio"<<
426  " "<<i<<":CS_lep_tautauhj_ratio"<<
427  " "<<i<<":CS_gg_hj_ratio"<<
428  " "<<i<<":CS_tev_tthj_ratio"<<
429  " "<<i<<":CS_lhc7_tthj_ratio"<<
430  " "<<i<<":CS_lhc8_tthj_ratio"<<
431  " "<<i<<":CS_lep_hjhi_ratio[0]"<<
432  " "<<i<<":CS_lep_hjhi_ratio[1]"<<
433  " "<<i<<":CS_lep_hjhi_ratio[2]"<<
434  " "<<i<<":BR_ss"<<
435  " "<<i<<":BR_cc"<<
436  " "<<i<<":BR_bb"<<
437  " "<<i<<":BR_mumu"<<
438  " "<<i<<":BR_tautau"<<
439  " "<<i<<":BR_WW"<<
440  " "<<i<<":BR_ZZ"<<
441  " "<<i<<":BR_Zga"<<
442  " "<<i<<":BR_gamgam"<<
443  " "<<i<<":BR_gg"<<
444  " "<<i<<":BR_invisible"<<
445  " "<<i<<":BR_hihi[0]"<<
446  " "<<i<<":BR_hihi[1]"<<
447  " "<<i<<":BR_hihi[2]";
448  f<<
449  " higgs index" <<
450  " "<<4<<"MHplus"<<
451  " "<<4<<":HpGammaTot"<<
452  " "<<4<<":CS_lep_HpjHmi_ratio"<<
453  " "<<4<<":BR_H+->cs"<<
454  " "<<4<<":BR_H+->cb"<<
455  " "<<4<<":BR_H+->taunu"<<
456  " "<<4<<":BR_t->W+b"<<
457  " "<<4<<":BR_t->H+b";
458  f << endl << std::setw(18) << result;
459  const int w = 24;
460  for (int i = 0; i < 3; i++)
461  {
462  f << std::setw(w) << i << std::setw(w) <<
463  ModelParam.CP[i] << std::setw(w) <<
464  ModelParam.Mh[i] << std::setw(w) <<
465  ModelParam.hGammaTot[i] << std::setw(w) <<
466  ModelParam.CS_lep_hjZ_ratio[i] << std::setw(w) <<
467  ModelParam.CS_tev_vbf_ratio[i] << std::setw(w) <<
468  ModelParam.CS_lep_bbhj_ratio[i] << std::setw(w) <<
469  ModelParam.CS_lep_tautauhj_ratio[i] << std::setw(w) <<
470  ModelParam.CS_gg_hj_ratio[i] << std::setw(w) <<
471  ModelParam.CS_tev_tthj_ratio[i] << std::setw(w) <<
472  ModelParam.CS_lhc7_tthj_ratio[i] << std::setw(w) <<
473  ModelParam.CS_lhc8_tthj_ratio[i];
474  for (int j = 0; j < 3; j++) f << std::setw(w) << ModelParam.CS_lep_hjhi_ratio[i][j];
475  f << std::setw(w) <<
476  ModelParam.BR_hjss[i] << std::setw(w) <<
477  ModelParam.BR_hjcc[i] << std::setw(w) <<
478  ModelParam.BR_hjbb[i] << std::setw(w) <<
479  ModelParam.BR_hjmumu[i] << std::setw(w) <<
480  ModelParam.BR_hjtautau[i] << std::setw(w) <<
481  ModelParam.BR_hjWW[i] << std::setw(w) <<
482  ModelParam.BR_hjZZ[i] << std::setw(w) <<
483  ModelParam.BR_hjZga[i] << std::setw(w) <<
484  ModelParam.BR_hjgaga[i] << std::setw(w) <<
485  ModelParam.BR_hjgg[i] << std::setw(w) <<
486  ModelParam.BR_hjinvisible[i];
487  for (int j = 0; j < 3; j++) f << std::setw(w) << ModelParam.BR_hjhihi[i][j];
488  }
489  f << std::setw(w) << 4 << std::setw(w) <<
490  ModelParam.MHplus[0] << std::setw(w) <<
491  ModelParam.HpGammaTot[0] << std::setw(w) <<
492  ModelParam.CS_lep_HpjHmi_ratio[0] << std::setw(w) <<
493  ModelParam.BR_Hpjcs[0] << std::setw(w) <<
494  ModelParam.BR_Hpjcb[0] << std::setw(w) <<
495  ModelParam.BR_Hptaunu[0] << std::setw(w) <<
496  ModelParam.BR_tWpb << std::setw(w) <<
497  ModelParam.BR_tHpjb[0];
498  f.close();
499  #endif
500 
501  }
void calc_HS_LHC_LogLike(double &result)
Get an LHC chisq from HiggsSignals.

◆ calc_LHC_LogLikes()

void Gambit::ColliderBit::calc_LHC_LogLikes ( map_str_AnalysisLogLikes result)

Loop over all analyses and fill a map of AnalysisLogLikes objects.

If no events have been generated (xsec veto) or too many events have failed, short-circut the loop and return delta log-likelihood = 0 for every SR in each analysis.

If (simplified) SR-correlation info is available, so use the covariance matrix to construct composite marginalised likelihood Despite initial thoughts, we can't just do independent LL calculations in a rotated basis, but have to sample from the covariance matrix.

Note
This means we can't use the nulike LL functions, which operate in 1D only. Also, log-normal sampling in the diagonal basis is not helpful, since the rotation will re-generate negative rates.
Note
How to correct negative rates? Discard (scales badly), set to epsilon (= discontinuous & unphysical pdf), transform to log-space (distorts the pdf quite badly), or something else (skew term)? We're using the "set to epsilon" version for now. Ben: I would vote for 'discard'. It can't be that inefficient, surely?

Definition at line 93 of file LHC_likelihoods.cpp.

References Gambit::ColliderBit::AnalysisData::analysis_name, Gambit::ColliderBit::SignalRegionData::background_sys, debug_prefix(), double, int, Gambit::invalid_point(), lnlike_marg_poisson_lognormal_error, Gambit::ColliderBit::SignalRegionData::n_background, Gambit::ColliderBit::SignalRegionData::n_observed, Gambit::ColliderBit::SignalRegionData::n_signal, Gambit::ColliderBit::SignalRegionData::n_signal_at_lumi, Gambit::invalid_point_exception::raise(), Gambit::Random::rng(), Gambit::ColliderBit::SignalRegionData::signal_sys, Gambit::ColliderBit::AnalysisData::size(), Gambit::ColliderBit::SignalRegionData::sr_label, and Gambit::ColliderBit::AnalysisData::srcov.

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  }
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
void calc_LHC_LogLikes(map_str_AnalysisLogLikes &result)
Loop over all analyses and fill a map of AnalysisLogLikes objects.
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
DS5_MSPCTM DS_INTDOF int
invalid_point_exception & invalid_point()
Invalid point exceptions.
str debug_prefix()
Debug prefix string giving thread number.
AnalysisDataPointers lnlike_marg_poisson_lognormal_error
Here is the call graph for this function:

◆ calc_LHC_signals()

void Gambit::ColliderBit::calc_LHC_signals ( map_str_dbl result)

Loop over all analyses and fill a map of predicted counts.

Definition at line 55 of file LHC_likelihoods.cpp.

References Gambit::ColliderBit::AnalysisData::analysis_name, Gambit::LogTags::debug, Gambit::EOM, Gambit::logger(), Gambit::ColliderBit::SignalRegionData::n_signal, Gambit::ColliderBit::SignalRegionData::n_signal_at_lumi, Gambit::ColliderBit::SignalRegionData::signal_sys, Gambit::ColliderBit::AnalysisData::size(), and Gambit::ColliderBit::SignalRegionData::sr_label.

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  }
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
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
void calc_LHC_signals(map_str_dbl &result)
Loop over all analyses and fill a map of predicted counts.
Here is the call graph for this function:

◆ calcMT()

double Gambit::ColliderBit::calcMT ( HEPUtils::P4  jetMom,
HEPUtils::P4  metMom 
)

Definition at line 41 of file Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp.

References Gambit::ColliderBit::ATLAS::applyElectronEff(), Gambit::ColliderBit::ATLAS::applyMuonEff(), b, colouring::combine(), DEFINE_ANALYSIS_FACTORY, get_jets(), mt2_bisect::mt2::get_mt2(), has_tag(), int, mt, run(), mt2_bisect::mt2::set_mn(), mt2_bisect::mt2::set_momenta(), sortByPT13(), and sortByPT13_sharedptr().

41  {
42 
43  //std::cout << "metMom.px() " << metMom.px() << " jetMom PT " << jetMom.pT() << std::endl;
44 
45  double met=sqrt(metMom.px()*metMom.px()+metMom.py()*metMom.py());
46  double dphi = metMom.deltaPhi(jetMom);
47  double mt=sqrt(2*jetMom.pT()*met*(1-std::cos(dphi)));
48 
49  return mt;
50 
51  }
const double mt
Definition: topness.h:39
Here is the call graph for this function:

◆ calcMT_1l()

double Gambit::ColliderBit::calcMT_1l ( HEPUtils::P4  jetMom,
HEPUtils::P4  metMom 
)

◆ cmpJetsByPt()

bool Gambit::ColliderBit::cmpJetsByPt ( const HEPUtils::Jet *  jet1,
const HEPUtils::Jet *  jet2 
)
inline

Comparison function to give a jets sorting order decreasing by pT.

Definition at line 292 of file Utils.hpp.

References sortBy(), and sortByPt().

292 { return jet1->pT() > jet2->pT(); }
Here is the call graph for this function:

◆ cmpParticlesByPt()

bool Gambit::ColliderBit::cmpParticlesByPt ( const HEPUtils::Particle *  lep1,
const HEPUtils::Particle *  lep2 
)
inline

Comparison function to give a particles sorting order decreasing by pT.

Definition at line 280 of file Utils.hpp.

Referenced by sortByPt().

280 { return lep1->pT() > lep2->pT(); }
Here is the caller graph for this function:

◆ CollectAnalyses()

void Gambit::ColliderBit::CollectAnalyses ( AnalysisDataPointers result)

Loop over all analyses and collect them in one place.

Definition at line 270 of file ColliderBit_eventloop.cpp.

References debug_prefix(), and LOCAL_INFO.

271  {
272  using namespace Pipes::CollectAnalyses;
273  static bool first = true;
274 
275  // Start with an empty vector
276  result.clear();
277 
278  #ifdef COLLIDERBIT_DEBUG
279  cout << debug_prefix() << "CollectAnalyses: Dep::ATLASAnalysisNumbers->size() = " << Dep::ATLASAnalysisNumbers->size() << endl;
280  cout << debug_prefix() << "CollectAnalyses: Dep::CMSAnalysisNumbers->size() = " << Dep::CMSAnalysisNumbers->size() << endl;
281  cout << debug_prefix() << "CollectAnalyses: Dep::IdentityAnalysisNumbers->size() = " << Dep::IdentityAnalysisNumbers->size() << endl;
282  #endif
283 
284  // Add results
285  if (Dep::ATLASAnalysisNumbers->size() != 0) result.insert(result.end(), Dep::ATLASAnalysisNumbers->begin(), Dep::ATLASAnalysisNumbers->end());
286  if (Dep::CMSAnalysisNumbers->size() != 0) result.insert(result.end(), Dep::CMSAnalysisNumbers->begin(), Dep::CMSAnalysisNumbers->end());
287  if (Dep::IdentityAnalysisNumbers->size() != 0) result.insert(result.end(), Dep::IdentityAnalysisNumbers->begin(), Dep::IdentityAnalysisNumbers->end());
288 
289  // When first called, check that all analyses contain at least one signal region.
290  if (first)
291  {
292  // Loop over all AnalysisData pointers
293  for (auto& adp : result)
294  {
295  if (adp->size() == 0)
296  {
297  str errmsg;
298  errmsg = "The analysis " + adp->analysis_name + " has no signal regions.";
299  ColliderBit_error().raise(LOCAL_INFO, errmsg);
300  }
301  }
302  first = false;
303  }
304 
305 
306  // #ifdef COLLIDERBIT_DEBUG
307  // cout << debug_prefix() << "CollectAnalyses: Current size of 'result': " << result.size() << endl;
308  // if (result.size() > 0)
309  // {
310  // cout << debug_prefix() << "CollectAnalyses: Will loop through 'result'..." << endl;
311  // for (auto& adp : result)
312  // {
313  // cout << debug_prefix() << "CollectAnalyses: 'result' contains AnalysisData pointer to " << adp << endl;
314  // cout << debug_prefix() << "CollectAnalyses: -- Will now loop over all signal regions in " << adp << endl;
315  // for (auto& sr : adp->srdata)
316  // {
317  // cout << debug_prefix() << "CollectAnalyses: -- " << adp << " contains signal region: " << sr.sr_label << ", n_signal = " << sr.n_signal << ", n_signal_at_lumi = " << n_signal_at_lumi << endl;
318  // }
319  // cout << debug_prefix() << "CollectAnalyses: -- Done looping over signal regions in " << adp << endl;
320  // }
321  // cout << debug_prefix() << "CollectAnalyses: ...Done looping through 'result'." << endl;
322  // }
323  // #endif
324  }
#define LOCAL_INFO
Definition: local_info.hpp:34
void CollectAnalyses(AnalysisDataPointers &result)
Loop over all analyses and collect them in one place.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
str debug_prefix()
Debug prefix string giving thread number.
Here is the call graph for this function:

◆ debug_prefix()

str Gambit::ColliderBit::debug_prefix ( )
inline

Debug prefix string giving thread number.

Definition at line 64 of file ColliderBit_eventloop.hpp.

Referenced by calc_combined_LHC_LogLike(), calc_LHC_LogLikes(), CollectAnalyses(), generateEventColliderPythia(), getColliderPythia(), getMCxsec(), operateLHCLoop(), runAnalyses(), and smearEvent().

65  {
66  std::stringstream ss;
67  ss << "DEBUG: OMP thread " << omp_get_thread_num() << ": ";
68  return ss.str();
69  }
Here is the caller graph for this function:

◆ DEFINE_ANALYSIS_FACTORY()

Gambit::ColliderBit::DEFINE_ANALYSIS_FACTORY ( ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb  )

Referenced by Gambit::ColliderBit::AnalysisData::pythonize_me(), and sortByPT_2lep().

Here is the caller graph for this function:

◆ FH_HiggsProd()

void Gambit::ColliderBit::FH_HiggsProd ( fh_HiggsProd &  result)

Higgs production cross-sections from FeynHiggs.

Definition at line 504 of file ColliderBit_Higgs.cpp.

References Gambit::LogTags::err, Gambit::invalid_point(), and Gambit::invalid_point_exception::raise().

505  {
506  using namespace Pipes::FH_HiggsProd;
507 
508  Farray<fh_real, 1,52> prodxs;
509 
510  fh_HiggsProd HiggsProd;
511  int error;
512  fh_real sqrts;
513 
514  // Tevatron
515  sqrts = 2.;
516  error = 1;
517  BEreq::FHHiggsProd(error, sqrts, prodxs);
518  if (error != 0)
519  {
520  std::ostringstream err;
521  err << "BEreq::FHHiggsProd raised error flag for Tevatron: " << error << ".";
522  invalid_point().raise(err.str());
523  }
524  for(int i = 0; i < 52; i++) HiggsProd.prodxs_Tev[i] = prodxs(i+1);
525  // LHC7
526  sqrts = 7.;
527  error = 1;
528  BEreq::FHHiggsProd(error, sqrts, prodxs);
529  if (error != 0)
530  {
531  std::ostringstream err;
532  err << "BEreq::FHHiggsProd raised error flag for LHC7: " << error << ".";
533  invalid_point().raise(err.str());
534  }
535  for(int i = 0; i < 52; i++) HiggsProd.prodxs_LHC7[i] = prodxs(i+1);
536  // LHC8
537  sqrts = 8.;
538  error = 1;
539  BEreq::FHHiggsProd(error, sqrts, prodxs);
540  if (error != 0)
541  {
542  std::ostringstream err;
543  err << "BEreq::FHHiggsProd raised error flag for LHC8: " << error << ".";
544  invalid_point().raise(err.str());
545  }
546  for(int i = 0; i < 52; i++) HiggsProd.prodxs_LHC8[i] = prodxs(i+1);
547 
548  // The ttbar production cross-sections for the (BSM,SM) model can be found at (prodxs_X[h+27], prodxs_X[h+30]),
549  // where h is the higgs index (0 = h0_1, 1 = h0_2, 2 = A0) and X is one of Tev, LHC7 or LHC8.
550  result = HiggsProd;
551 
552  }
void FH_HiggsProd(fh_HiggsProd &result)
Higgs production cross-sections from FeynHiggs.
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
invalid_point_exception & invalid_point()
Invalid point exceptions.
Here is the call graph for this function:

◆ filter_reject() [1/2]

ParticlePtrs Gambit::ColliderBit::filter_reject ( const ParticlePtrs particles,
std::function< bool(Particle *)>  rejfn,
bool  do_delete = true 
)
inline

Filter a supplied particle vector by rejecting those which fail a supplied cut.

Definition at line 71 of file Utils.hpp.

References ifilter_reject().

Referenced by filter_select().

72  {
73  ParticlePtrs rtn = particles;
74  ifilter_reject(rtn, rejfn, do_delete);
75  return rtn;
76  }
std::vector< HEPUtils::Particle * > ParticlePtrs
Typedef for a vector of Particle pointers.
Definition: Utils.hpp:31
void ifilter_reject(JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
In-place filter a supplied jet vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:92
Here is the call graph for this function:
Here is the caller graph for this function:

◆ filter_reject() [2/2]

JetPtrs Gambit::ColliderBit::filter_reject ( const JetPtrs jets,
std::function< bool(Jet *)>  rejfn,
bool  do_delete = true 
)
inline

Filter a supplied particle vector by rejecting those which fail a supplied cut.

Definition at line 110 of file Utils.hpp.

References ifilter_reject().

111  {
112  JetPtrs rtn = jets;
113  ifilter_reject(rtn, rejfn, do_delete);
114  return rtn;
115  }
void ifilter_reject(JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
In-place filter a supplied jet vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:92
std::vector< HEPUtils::Jet * > JetPtrs
Typedef for a vector of Jet pointers.
Definition: Utils.hpp:37
Here is the call graph for this function:

◆ filter_select() [1/2]

ParticlePtrs Gambit::ColliderBit::filter_select ( const ParticlePtrs particles,
std::function< bool(Particle *)>  selfn,
bool  do_delete = true 
)
inline

Filter a supplied particle vector by keeping those which pass a supplied cut.

Definition at line 79 of file Utils.hpp.

References filter_reject().

80  {
81  return filter_reject(particles, [&](Particle* p) { return !selfn(p); }, do_delete);
82  }
JetPtrs filter_reject(const JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
Filter a supplied particle vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:110
Here is the call graph for this function:

◆ filter_select() [2/2]

JetPtrs Gambit::ColliderBit::filter_select ( const JetPtrs jets,
std::function< bool(Jet *)>  selfn,
bool  do_delete = true 
)
inline

Filter a supplied particle vector by keeping those which pass a supplied cut.

Definition at line 118 of file Utils.hpp.

References filter_reject(), and random_bool().

119  {
120  return filter_reject(jets, [&](Jet* j) { return !selfn(j); }, do_delete);
121  }
JetPtrs filter_reject(const JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
Filter a supplied particle vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:110
Here is the call graph for this function:

◆ filtereff() [1/2]

void Gambit::ColliderBit::filtereff ( std::vector< HEPUtils::Particle *> &  particles,
double  eff,
bool  do_delete = false 
)

Utility function for filtering a supplied particle vector by sampling wrt an efficiency scalar.

Definition at line 16 of file Utils.cpp.

References random_bool().

Referenced by Gambit::ColliderBit::CMS::applyTauEfficiency(), Gambit::ColliderBit::ATLAS::applyTauEfficiencyR1(), and random_bool().

16  {
17  if (particles.empty()) return;
18  auto keptParticlesEnd = std::remove_if(particles.begin(), particles.end(),
19  [&](HEPUtils::Particle* p) {
20  const bool rm = !random_bool(eff);
21  if (do_delete && rm) delete p;
22  return rm;
23  } );
24  particles.erase(keptParticlesEnd, particles.end());
25  }
bool random_bool(double eff)
Return a random true/false at a success rate given by a number.
Definition: Utils.cpp:10
Here is the call graph for this function:
Here is the caller graph for this function:

◆ filtereff() [2/2]

void Gambit::ColliderBit::filtereff ( std::vector< HEPUtils::Particle *> &  particles,
std::function< double(HEPUtils::Particle *)>  eff_fn,
bool  do_delete = false 
)

Utility function for filtering a supplied particle vector by sampling an efficiency returned by a provided function object.

Utility function for filtering a supplied particle vector by sampling wrt a binned 1D efficiency map in pT.

Definition at line 29 of file Utils.cpp.

References random_bool().

29  {
30  if (particles.empty()) return;
31  auto keptParticlesEnd = std::remove_if(particles.begin(), particles.end(),
32  [&](HEPUtils::Particle* p) {
33  const double eff = eff_fn(p);
34  const bool rm = !random_bool(eff);
35  if (do_delete && rm) delete p;
36  return rm;
37  } );
38  particles.erase(keptParticlesEnd, particles.end());
39  }
bool random_bool(double eff)
Return a random true/false at a success rate given by a number.
Definition: Utils.cpp:10
Here is the call graph for this function:

◆ filtereff_etapt()

void Gambit::ColliderBit::filtereff_etapt ( std::vector< HEPUtils::Particle *> &  particles,
const HEPUtils::BinnedFn2D< double > &  eff_etapt,
bool  do_delete = false 
)

Utility function for filtering a supplied particle vector by sampling wrt a binned 2D efficiency map in |eta| and pT.

Definition at line 56 of file Utils.cpp.

References random_bool().

Referenced by Gambit::ColliderBit::ATLAS::applyElectronEff(), Gambit::ColliderBit::CMS::applyElectronEff(), Gambit::ColliderBit::CMS::applyElectronTrackingEff(), Gambit::ColliderBit::ATLAS::applyMuonEff(), Gambit::ColliderBit::ATLAS::applyMuonEffR2(), Gambit::ColliderBit::CMS::applyMuonTrackEff(), Gambit::ColliderBit::ATLAS::applyPhotonEfficiencyR2(), and random_bool().

56  {
57  if (particles.empty()) return;
58  auto keptParticlesEnd = std::remove_if(particles.begin(), particles.end(),
59  [&](const HEPUtils::Particle* p) {
60  const bool rm = !random_bool(eff_etapt, p->abseta(), p->pT());
61  if (do_delete && rm) delete p;
62  return rm;
63  } );
64  particles.erase(keptParticlesEnd, particles.end());
65  }
bool random_bool(double eff)
Return a random true/false at a success rate given by a number.
Definition: Utils.cpp:10
Here is the call graph for this function:
Here is the caller graph for this function:

◆ filtereff_pt()

void Gambit::ColliderBit::filtereff_pt ( std::vector< HEPUtils::Particle *> &  particles,
const HEPUtils::BinnedFn1D< double > &  eff_pt,
bool  do_delete = false 
)

Utility function for filtering a supplied particle vector by sampling wrt a binned 1D efficiency map in pT.

Definition at line 43 of file Utils.cpp.

References random_bool().

Referenced by Gambit::ColliderBit::ATLAS::applyTauEfficiencyR2(), and random_bool().

43  {
44  if (particles.empty()) return;
45  auto keptParticlesEnd = std::remove_if(particles.begin(), particles.end(),
46  [&](const HEPUtils::Particle* p) {
47  const bool rm = !random_bool(eff_pt, p->pT());
48  if (do_delete && rm) delete p;
49  return rm;
50  } );
51  particles.erase(keptParticlesEnd, particles.end());
52  }
bool random_bool(double eff)
Return a random true/false at a success rate given by a number.
Definition: Utils.cpp:10
Here is the call graph for this function:
Here is the caller graph for this function:

◆ fromBottom()

template<typename EventT >
bool Gambit::ColliderBit::fromBottom ( int  n,
const EventT &  evt 
)
inline

Definition at line 57 of file Py8Utils.hpp.

58  {
59  // Root particle is invalid
60  if (n == 0) return false;
61  const auto& p = evt[n];
62  if (abs(p.id()) == 5 || MCUtils::PID::hasBottom(p.id())) return true;
64  if (p.isParton()) return false; // stop the walking at hadron level
65  for (int m : p.motherList()) {
66  if (fromBottom(m, evt)) return true;
67  }
68  return false;
69  }
bool fromBottom(int n, const EventT &evt)
Definition: Py8Utils.hpp:57

◆ fromHadron()

template<typename EventT >
bool Gambit::ColliderBit::fromHadron ( int  n,
const EventT &  evt 
)
inline

Definition at line 90 of file Py8Utils.hpp.

Referenced by Gambit::ColliderBit::BuckFast< EventT >::convertParticleEvent().

91  {
92  // Root particle is invalid
93  if (n == 0) return false;
94  const auto& p = evt[n];
95  if (p.isHadron()) return true;
96  if (p.isParton()) return false; // stop the walking at the end of the hadron level
97  for (int m : p.motherList()) {
98  if (fromHadron(m, evt)) return true;
99  }
100  return false;
101  }
bool fromHadron(int n, const EventT &evt)
Definition: Py8Utils.hpp:90
Here is the caller graph for this function:

◆ fromTau()

template<typename EventT >
bool Gambit::ColliderBit::fromTau ( int  n,
const EventT &  evt 
)
inline

Definition at line 74 of file Py8Utils.hpp.

75  {
76  // Root particle is invalid
77  if (n == 0) return false;
78  const auto& p = evt[n];
79  if (abs(p.id()) == 15) return true;
80  if (p.isParton()) return false; // stop the walking at the end of the hadron level
81  for (int m : p.motherList()) {
82  if (fromTau(m, evt)) return true;
83  }
84  return false;
85  }
bool fromTau(int n, const EventT &evt)
Definition: Py8Utils.hpp:74

◆ generateEventColliderPythia()

template<typename PythiaT , typename EventT >
void Gambit::ColliderBit::generateEventColliderPythia ( EventT &  result,
const MCLoopInfo RunMC,
const ColliderPythia< PythiaT, EventT > &  HardScatteringSim,
const int  iteration,
void(*)()  wrapup 
)

Generate a hard scattering event with Pythia.

Definition at line 51 of file generateEventColliderPythia.hpp.

References BASE_INIT, Gambit::ColliderBit::MCLoopInfo::current_maxFailedEvents(), Gambit::LogTags::debug, debug_prefix(), Gambit::EOM, LOCAL_INFO, Gambit::logger(), Gambit::ColliderBit::ColliderPythia< PythiaT, EventT >::nextEvent(), Gambit::piped_warnings, and Gambit::Piped_exceptions::request().

56  {
57  static int nFailedEvents;
58 
59  // If the event loop has not been entered yet, reset the counter for the number of failed events
60  if (iteration == BASE_INIT)
61  {
62  // Although BASE_INIT should never be called in parallel, we do this in omp_critical just in case.
63  #pragma omp critical (pythia_event_failure)
64  {
65  nFailedEvents = 0;
66  }
67  return;
68  }
69 
70  // If in any other special iteration, do nothing
71  if (iteration < BASE_INIT) return;
72 
73  // Reset the event
74  result.clear();
75 
76  // Attempt (possibly repeatedly) to generate an event
77  while(nFailedEvents <= RunMC.current_maxFailedEvents())
78  {
79  try
80  {
81  HardScatteringSim.nextEvent(result);
82  break;
83  }
84  catch (typename ColliderPythia<PythiaT,EventT>::EventGenerationError& e)
85  {
86  #ifdef COLLIDERBIT_DEBUG
87  cout << debug_prefix() << "ColliderPythia::EventGenerationError caught in generatePythia8Event. Check the ColliderBit log for event details." << endl;
88  #endif
89  #pragma omp critical (pythia_event_failure)
90  {
91  // Update global counter
92  nFailedEvents += 1;
93  // Store Pythia event record in the logs
94  std::stringstream ss;
95  result.list(ss, 1);
96  logger() << LogTags::debug << "ColliderPythia::EventGenerationError error caught in generatePythia8Event. Pythia record for event that failed:\n" << ss.str() << EOM;
97  }
98  }
99  }
100  // Wrap up event loop if too many events fail.
101  if(nFailedEvents > RunMC.current_maxFailedEvents())
102  {
103  piped_warnings.request(LOCAL_INFO,"exceeded maxFailedEvents");
104  wrapup();
105  return;
106  }
107  }
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:527
#define LOCAL_INFO
Definition: local_info.hpp:34
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
str debug_prefix()
Debug prefix string giving thread number.
Here is the call graph for this function:

◆ GET_ANALYSIS_CONTAINER()

Gambit::ColliderBit::GET_ANALYSIS_CONTAINER ( getATLASAnalysisContainer  ,
ATLAS   
)

◆ GET_BUCKFAST_AS_BASE_DETECTOR() [1/4]

Pythia_EM_default::Pythia8::Event Gambit::ColliderBit::GET_BUCKFAST_AS_BASE_DETECTOR ( getBuckFastATLASPythia_EM  ,
Pythia_EM_default::Pythia8::Event  ,
ATLAS   
)

◆ GET_BUCKFAST_AS_BASE_DETECTOR() [2/4]

Pythia_default::Pythia8::Event Gambit::ColliderBit::GET_BUCKFAST_AS_BASE_DETECTOR ( getBuckFastATLASPythia  ,
Pythia_default::Pythia8::Event  ,
ATLAS   
)

◆ GET_BUCKFAST_AS_BASE_DETECTOR() [3/4]

Pythia_EM_default::Pythia8::Event CMS Gambit::ColliderBit::GET_BUCKFAST_AS_BASE_DETECTOR ( getBuckFastIdentityPythia_EM  ,
Pythia_EM_default::Pythia8::Event  ,
Identity   
)

◆ GET_BUCKFAST_AS_BASE_DETECTOR() [4/4]

Pythia_default::Pythia8::Event CMS Gambit::ColliderBit::GET_BUCKFAST_AS_BASE_DETECTOR ( getBuckFastIdentityPythia  ,
Pythia_default::Pythia8::Event  ,
Identity   
)

◆ get_jets()

template<typename MOM >
std::vector<std::shared_ptr<HEPUtils::Jet> > Gambit::ColliderBit::get_jets ( const std::vector< MOM *> &  moms,
double  R,
double  ptmin = 0*GeV,
FJNS::JetAlgorithm  alg = FJNS::antikt_algorithm 
)
inline

Run jet clustering from any P4-compatible momentum type.

Definition at line 205 of file Utils.hpp.

References mk_p4(), and mk_pseudojet().

Referenced by calcMT().

206  {
207  // Make PseudoJets
208  std::vector<FJNS::PseudoJet> constituents;
209  for (const MOM* p : moms) constituents.push_back(HEPUtils::mk_pseudojet(*p));
210  // Run clustering
211  std::vector<FJNS::PseudoJet> jets = HEPUtils::get_jets(constituents, R, ptmin, alg);
212  // Make newly-allocated Jets
213  std::vector<std::shared_ptr<HEPUtils::Jet>> rtn;
214  for (const FJNS::PseudoJet& j : jets) rtn.push_back(std::make_shared<HEPUtils::Jet>(HEPUtils::mk_p4(j)));
215  return rtn;
216  }
FJNS::PseudoJet mk_pseudojet(const Vec4T &p)
Definition: Py8Utils.hpp:35
HEPUtils::P4 mk_p4(const Vec4T &p)
Definition: Py8Utils.hpp:41
std::vector< std::shared_ptr< HEPUtils::Jet > > get_jets(const std::vector< MOM *> &moms, double R, double ptmin=0 *GeV, FJNS::JetAlgorithm alg=FJNS::antikt_algorithm)
Run jet clustering from any P4-compatible momentum type.
Definition: Utils.hpp:205
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_LHC_LogLike_per_analysis()

void Gambit::ColliderBit::get_LHC_LogLike_per_analysis ( map_str_dbl result)

Extract the combined log likelihood for each analysis.

Definition at line 597 of file LHC_likelihoods.cpp.

References Gambit::ColliderBit::AnalysisLogLikes::combination_loglike, Gambit::LogTags::debug, Gambit::EOM, and Gambit::logger().

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  }
void get_LHC_LogLike_per_analysis(map_str_dbl &result)
Extract the combined log likelihood for each analysis.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
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
Here is the call graph for this function:

◆ get_LHC_LogLike_per_SR()

void Gambit::ColliderBit::get_LHC_LogLike_per_SR ( map_str_dbl result)

Extract the log likelihood for each SR.

Definition at line 618 of file LHC_likelihoods.cpp.

References Gambit::ColliderBit::AnalysisLogLikes::combination_loglike, Gambit::LogTags::debug, Gambit::EOM, Gambit::logger(), Gambit::ColliderBit::AnalysisLogLikes::sr_indices, and Gambit::ColliderBit::AnalysisLogLikes::sr_loglikes.

Referenced by get_LHC_LogLike_SR_indices(), and get_LHC_LogLike_SR_labels().

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  }
void get_LHC_LogLike_per_SR(map_str_dbl &result)
Extract the log likelihood for each SR.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_LHC_LogLike_SR_indices()

void Gambit::ColliderBit::get_LHC_LogLike_SR_indices ( map_str_dbl result)

Extract the indices for the SRs used in the analysis loglikes.

Definition at line 668 of file LHC_likelihoods.cpp.

References Gambit::ColliderBit::AnalysisLogLikes::combination_sr_index, Gambit::LogTags::debug, double, Gambit::EOM, get_LHC_LogLike_per_SR(), and Gambit::logger().

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  }
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
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
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
Here is the call graph for this function:

◆ get_LHC_LogLike_SR_labels()

void Gambit::ColliderBit::get_LHC_LogLike_SR_labels ( map_str_str result)

Extract the labels for the SRs used in the analysis loglikes.

Definition at line 653 of file LHC_likelihoods.cpp.

References Gambit::ColliderBit::AnalysisLogLikes::combination_sr_label, and get_LHC_LogLike_per_SR().

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  }
void get_LHC_LogLike_per_SR(map_str_dbl &result)
Extract the log likelihood for each SR.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
Here is the call graph for this function:

◆ get_sigma_ee_chi00()

void Gambit::ColliderBit::get_sigma_ee_chi00 ( triplet< double > &  result,
const double  sqrts,
const int  chi_first,
const int  chi_second,
const double  tol,
const bool  pt_error,
const Spectrum spec,
const double  gammaZ 
)

Retrieve the production cross-section at an e+e- collider for neutralino pairs.

Definition at line 151 of file lep_mssm_xsecs.cpp.

References alpha, Gambit::triplet< TYPE >::central, Gambit::Par::dimensionless, Gambit::Spectrum::get_HE(), LOCAL_INFO, Gambit::triplet< TYPE >::lower, Gambit::slhahelp::mass_es_from_gauge_es(), Gambit::pi, Gambit::Par::Pole_Mass, Gambit::Par::Pole_Mass_1srd_high, Gambit::Par::Pole_Mass_1srd_low, Gambit::Par::Pole_Mixing, Gambit::Spectrum::safeget(), Gambit::SubSpectrum::safeget(), SLHA2BFM_NN(), Gambit::triplet< TYPE >::upper, and xsec_neuineuj().

Referenced by LEP188_SLHA1_convention_xsec_chi00_11(), LEP188_SLHA1_convention_xsec_chi00_12(), LEP188_SLHA1_convention_xsec_chi00_13(), LEP188_SLHA1_convention_xsec_chi00_14(), LEP188_SLHA1_convention_xsec_chi00_22(), LEP188_SLHA1_convention_xsec_chi00_23(), LEP188_SLHA1_convention_xsec_chi00_24(), LEP188_SLHA1_convention_xsec_chi00_33(), LEP188_SLHA1_convention_xsec_chi00_34(), LEP188_SLHA1_convention_xsec_chi00_44(), LEP205_SLHA1_convention_xsec_chi00_11(), LEP205_SLHA1_convention_xsec_chi00_12(), LEP205_SLHA1_convention_xsec_chi00_13(), LEP205_SLHA1_convention_xsec_chi00_14(), LEP205_SLHA1_convention_xsec_chi00_22(), LEP205_SLHA1_convention_xsec_chi00_23(), LEP205_SLHA1_convention_xsec_chi00_24(), LEP205_SLHA1_convention_xsec_chi00_33(), LEP205_SLHA1_convention_xsec_chi00_34(), LEP205_SLHA1_convention_xsec_chi00_44(), LEP207_SLHA1_convention_xsec_chi00_11(), LEP208_SLHA1_convention_xsec_chi00_11(), LEP208_SLHA1_convention_xsec_chi00_12(), LEP208_SLHA1_convention_xsec_chi00_13(), LEP208_SLHA1_convention_xsec_chi00_14(), LEP208_SLHA1_convention_xsec_chi00_22(), LEP208_SLHA1_convention_xsec_chi00_23(), LEP208_SLHA1_convention_xsec_chi00_24(), LEP208_SLHA1_convention_xsec_chi00_33(), LEP208_SLHA1_convention_xsec_chi00_34(), and LEP208_SLHA1_convention_xsec_chi00_44().

153  {
154  // Subspectrum
155  const SubSpectrum& mssm = spec.get_HE();
156 
157  // PDG codes
158  const int id1 = 1000021 + chi_first + (chi_first > 2 ? 1 + (chi_first -3)*9 : 0);
159  const int id2 = 1000021 + chi_second + (chi_second > 2 ? 1 + (chi_second-3)*9 : 0);
160 
161  // SM parameters
162  const double mZ = spec.safeget(Par::Pole_Mass,23,0);
163  const double g2 = mssm.safeget(Par::dimensionless,"g2");
164  const double sinW2 = mssm.safeget(Par::dimensionless,"sinW2");
165  const double alpha = 0.25*sinW2*g2*g2/pi;
166 
167  // MSSM parameters
168  const double tanb = mssm.safeget(Par::dimensionless,"tanbeta");
169  // Get the mass eigenstates best corresponding to ~eL and ~eR.
170  const str mass_esL = slhahelp::mass_es_from_gauge_es("~e_L", mssm, tol, LOCAL_INFO, pt_error);
171  const str mass_esR = slhahelp::mass_es_from_gauge_es("~e_R", mssm, tol, LOCAL_INFO, pt_error);
172  // Get the slepton masses
173  const double mS[2] = {spec.safeget(Par::Pole_Mass,mass_esL), spec.safeget(Par::Pole_Mass,mass_esR)};
174  // Get the neutralino masses
175  const double m1 = spec.safeget(Par::Pole_Mass,id1,0);
176  const double m2 = spec.safeget(Par::Pole_Mass,id2,0);
177  std::pair<double,double> m1_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id1, 0),
178  mssm.safeget(Par::Pole_Mass_1srd_low, id1, 0));
179  std::pair<double,double> m2_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id2, 0),
180  mssm.safeget(Par::Pole_Mass_1srd_low, id2, 0));
181 
182  // Just return zero if the final state is kinematically inaccessible
183  // *even* if both masses are 2simga lower than their central values
184  if (std::abs(m1)*(1.0-2.0*m1_uncerts.second) + std::abs(m2)*(1.0-2.0*m2_uncerts.second) > sqrts)
185  {
186  result.central = 0.0;
187  result.upper = 0.0;
188  result.lower = 0.0;
189  return;
190  }
191 
192  // Get the 4x4 neutralino mixing matrix
193  MixMatrix neutmix(4,std::vector<double>(4));
194  for (int i=0; i<4; i++) for (int j=0; j<4; j++) neutmix[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi0",i+1,j+1);
195 
196  // Convert neutralino mixing matrix to BFM convention
197  SLHA2BFM_NN(neutmix, tanb, sinW2);
198 
199  // Calculate the cross-section
200  result.central = xsec_neuineuj(id1, id2, sqrts, m1, m2, neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2);
201 
202  // Calculate the uncertainty on the cross-section due to final state masses varying by +/- 1 sigma
203  std::vector<double> xsecs;
204  xsecs.push_back(result.central);
205  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.+m2_uncerts.first),
206  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
207  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.-m2_uncerts.second),
208  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
209  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.+m2_uncerts.first),
210  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
211  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.-m2_uncerts.second),
212  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
213  result.upper = *std::max_element(xsecs.begin(), xsecs.end());
214  result.lower = *std::min_element(xsecs.begin(), xsecs.end());
215 
216  }
void SLHA2BFM_NN(MixMatrix &NN, double tanb, double sin2thetaW)
Conversion between SLHA and BFM conventions.
#define LOCAL_INFO
Definition: local_info.hpp:34
START_MODEL alpha
const double pi
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
double xsec_neuineuj(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix N, const double mS[2], double tanb, double alpha, double mZ, double gZ, double sin2thetaW)
Cross section [pb] for Masses mi and mj for the neutralinos are signed.
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
std::vector< std::vector< double > > MixMatrix
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_sigma_ee_chipm()

void Gambit::ColliderBit::get_sigma_ee_chipm ( triplet< double > &  result,
const double  sqrts,
const int  chi_plus,
const int  chi_minus,
const double  tol,
const bool  pt_error,
const Spectrum spec,
const double  gammaZ 
)

Retrieve the production cross-section at an e+e- collider for chargino pairs.

Definition at line 219 of file lep_mssm_xsecs.cpp.

References alpha, Gambit::triplet< TYPE >::central, Gambit::Par::dimensionless, Gambit::Spectrum::get_HE(), LOCAL_INFO, Gambit::triplet< TYPE >::lower, Gambit::slhahelp::mass_es_from_gauge_es(), Gambit::pi, Gambit::Par::Pole_Mass, Gambit::Par::Pole_Mass_1srd_high, Gambit::Par::Pole_Mass_1srd_low, Gambit::Par::Pole_Mixing, Gambit::Spectrum::safeget(), Gambit::SubSpectrum::safeget(), SLHA2BFM_VV(), Gambit::triplet< TYPE >::upper, and xsec_chaichaj().

Referenced by LEP188_SLHA1_convention_xsec_chipm_11(), LEP188_SLHA1_convention_xsec_chipm_12(), LEP188_SLHA1_convention_xsec_chipm_22(), LEP205_SLHA1_convention_xsec_chipm_11(), LEP205_SLHA1_convention_xsec_chipm_12(), LEP205_SLHA1_convention_xsec_chipm_22(), LEP208_SLHA1_convention_xsec_chipm_11(), LEP208_SLHA1_convention_xsec_chipm_12(), and LEP208_SLHA1_convention_xsec_chipm_22().

221  {
222  // Subspectrum
223  const SubSpectrum& mssm = spec.get_HE();
224 
225  // PDG codes
226  const int id1 = 1000023 + chi_plus + (chi_plus - 1)*12;
227  const int id2 = -(1000023 + chi_minus + (chi_minus - 1)*12);
228 
229  // SM parameters
230  const double mZ = spec.safeget(Par::Pole_Mass,23,0);
231  const double g2 = mssm.safeget(Par::dimensionless,"g2");
232  const double sinW2 = mssm.safeget(Par::dimensionless,"sinW2");
233  const double alpha = 0.25*sinW2*g2*g2/pi;
234 
235  // MSSM parameters
236  // Get the mass eigenstate best corresponding to ~nu_e_L.
237  const str mass_snue = slhahelp::mass_es_from_gauge_es("~nu_e_L", mssm, tol, LOCAL_INFO, pt_error);
238  // Get the electron sneutrino mass
239  const double msn = spec.safeget(Par::Pole_Mass,mass_snue);
240  // Get the chargino masses
241  const double m1 = spec.safeget(Par::Pole_Mass,id1,0);
242  const double m2 = spec.safeget(Par::Pole_Mass,id2,0);
243  std::pair<double,double> m1_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id1, 0),
244  mssm.safeget(Par::Pole_Mass_1srd_low, id1, 0));
245  std::pair<double,double> m2_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id2, 0),
246  mssm.safeget(Par::Pole_Mass_1srd_low, id2, 0));
247 
248  // Just return zero if the final state is kinematically inaccessible
249  // *even* if both masses are 2simga lower than their central values
250  if (std::abs(m1)*(1.0-2.0*m1_uncerts.second) + std::abs(m2)*(1.0-2.0*m2_uncerts.second) > sqrts)
251  {
252  result.central = 0.0;
253  result.upper = 0.0;
254  result.lower = 0.0;
255  return;
256  }
257 
258  // Get the 2x2 chargino mixing matrices
259  MixMatrix charginomixV(2,std::vector<double>(2));
260  MixMatrix charginomixU(2,std::vector<double>(2));
261  for (int i=0; i<2; i++) for (int j=0; j<2; j++)
262  {
263  charginomixV[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi+",i+1,j+1);
264  charginomixU[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi-",i+1,j+1);
265  }
266 
267  // Convert chargino mixing matrices to BFM convention
268  SLHA2BFM_VV(charginomixV);
269  SLHA2BFM_VV(charginomixU);
270 
271  // Calculate the cross-section
272  result.central = xsec_chaichaj(id1, id2, sqrts, m1, m2, charginomixV, charginomixU,
273  msn, alpha, mZ, gammaZ, sinW2);
274 
275  // Calculate the uncertainty on the cross-section due to final state masses varying by +/- 1 sigma
276  std::vector<double> xsecs;
277  xsecs.push_back(result.central);
278  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.+m2_uncerts.first), charginomixV, charginomixU,
279  msn, alpha, mZ, gammaZ, sinW2));
280  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.-m2_uncerts.second), charginomixV, charginomixU,
281  msn, alpha, mZ, gammaZ, sinW2));
282  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.+m2_uncerts.first), charginomixV, charginomixU,
283  msn, alpha, mZ, gammaZ, sinW2));
284  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.-m2_uncerts.second), charginomixV, charginomixU,
285  msn, alpha, mZ, gammaZ, sinW2));
286  result.upper = *std::max_element(xsecs.begin(), xsecs.end());
287  result.lower = *std::min_element(xsecs.begin(), xsecs.end());
288 
289  }
double xsec_chaichaj(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix V, MixMatrix U, double msn, double alpha, double mZ, double gZ, double sin2thetaW)
Cross section [pb] for Masses mi and mj for the charginos are signed.
#define LOCAL_INFO
Definition: local_info.hpp:34
START_MODEL alpha
const double pi
void SLHA2BFM_VV(MixMatrix &VV)
Converts the chargino mixing matrix V in SLHA conventions to BFM conventions.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
std::vector< std::vector< double > > MixMatrix
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_sigma_ee_ll()

void Gambit::ColliderBit::get_sigma_ee_ll ( triplet< double > &  result,
const double  sqrts,
const int  generation,
const int  l_chirality,
const int  lbar_chirality,
const double  gtol,
const double  ftol,
const bool  gpt_error,
const bool  fpt_error,
const Spectrum spec,
const double  gammaZ,
const bool  l_are_gauge_es 
)

High-level cross section routines.

Retrieve the production cross-section at an e+e- collider for slepton pairs.

Retrieve the production cross-section at an e+e- collider for slepton pairs. If l_are_gauge_es = T, then l(bar)_chirality = 1 => (anti-)left-type slepton = 2 => (anti-)right-type slepton If l_are_gauge_es = F, then l(bar)_chirality = 1 => (anti-)slepton is lightest family state = 2 => (anti-)slepton is heaviest family state

If l_are_gauge_es = T, then l(bar)_chirality = 1 => (anti-)left-type slepton = 2 => (anti-)right-type slepton If l_are_gauge_es = F, then l(bar)_chirality = 1 => (anti-)slepton is lightest family state = 2 => (anti-)slepton is heaviest family state

Definition at line 49 of file lep_mssm_xsecs.cpp.

References alpha, Gambit::triplet< TYPE >::central, Gambit::Par::dimensionless, Gambit::slhahelp::family_state_mix_matrix(), Gambit::Spectrum::get_HE(), LOCAL_INFO, Gambit::triplet< TYPE >::lower, Gambit::slhahelp::mass_es_from_gauge_es(), Gambit::pi, Gambit::Par::Pole_Mass, Gambit::Par::Pole_Mass_1srd_high, Gambit::Par::Pole_Mass_1srd_low, Gambit::Par::Pole_Mixing, Gambit::Spectrum::safeget(), Gambit::SubSpectrum::safeget(), SLHA2BFM_NN(), Gambit::triplet< TYPE >::upper, and xsec_sleislej().

Referenced by LEP188_SLHA1_convention_xsec_se1se1bar(), LEP188_SLHA1_convention_xsec_se1se2bar(), LEP188_SLHA1_convention_xsec_se2se2bar(), LEP188_SLHA1_convention_xsec_selselbar(), LEP188_SLHA1_convention_xsec_selserbar(), LEP188_SLHA1_convention_xsec_serserbar(), LEP188_SLHA1_convention_xsec_smu1smu1bar(), LEP188_SLHA1_convention_xsec_smu1smu2bar(), LEP188_SLHA1_convention_xsec_smu2smu2bar(), LEP188_SLHA1_convention_xsec_smulsmulbar(), LEP188_SLHA1_convention_xsec_smulsmurbar(), LEP188_SLHA1_convention_xsec_smursmurbar(), LEP188_SLHA1_convention_xsec_stau1stau1bar(), LEP188_SLHA1_convention_xsec_stau1stau2bar(), LEP188_SLHA1_convention_xsec_stau2stau2bar(), LEP188_SLHA1_convention_xsec_staulstaulbar(), LEP188_SLHA1_convention_xsec_staulstaurbar(), LEP188_SLHA1_convention_xsec_staurstaurbar(), LEP205_SLHA1_convention_xsec_se1se1bar(), LEP205_SLHA1_convention_xsec_se1se2bar(), LEP205_SLHA1_convention_xsec_se2se2bar(), LEP205_SLHA1_convention_xsec_selselbar(), LEP205_SLHA1_convention_xsec_selserbar(), LEP205_SLHA1_convention_xsec_serserbar(), LEP205_SLHA1_convention_xsec_smu1smu1bar(), LEP205_SLHA1_convention_xsec_smu1smu2bar(), LEP205_SLHA1_convention_xsec_smu2smu2bar(), LEP205_SLHA1_convention_xsec_smulsmulbar(), LEP205_SLHA1_convention_xsec_smulsmurbar(), LEP205_SLHA1_convention_xsec_smursmurbar(), LEP205_SLHA1_convention_xsec_stau1stau1bar(), LEP205_SLHA1_convention_xsec_stau1stau2bar(), LEP205_SLHA1_convention_xsec_stau2stau2bar(), LEP205_SLHA1_convention_xsec_staulstaulbar(), LEP205_SLHA1_convention_xsec_staulstaurbar(), LEP205_SLHA1_convention_xsec_staurstaurbar(), LEP208_SLHA1_convention_xsec_se1se1bar(), LEP208_SLHA1_convention_xsec_se1se2bar(), LEP208_SLHA1_convention_xsec_se2se2bar(), LEP208_SLHA1_convention_xsec_selselbar(), LEP208_SLHA1_convention_xsec_selserbar(), LEP208_SLHA1_convention_xsec_serserbar(), LEP208_SLHA1_convention_xsec_smu1smu1bar(), LEP208_SLHA1_convention_xsec_smu1smu2bar(), LEP208_SLHA1_convention_xsec_smu2smu2bar(), LEP208_SLHA1_convention_xsec_smulsmulbar(), LEP208_SLHA1_convention_xsec_smulsmurbar(), LEP208_SLHA1_convention_xsec_smursmurbar(), LEP208_SLHA1_convention_xsec_stau1stau1bar(), LEP208_SLHA1_convention_xsec_stau1stau2bar(), LEP208_SLHA1_convention_xsec_stau2stau2bar(), LEP208_SLHA1_convention_xsec_staulstaulbar(), LEP208_SLHA1_convention_xsec_staulstaurbar(), and LEP208_SLHA1_convention_xsec_staurstaurbar().

52  {
53  static const str genmap[3][2] =
54  {
55  {"~e_L", "~e_R" },
56  {"~mu_L", "~mu_R" },
57  {"~tau_L", "~tau_R"}
58  };
59 
60  // Subspectrum
61  const SubSpectrum& mssm = spec.get_HE();
62  // PDG codes
63  const int id1 = 1000000*l_chirality + 11 +2*(generation-1);
64  const int id2 = -(1000000*lbar_chirality + 11 +2*(generation-1));
65 
66  // SM parameters
67  const double mZ = spec.safeget(Par::Pole_Mass,23,0);
68  const double g2 = mssm.safeget(Par::dimensionless,"g2");
69  const double sinW2 = mssm.safeget(Par::dimensionless,"sinW2");
70  const double alpha = 0.25*sinW2*g2*g2/pi;
71  // MSSM parameters
72  const double tanb = mssm.safeget(Par::dimensionless,"tanbeta");
73 
74  // Get the mass eigenstate strings and 2x2 slepton generation mass mixing matrix
75  str mass_es1, mass_es2;
76  MixMatrix sleptonmix(2,std::vector<double>(2));
77 
78  if (l_are_gauge_es)
79  {
80  // Requested final states are gauge eigenstates. Pass diagonal mixing matrix to low-level routine.
81  sleptonmix[0][0] = 1.0;
82  sleptonmix[0][1] = 0.0;
83  sleptonmix[1][0] = 0.0;
84  sleptonmix[1][1] = 1.0;
85  mass_es1 = slhahelp::mass_es_from_gauge_es(genmap[generation-1][l_chirality-1], mssm, gtol, LOCAL_INFO, gpt_error);
86  mass_es2 = slhahelp::mass_es_from_gauge_es(genmap[generation-1][lbar_chirality-1], mssm, gtol, LOCAL_INFO, gpt_error);
87  }
88  else
89  {
90  // Requested final states are family mass eigenstates. Pass 2x2 family mass mixing matrix to low-level routine.
91  str m_light, m_heavy;
92  std::vector<double> slepton4vec = slhahelp::family_state_mix_matrix("~e-", generation, m_light, m_heavy, mssm, ftol, LOCAL_INFO, fpt_error);
93  mass_es1 = (l_chirality == 1) ? m_light : m_heavy;
94  mass_es2 = (lbar_chirality == 1) ? m_light : m_heavy;
95  sleptonmix[0][0] = slepton4vec[0];
96  sleptonmix[0][1] = slepton4vec[1];
97  sleptonmix[1][0] = slepton4vec[2];
98  sleptonmix[1][1] = slepton4vec[3];
99  }
100  const double m1 = spec.safeget(Par::Pole_Mass,mass_es1);
101  // FIXME when spectrum object has separate pole mass getters for antiparticles
102  //const double m2 = spec.safeget(Par::Pole_Mass,Models::ParticleDB().get_antiparticle(mass_es2));
103  // until then
104  const double m2 = spec.safeget(Par::Pole_Mass,mass_es2);
105  std::pair<double,double> m1_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, mass_es1),
106  mssm.safeget(Par::Pole_Mass_1srd_low, mass_es1));
107  std::pair<double,double> m2_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, mass_es2),
108  mssm.safeget(Par::Pole_Mass_1srd_low, mass_es2));
109 
110  // If the final state is kinematically inaccessible *even* if both masses
111  // are 2simga lower than their central values, then return zero.
112  if (m1*(1.0-2.0*m1_uncerts.second) + m2*(1.0-2.0*m2_uncerts.second) > sqrts)
113  {
114  result.central = 0.0;
115  result.upper = 0.0;
116  result.lower = 0.0;
117  return;
118  }
119 
120  // Get the neutralino masses
121  const double neutmass[4] = { spec.safeget(Par::Pole_Mass,1000022,0), spec.safeget(Par::Pole_Mass,1000023,0),
122  spec.safeget(Par::Pole_Mass,1000025,0), spec.safeget(Par::Pole_Mass,1000035,0) };
123  // Get the 4x4 neutralino mixing matrix
124  MixMatrix neutmix(4,std::vector<double>(4));
125  for (int i=0; i<4; i++) for (int j=0; j<4; j++) neutmix[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi0",i+1,j+1);
126 
127  // Convert neutralino mixing matrix to BFM convention
128  SLHA2BFM_NN(neutmix, tanb, sinW2);
129 
130  // Calculate the cross-section
131  result.central = xsec_sleislej(id1, id2, sqrts, m1, m2, sleptonmix, neutmix, neutmass, alpha, mZ, gammaZ, sinW2);
132 
133  // Calculate the uncertainty on the cross-section due to final state masses varying by +/- 1 sigma
134  std::vector<double> xsecs;
135  xsecs.push_back(result.central);
136  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.+m2_uncerts.first), sleptonmix, neutmix,
137  neutmass, alpha, mZ, gammaZ, sinW2, false));
138  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.+m2_uncerts.first), sleptonmix, neutmix,
139  neutmass, alpha, mZ, gammaZ, sinW2, false));
140  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.-m2_uncerts.second), sleptonmix, neutmix,
141  neutmass, alpha, mZ, gammaZ, sinW2, false));
142  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.-m2_uncerts.second), sleptonmix, neutmix,
143  neutmass, alpha, mZ, gammaZ, sinW2, false));
144  result.upper = *std::max_element(xsecs.begin(), xsecs.end());
145  result.lower = *std::min_element(xsecs.begin(), xsecs.end());
146 
147  }
std::vector< double > family_state_mix_matrix(str type, int generation, str &mass_es1, str &mass_es2, const SubSpectrum &mssm, double tol, str context, bool pterror)
Get the family mixing matrix and corresponding mass eigenstates, then check for interfamily mixing...
void SLHA2BFM_NN(MixMatrix &NN, double tanb, double sin2thetaW)
Conversion between SLHA and BFM conventions.
#define LOCAL_INFO
Definition: local_info.hpp:34
START_MODEL alpha
double xsec_sleislej(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix F, MixMatrix N, const double mN[4], double alpha, double mZ, double gZ, double sin2thetaW, bool warn_on_CP_violating_masses=true)
Low-level cross section routines.
const double pi
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
std::vector< std::vector< double > > MixMatrix
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GET_SPECIFIC_PYTHIA() [1/2]

Gambit::ColliderBit::GET_SPECIFIC_PYTHIA ( getPythia_EM  ,
Pythia_EM_default  ,
EM_spectrum  ,
_EM  ,
NOT_SUSY   
)

◆ GET_SPECIFIC_PYTHIA() [2/2]

Gambit::ColliderBit::GET_SPECIFIC_PYTHIA ( getPythia  ,
Pythia_default  ,
MSSM_spectrum  ,
IS_SUSY   
)

◆ getAnalysisContainer()

void Gambit::ColliderBit::getAnalysisContainer ( AnalysisContainer result,
const str detname,
const MCLoopInfo RunMC,
const xsec CrossSection,
int  iteration 
)

Retrieve an analysis container for a specific detector.

Definition at line 49 of file getAnalysisContainer.cpp.

References Gambit::ColliderBit::MCLoopInfo::analyses, BASE_INIT, Gambit::ColliderBit::AnalysisContainer::collect_and_add_signal(), COLLIDER_FINALIZE, Gambit::ColliderBit::MCLoopInfo::current_analyses_exist_for(), Gambit::ColliderBit::MCLoopInfo::current_analyses_for(), Gambit::ColliderBit::MCLoopInfo::current_collider(), Gambit::ColliderBit::AnalysisContainer::has_analyses(), Gambit::ColliderBit::AnalysisContainer::init(), LOCAL_INFO, Gambit::piped_errors, Gambit::ColliderBit::AnalysisContainer::register_thread(), Gambit::Piped_exceptions::request(), Gambit::ColliderBit::AnalysisContainer::reset(), Gambit::ColliderBit::AnalysisContainer::scale(), Gambit::ColliderBit::AnalysisContainer::set_current_collider(), START_SUBPROCESS, and Gambit::ColliderBit::xsec::xsec_per_event().

54  {
55  if (RunMC.analyses.empty() or iteration == BASE_INIT) return;
56  if (not RunMC.current_analyses_exist_for(detname)) return;
57 
58  if (iteration == START_SUBPROCESS)
59  {
60  // Register analysis container
61  result.register_thread(detname+"AnalysisContainer");
62 
63  // Set current collider
64  result.set_current_collider(RunMC.current_collider());
65 
66  // Initialize analysis container or reset all the contained analyses
67  if (!result.has_analyses())
68  {
69  try
70  {
71  result.init(RunMC.current_analyses_for(detname));
72  }
73  catch (std::runtime_error& e)
74  {
75  piped_errors.request(LOCAL_INFO, e.what());
76  }
77  }
78  else result.reset();
79  }
80 
81  if (iteration == COLLIDER_FINALIZE)
82  {
83  result.collect_and_add_signal();
84  result.scale(CrossSection.xsec_per_event());
85  }
86 
87  }
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:527
#define LOCAL_INFO
Definition: local_info.hpp:34
Here is the call graph for this function:

◆ getBuckFast()

template<typename EventT >
BaseDetector<EventT>* Gambit::ColliderBit::getBuckFast ( const str detname,
const MCLoopInfo RunMC,
int  iteration,
const Options runOptions 
)

Get a BuckFast detector simulation.

Definition at line 58 of file getBuckFast.hpp.

References COLLIDER_INIT, Gambit::ColliderBit::MCLoopInfo::current_collider(), Gambit::Options::getValue(), Gambit::Options::getValueOrDef(), Gambit::Options::hasKey(), LOCAL_INFO, Gambit::ColliderBit::CMS::smearElectronEnergy(), Gambit::ColliderBit::ATLAS::smearElectronEnergy(), Gambit::ColliderBit::CMS::smearJets(), Gambit::ColliderBit::ATLAS::smearJets(), Gambit::ColliderBit::CMS::smearMuonMomentum(), Gambit::ColliderBit::ATLAS::smearMuonMomentum(), Gambit::ColliderBit::CMS::smearTaus(), Gambit::ColliderBit::ATLAS::smearTaus(), and START_SUBPROCESS.

62  {
63  static bool partonOnly;
64  static double antiktR;
65 
66  // Where the real action is
67  static std::unique_ptr<BuckFast<EventT>[]> bucky(new BuckFast<EventT>[omp_get_max_threads()]);
68  int mine = omp_get_thread_num();
69 
70  if (iteration == COLLIDER_INIT)
71  {
72  // The default values for partonOnly and antiktR
73  bool default_partonOnly = false;
74  double default_antiktR = 0.4;
75  // Retrieve any yaml options specifying partonOnly and antikt R for this collider
76  if (runOptions.hasKey(RunMC.current_collider()))
77  {
78  Options colOptions(runOptions.getValue<YAML::Node>(RunMC.current_collider()));
79  partonOnly = colOptions.getValueOrDef<bool>(default_partonOnly, "partonOnly");
80  antiktR = colOptions.getValueOrDef<double>(default_antiktR, "antiktR");
81  }
82  else
83  {
84  partonOnly = default_partonOnly;
85  antiktR = default_antiktR;
86  }
87  }
88 
89  if (iteration == START_SUBPROCESS)
90  {
91  // Each thread gets its own copy of the detector sim, so it is initialised *after* COLLIDER_INIT, within omp parallel.
92  bucky[mine].init(partonOnly, antiktR);
93  // Assign detector functions
94  if (detname == "ATLAS")
95  {
96  bucky[mine].smearElectronEnergy = &ATLAS::smearElectronEnergy;
97  bucky[mine].smearMuonMomentum = &ATLAS::smearMuonMomentum ;
98  bucky[mine].smearTaus = &ATLAS::smearTaus;
99  bucky[mine].smearJets = &ATLAS::smearJets;
100  }
101  else if (detname == "CMS")
102  {
103  bucky[mine].smearElectronEnergy = &CMS::smearElectronEnergy;
104  bucky[mine].smearMuonMomentum = &CMS::smearMuonMomentum;
105  bucky[mine].smearTaus = &CMS::smearTaus;
106  bucky[mine].smearJets = &CMS::smearJets;
107  }
108  else if (detname == "Identity") { /* relax */ }
109  else
110  {
111  ColliderBit_error().raise(LOCAL_INFO, "Unrecognised detector name.");
112  }
113  }
114 
115  // Paper-bag it
116  return &bucky[mine];
117 
118  }
#define LOCAL_INFO
Definition: local_info.hpp:34
void smearElectronEnergy(std::vector< HEPUtils::Particle *> &electrons)
Randomly smear the supplied electrons&#39; momenta by parameterised resolutions.
void smearMuonMomentum(std::vector< HEPUtils::Particle *> &muons)
Randomly smear the supplied muons&#39; momenta by parameterised resolutions.
void smearJets(std::vector< HEPUtils::Jet *> &jets)
Randomly smear the supplied jets&#39; momenta by parameterised resolutions.
void smearTaus(std::vector< HEPUtils::Particle *> &taus)
Randomly smear the supplied taus&#39; momenta by parameterised resolutions.
Here is the call graph for this function:

◆ getColliderPythia()

template<typename PythiaT , typename EventT >
void Gambit::ColliderBit::getColliderPythia ( ColliderPythia< PythiaT, EventT > &  result,
const MCLoopInfo RunMC,
const Spectrum spectrum,
const DecayTable decay_rates,
const str  model_suffix,
const int  iteration,
void(*)()  wrapup,
const Options runOptions,
bool  is_SUSY 
)

Retrieve a Pythia hard-scattering Monte Carlo simulation.

Definition at line 51 of file getColliderPythia.hpp.

References Gambit::ColliderBit::ColliderPythia< PythiaT, EventT >::banner(), BASE_FINALIZE, BASE_INIT, Gambit::ColliderBit::ColliderPythia< PythiaT, EventT >::clear(), COLLIDER_INIT, Gambit::ColliderBit::MCLoopInfo::current_collider(), debug_prefix(), Gambit::Random::draw(), Gambit::EOM, Gambit::DecayTable::getSLHAea(), Gambit::Spectrum::getSLHAea(), Gambit::Options::getValue(), Gambit::Options::hasKey(), Gambit::ColliderBit::ColliderPythia< PythiaT, EventT >::init(), int, Gambit::invalid_point(), Gambit::logger(), Gambit::piped_invalid_point, Gambit::invalid_point_exception::raise(), Gambit::Piped_invalid_point::request(), and START_SUBPROCESS.

60  {
61  static bool first = true;
62  static std::vector<str> filenames;
63  static str pythia_doc_path;
64  static std::vector<str> pythiaCommonOptions;
65  static SLHAstruct slha;
66  static SLHAstruct slha_spectrum;
67  static double xsec_veto_fb;
68  static unsigned int fileCounter = 0;
69 
70  if (iteration == BASE_INIT)
71  {
72  // Setup the Pythia documentation path and print the banner once
73  if (first)
74  {
75  const str be = "Pythia" + model_suffix;
76  const str ver = Backends::backendInfo().default_version(be);
77  pythia_doc_path = Backends::backendInfo().path_dir(be, ver) + "/../share/Pythia8/xmldoc/";
78  result.banner(pythia_doc_path);
79  if (runOptions.hasKey("SLHA_filenames"))
80  {
81  filenames = runOptions.getValue<std::vector<str> >("SLHA_filenames");
82  }
83  first = false;
84  }
85 
86  if (filenames.empty())
87  {
88  // SLHAea object constructed from dependencies on the spectrum and decays.
89  slha.clear();
90  slha_spectrum.clear();
91  slha = decay_rates.getSLHAea(2);
92  // SLHAea in SLHA2 format, please.
93  slha_spectrum = spectrum.getSLHAea(2);
94  slha.insert(slha.begin(), slha_spectrum.begin(), slha_spectrum.end());
95  if (is_SUSY)
96  {
97  SLHAea::Block block("MODSEL");
98  block.push_back("BLOCK MODSEL # Model selection");
99  SLHAea::Line line;
100  line << 1 << 0 << "# Tell Pythia that this is a SUSY model.";
101  block.push_back(line);
102  slha.push_front(block);
103  }
104  }
105  else
106  {
107  if (filenames.size() <= fileCounter) invalid_point().raise("No more SLHA files. My work is done.");
108  }
109 
110  }
111 
112  else if (iteration == COLLIDER_INIT)
113  {
114  // Collect Pythia options that are common across all OMP threads
115  pythiaCommonOptions.clear();
116 
117  // By default we tell Pythia to be quiet. (Can be overridden from yaml settings)
118  pythiaCommonOptions.push_back("Print:quiet = on");
119  pythiaCommonOptions.push_back("SLHA:verbose = 0");
120 
121  // Get options from yaml file.
122  double xsec_veto_default = 0.0;
123  if (runOptions.hasKey(RunMC.current_collider()))
124  {
125  YAML::Node colNode = runOptions.getValue<YAML::Node>(RunMC.current_collider());
126  Options colOptions(colNode);
127  xsec_veto_fb = colOptions.getValueOrDef<double>(xsec_veto_default, "xsec_veto");
128 
129  if (colOptions.hasKey("pythia_settings"))
130  {
131  std::vector<str> addPythiaOptions = colNode["pythia_settings"].as<std::vector<str> >();
132  pythiaCommonOptions.insert(pythiaCommonOptions.end(), addPythiaOptions.begin(), addPythiaOptions.end());
133  }
134  }
135  else xsec_veto_fb = xsec_veto_default;
136 
137  // We need showProcesses for the xsec veto.
138  pythiaCommonOptions.push_back("Init:showProcesses = on");
139 
140  // We need "SLHA:file = slhaea" for the SLHAea interface, and the filename for the SLHA interface.
141  str slha_string = (filenames.empty() ? "slhaea" : filenames.at(fileCounter));
142  pythiaCommonOptions.push_back("SLHA:file = " + slha_string);
143  }
144 
145  else if (iteration == START_SUBPROCESS)
146  {
147  // Variables needed for the xsec veto
148  std::stringstream processLevelOutput;
149  str _junk, readline;
150  int code, nxsec;
151  double xsec, totalxsec;
152 
153  // Each thread needs an independent Pythia instance at the start
154  // of each event generation loop.
155  // Thus, the actual Pythia initialization is
156  // *after* COLLIDER_INIT, within omp parallel.
157 
158  result.clear();
159 
160  if (not filenames.empty() and omp_get_thread_num() == 0)
161  logger() << "Reading SLHA file: " << filenames.at(fileCounter) << EOM;
162 
163  // Get the Pythia options that are common across all OMP threads ('pythiaCommonOptions')
164  // and then add the thread-specific seed
165  std::vector<str> pythiaOptions = pythiaCommonOptions;
166  str seed = std::to_string(int(Random::draw() * 899990000.));
167  pythiaOptions.push_back("Random:seed = " + seed);
168 
169  #ifdef COLLIDERBIT_DEBUG
170  cout << debug_prefix() << "getPythia"+model_suffix+": My Pythia seed is: " << seed << endl;
171  #endif
172 
173  try
174  {
175  if (filenames.empty())
176  {
177  result.init(pythia_doc_path, pythiaOptions, &slha, processLevelOutput);
178  }
179  else
180  {
181  result.init(pythia_doc_path, pythiaOptions, processLevelOutput);
182  }
183  }
184  catch (typename ColliderPythia<PythiaT,EventT>::InitializationError& e)
185  {
186  // Append new seed to override the previous one
187  int newSeedBase = int(Random::draw() * 899990000.);
188  pythiaOptions.push_back("Random:seed = " + std::to_string(newSeedBase));
189  try
190  {
191  if (filenames.empty())
192  {
193  result.init(pythia_doc_path, pythiaOptions, &slha, processLevelOutput);
194  }
195  else
196  {
197  result.init(pythia_doc_path, pythiaOptions, processLevelOutput);
198  }
199  }
200  catch (typename ColliderPythia<PythiaT,EventT>::InitializationError& e)
201  {
202  #ifdef COLLIDERBIT_DEBUG
203  cout << debug_prefix() << "ColliderPythia::InitializationError caught in getColliderPythia. Will discard this point." << endl;
204  #endif
205  piped_invalid_point.request("Bad point: Pythia can't initialize");
206  wrapup();
207  return;
208  }
209  }
210 
211  // Should we apply the xsec veto and skip event generation?
212 
213  // - Get the upper limt xsec as estimated by Pythia
214  code = -1;
215  nxsec = 0;
216  totalxsec = 0.;
217  while(true)
218  {
219  std::getline(processLevelOutput, readline);
220  std::istringstream issPtr(readline);
221  issPtr.seekg(47, issPtr.beg);
222  issPtr >> code;
223  if (!issPtr.good() && nxsec > 0) break;
224  issPtr >> _junk >> xsec;
225  if (issPtr.good())
226  {
227  totalxsec += xsec;
228  nxsec++;
229  }
230  }
231 
232  #ifdef COLLIDERBIT_DEBUG
233  cout << debug_prefix() << "totalxsec [fb] = " << totalxsec * 1e12 << ", veto limit [fb] = " << xsec_veto_fb << endl;
234  #endif
235 
236  // - Check for NaN xsec
237  if (Utils::isnan(totalxsec))
238  {
239  #ifdef COLLIDERBIT_DEBUG
240  cout << debug_prefix() << "Got NaN cross-section estimate from Pythia." << endl;
241  #endif
242  piped_invalid_point.request("Got NaN cross-section estimate from Pythia.");
243  wrapup();
244  return;
245  }
246 
247  // - Wrap up loop if veto applies
248  if (totalxsec * 1e12 < xsec_veto_fb)
249  {
250  #ifdef COLLIDERBIT_DEBUG
251  cout << debug_prefix() << "Cross-section veto applies. Will now call Loop::wrapup() to skip event generation for this collider." << endl;
252  #endif
253  wrapup();
254  }
255 
256  }
257 
258  else if (iteration == BASE_FINALIZE and not filenames.empty()) fileCounter++;
259 
260  }
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
SLHAea::Coll SLHAstruct
Less confusing name for SLHAea container class.
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
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
DS5_MSPCTM DS_INTDOF int
void request(std::string message)
Request an exception.
Definition: exceptions.cpp:463
invalid_point_exception & invalid_point()
Invalid point exceptions.
str debug_prefix()
Debug prefix string giving thread number.
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 DecayTable::Entry DecayTable::Entry DecayTable::Entry decay_rates
Here is the call graph for this function:

◆ getDetector()

str Gambit::ColliderBit::getDetector ( const str name)

Return the detector to be used for a given analysis name, checking that the analysis exists.

Return the detector to be used for a given analysis name (and check that the analysis exists).

Definition at line 154 of file AnalysisContainer.cpp.

References IF_X_RTN_DETECTOR, Gambit::ColliderBit::AnalysisContainer::instances_map, LOCAL_INFO, MAP_ANALYSES, MAP_ANALYSES_WITH_ROOT, MAP_ANALYSES_WITH_ROOT_RESTFRAMES(), and Gambit::utils_error().

Referenced by operateLHCLoop().

155  {
156  #ifndef EXCLUDE_ROOT
157  #ifndef EXCLUDE_RESTFRAMES
159  #endif
161  #endif
163 
164  // If we end up here the analysis has not been found
165  utils_error().raise(LOCAL_INFO, "The analysis " + name + " is not a known ColliderBit analysis.");
166  return "";
167  }
EXPORT_SYMBOLS error & utils_error()
Utility errors.
#define LOCAL_INFO
Definition: local_info.hpp:34
#define MAP_ANALYSES_WITH_ROOT(F)
#define MAP_ANALYSES(F)
#define MAP_ANALYSES_WITH_ROOT_RESTFRAMES(F)
#define IF_X_RTN_DETECTOR(A)
For the string-based analysis checker and detector retriever getDetector.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ getDummyColliderObservable()

void Gambit::ColliderBit::getDummyColliderObservable ( double result)

Dummy observable that creates a dependency on TestModel1D.

This is used to satisfy the normal GAMBIT model requrements in a minimal way. This is useful in the case where we just want to run ColliderBit on a single point with a custom Pythia version, using Pythia's SLHA interface.

Definition at line 38 of file ColliderBit_dummy.cpp.

39  {
40  result = 0.0;
41  }

◆ getLHCEventLoopInfo()

void Gambit::ColliderBit::getLHCEventLoopInfo ( map_str_dbl result)

Store some information about the event generation.

Definition at line 256 of file ColliderBit_eventloop.cpp.

References double.

257  {
258  using namespace Pipes::getLHCEventLoopInfo;
259  result.clear();
260  result["did_event_generation"] = double(Dep::RunMC->event_generation_began);
261  result["too_many_failed_events"] = double(Dep::RunMC->exceeded_maxFailedEvents);
262  for (auto& name : Dep::RunMC->collider_names)
263  {
264  result["event_count_" + name] = Dep::RunMC->event_count.at(name);
265  }
266  }
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
void getLHCEventLoopInfo(map_str_dbl &result)
Store some information about the event generation.

◆ getMCxsec()

void Gambit::ColliderBit::getMCxsec ( xsec result)

Compute a cross-section from Monte Carlo.

Definition at line 28 of file getxsec.cpp.

References COLLIDER_FINALIZE, debug_prefix(), END_SUBPROCESS, Gambit::ColliderBit::xsec::gather_xsecs(), Gambit::ColliderBit::xsec::log_event(), Gambit::ColliderBit::xsec::reset(), Gambit::ColliderBit::xsec::set_xsec(), and START_SUBPROCESS.

29  {
30  using namespace Pipes::getMCxsec;
31 
32  // Don't bother if there are no analyses that will use this.
33  if (Dep::RunMC->analyses.empty()) return;
34 
35  // Reset the xsec objects on all threads
36  if (*Loop::iteration == START_SUBPROCESS)
37  {
38  result.reset();
39  }
40 
41  // If we are in the main event loop, count the event towards cross-section normalisation on this thread
42  if (*Loop::iteration > 0)
43  {
44  result.log_event();
45  }
46 
47  // Extract the xsecs from the MC on each thread
48  if (*Loop::iteration == END_SUBPROCESS && Dep::RunMC->event_generation_began)
49  {
50  if (not Dep::RunMC->exceeded_maxFailedEvents)
51  {
52  const double xs_fb = (*Dep::HardScatteringSim)->xsec_pb() * 1000.;
53  const double xserr_fb = (*Dep::HardScatteringSim)->xsecErr_pb() * 1000.;
54  result.set_xsec(xs_fb, xserr_fb);
55  #ifdef COLLIDERBIT_DEBUG
56  cout << debug_prefix() << "xs_fb = " << xs_fb << " +/- " << xserr_fb << endl;
57  #endif
58  }
59  }
60 
61  // Gather the xsecs from all threads into one
62  if (*Loop::iteration == COLLIDER_FINALIZE)
63  {
64  result.gather_xsecs();
65  }
66 
67  }
void getMCxsec(xsec &result)
Compute a cross-section from Monte Carlo.
Definition: getxsec.cpp:28
str debug_prefix()
Debug prefix string giving thread number.
Here is the call graph for this function:

◆ getNLLFastxsec()

void Gambit::ColliderBit::getNLLFastxsec ( xsec result)

Get a cross-section from NLL-FAST.

Definition at line 70 of file getxsec.cpp.

References COLLIDER_FINALIZE, COLLIDER_INIT, Gambit::ColliderBit::xsec::log_event(), Gambit::ColliderBit::xsec::reset(), and Gambit::ColliderBit::xsec::set_xsec().

71  {
72  using namespace Pipes::getNLLFastxsec;
73 
74  // Don't bother if there are no analyses that will use this.
75  if (Dep::RunMC->analyses.empty()) return;
76 
77  // Reset the xsec object on the main thread (other threads do not matter)
78  if (*Loop::iteration == COLLIDER_INIT)
79  {
80  result.reset();
81  }
82 
83  // If we are in the main event loop, count the event towards cross-section normalisation on this thread
84  if (*Loop::iteration > 0)
85  {
86  result.log_event();
87  }
88 
89  // Set the xsec and its error
90  if (*Loop::iteration == COLLIDER_FINALIZE)
91  {
92  double xs_fb = 0.1; // replace with xsec from NLL-Fast
93  double xserr_fb = 0.1 * xs_fb; // or whatever
94  result.set_xsec(xs_fb, xserr_fb);
95  }
96 
97  }
void getNLLFastxsec(xsec &result)
Get a cross-section from NLL-FAST.
Definition: getxsec.cpp:70
Here is the call graph for this function:

◆ getOSpairs()

std::vector< std::vector< HEPUtils::Particle * > > Gambit::ColliderBit::getOSpairs ( std::vector< HEPUtils::Particle *>  particles)

Utility function for returning a collection of oppsosite-sign particle pairs.

Definition at line 86 of file Utils.cpp.

References combine_hdf5::pid.

Referenced by none_of().

86  {
87  std::vector<std::vector<HEPUtils::Particle*>> OSpair_container;
88  for (size_t ip1=0;ip1<particles.size();ip1++) {
89  for (size_t ip2=ip1+1; ip2<particles.size(); ip2++) {
90  if (particles[ip1]->pid()*particles[ip2]->pid()<0.) {
91  std::vector<HEPUtils::Particle*> OSpair;
92  OSpair.push_back(particles[ip1]);
93  OSpair.push_back(particles[ip2]);
94  OSpair_container.push_back(OSpair);
95  }
96  }
97  }
98  return OSpair_container;
99  }
Here is the caller graph for this function:

◆ getSFOSpairs()

std::vector< std::vector< HEPUtils::Particle * > > Gambit::ColliderBit::getSFOSpairs ( std::vector< HEPUtils::Particle *>  particles)

Utility function for returning a collection of same-flavour, oppsosite-sign particle pairs.

Definition at line 69 of file Utils.cpp.

Referenced by none_of().

69  {
70  std::vector<std::vector<HEPUtils::Particle*>> SFOSpair_container;
71  for (size_t ip1=0; ip1<particles.size(); ip1++) {
72  for (size_t ip2=ip1+1; ip2<particles.size(); ip2++) {
73  if (particles[ip1]->abspid()==particles[ip2]->abspid() && particles[ip1]->pid()!=particles[ip2]->pid()) {
74  std::vector<HEPUtils::Particle*> SFOSpair;
75  SFOSpair.push_back(particles[ip1]);
76  SFOSpair.push_back(particles[ip2]);
77  SFOSpair_container.push_back(SFOSpair);
78  }
79  }
80  }
81  return SFOSpair_container;
82  }
Here is the caller graph for this function:

◆ getSSpairs()

std::vector< std::vector< HEPUtils::Particle * > > Gambit::ColliderBit::getSSpairs ( std::vector< HEPUtils::Particle *>  particles)

Utility function for returning a collection of same-sign particle pairs.

Definition at line 103 of file Utils.cpp.

References combine_hdf5::pid.

Referenced by none_of().

103  {
104  std::vector<std::vector<HEPUtils::Particle*>> SSpair_container;
105  for (size_t ip1=0;ip1<particles.size();ip1++) {
106  for (size_t ip2=ip1+1; ip2<particles.size(); ip2++) {
107  if (particles[ip1]->pid()*particles[ip2]->pid()>0.) {
108  std::vector<HEPUtils::Particle*> SSpair;
109  SSpair.push_back(particles[ip1]);
110  SSpair.push_back(particles[ip2]);
111  SSpair_container.push_back(SSpair);
112  }
113  }
114  }
115  return SSpair_container;
116  }
Here is the caller graph for this function:

◆ has_tag()

bool Gambit::ColliderBit::has_tag ( const HEPUtils::BinnedFn2D< double > &  effmap,
double  eta,
double  pt 
)
inline

Randomly get a tag result (can be anything) from a 2D |eta|-pT efficiency map.

Definition at line 173 of file Utils.hpp.

References random_bool().

Referenced by _Phi_mpi_pi(), calcMT(), calcMT_1l(), sortByPT(), sortByPT_2lep(), sortByPT_l(), sortByPT_lep(), and SortJets().

173  {
174  try {
175  return random_bool(effmap, fabs(eta), pt);
176  } catch (...) {
177  return false; // No tag if outside lookup range... be careful!
178  }
179  }
bool random_bool(const HEPUtils::BinnedFn2D< double > &effmap, double x, double y)
Return a random true/false at a success rate given by a 2D efficiency map.
Definition: Utils.hpp:143
Here is the call graph for this function:
Here is the caller graph for this function:

◆ I1()

double Gambit::ColliderBit::I1 ( double  s,
double  m1,
double  m2,
double  mk,
double  ml 
)

Integrals for t-channel neutralino diagrams m1 and m2 are masses of final state sleptons mk and ml are neutralino masses.

Definition at line 295 of file lep_mssm_xsecs.cpp.

References pow2.

Referenced by Gambit::DecayBit::stau_1_decays_smallsplit(), and xsec_sleislej().

296  {
297  double S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
298  double m1sq = pow2(m1);
299  double m2sq = pow2(m2);
300  double mksq = pow2(mk);
301  double mlsq = pow2(ml);
302 
303  double I1 = 0;
304  // Careful with degenerate masses!
305  if( fabs(mksq-mlsq) < 0.1 ){
306  I1 = (m1sq+m2sq-2.*mksq-s)*log((m1sq+m2sq-2.*mksq-s+S)/(m1sq+m2sq-2.*mksq-s-S))-4*S*((m1sq-mksq)*(m2sq-mksq)+mksq*s)/(m1sq+m2sq-2.*mksq-s-S)/(m1sq+m2sq-2.*mksq-s+S)-S;
307  }
308  else{
309  I1 = (m1sq*(m2sq-mksq)+mksq*(mksq-m2sq+s))/(mksq-mlsq)*log((m1sq+m2sq-2.*mksq-s-S)/(m1sq+m2sq-2.*mksq-s+S));
310  I1 += (m1sq*(m2sq-mlsq)+mlsq*(mlsq-m2sq+s))/(mlsq-mksq)*log((m1sq+m2sq-2.*mlsq-s-S)/(m1sq+m2sq-2.*mlsq-s+S));
311  I1 -= S;
312  }
313  return I1;
314  }
double I1(double s, double m1, double m2, double mk, double ml)
Integrals for t-channel neutralino diagrams m1 and m2 are masses of final state sleptons mk and ml ar...
#define pow2(a)
Here is the caller graph for this function:

◆ I2()

double Gambit::ColliderBit::I2 ( double  s,
double  m1,
double  m2,
double  mk,
double  ml 
)

Definition at line 315 of file lep_mssm_xsecs.cpp.

References pow2.

Referenced by Gambit::DecayBit::stau_1_decays_smallsplit(), and xsec_sleislej().

316  {
317  double S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
318  double m1sq = pow2(m1);
319  double m2sq = pow2(m2);
320  double mksq = pow2(mk);
321  double mlsq = pow2(ml);
322 
323  double I2 = 0;
324  // Careful with degenerate masses!
325  if( fabs(mksq-mlsq) < 0.1 )
326  {
327  I2 = S/(m1sq*(m2sq-mksq)+mksq*(-m2sq+mksq+s));
328  }
329  else
330  {
331  I2 = log((m1sq+m2sq-2.*mksq-(s-S))/(m1sq+m2sq-2.*mksq-(s+S)));
332  I2 += log((m1sq+m2sq-2.*mlsq-(s+S))/(m1sq+m2sq-2.*mlsq-(s-S)));
333  I2 *= 1./(mksq-mlsq);
334  }
335  return I2;
336  }
double I2(double s, double m1, double m2, double mk, double ml)
#define pow2(a)
Here is the caller graph for this function:

◆ I3()

double Gambit::ColliderBit::I3 ( double  s,
double  m1,
double  m2,
double  mk 
)

Definition at line 337 of file lep_mssm_xsecs.cpp.

References pow2.

Referenced by xsec_sleislej().

338  {
339  double S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
340  double m1sq = pow2(m1);
341  double m2sq = pow2(m2);
342  double mksq = pow2(mk);
343 
344  double I3 = 0;
345  I3 = log((m1sq+m2sq-2.*mksq-(s+S))/(m1sq+m2sq-2.*mksq-(s-S)));
346  I3 *= m1sq*(m2sq-mksq)+mksq*(-m2sq+mksq+s);
347  I3 += (m1sq+m2sq-2.*mksq-s)*S/2.;
348  return I3;
349  }
#define pow2(a)
double I3(double s, double m1, double m2, double mk)
Here is the caller graph for this function:

◆ ifilter_reject() [1/2]

void Gambit::ColliderBit::ifilter_reject ( ParticlePtrs particles,
std::function< bool(Particle *)>  rejfn,
bool  do_delete = true 
)
inline

In-place filter a supplied particle vector by rejecting those which fail a supplied cut.

Definition at line 53 of file Utils.hpp.

References iremoveerase().

Referenced by filter_reject(), ifilter_select(), and removeOverlap().

54  {
55  iremoveerase(particles, [&](Particle* p) {
56  const bool rm = rejfn(p);
57  if (rm && do_delete) delete p;
58  return rm;
59  });
60  }
void iremoveerase(CONTAINER &c, const RMFN &fn)
Convenience combination of remove_if and erase.
Definition: Utils.hpp:47
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ifilter_reject() [2/2]

void Gambit::ColliderBit::ifilter_reject ( JetPtrs jets,
std::function< bool(Jet *)>  rejfn,
bool  do_delete = true 
)
inline

In-place filter a supplied jet vector by rejecting those which fail a supplied cut.

Definition at line 92 of file Utils.hpp.

References iremoveerase().

93  {
94  iremoveerase(jets, [&](Jet* j) {
95  const bool rm = rejfn(j);
96  if (rm && do_delete) delete j;
97  return rm;
98  });
99  }
void iremoveerase(CONTAINER &c, const RMFN &fn)
Convenience combination of remove_if and erase.
Definition: Utils.hpp:47
Here is the call graph for this function:

◆ ifilter_select() [1/2]

void Gambit::ColliderBit::ifilter_select ( ParticlePtrs particles,
std::function< bool(Particle *)>  selfn,
bool  do_delete = true 
)
inline

In-place filter a supplied particle vector by keeping those which pass a supplied cut.

Definition at line 63 of file Utils.hpp.

References ifilter_reject().

64  {
65  ifilter_reject(particles, [&](Particle* p) { return !selfn(p); }, do_delete);
66  }
void ifilter_reject(JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
In-place filter a supplied jet vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:92
Here is the call graph for this function:

◆ ifilter_select() [2/2]

void Gambit::ColliderBit::ifilter_select ( JetPtrs jets,
std::function< bool(Jet *)>  selfn,
bool  do_delete = true 
)
inline

In-place filter a supplied jet vector by keeping those which pass a supplied cut.

Definition at line 102 of file Utils.hpp.

References ifilter_reject().

103  {
104  ifilter_reject(jets, [&](Jet* j) { return !selfn(j); }, do_delete);
105  }
void ifilter_reject(JetPtrs &jets, std::function< bool(Jet *)> rejfn, bool do_delete=true)
In-place filter a supplied jet vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:92
Here is the call graph for this function:

◆ iremoveerase()

template<typename CONTAINER , typename RMFN >
void Gambit::ColliderBit::iremoveerase ( CONTAINER &  c,
const RMFN &  fn 
)
inline

Convenience combination of remove_if and erase.

Definition at line 47 of file Utils.hpp.

Referenced by ifilter_reject().

47  {
48  auto newend = std::remove_if(c.begin(), c.end(), fn);
49  c.erase(newend, c.end());
50  }
Here is the caller graph for this function:

◆ is_xsec_sane()

bool Gambit::ColliderBit::is_xsec_sane ( const triplet< double > &  xsecWithError)

LEP limit debugging function.

Definition at line 96 of file ColliderBit_LEP.cpp.

References Gambit::triplet< TYPE >::central, Gambit::triplet< TYPE >::lower, and Gambit::triplet< TYPE >::upper.

Referenced by LEP188_SLHA1_convention_xsec_chi00_11(), LEP188_SLHA1_convention_xsec_chi00_12(), LEP188_SLHA1_convention_xsec_chi00_13(), LEP188_SLHA1_convention_xsec_chi00_14(), LEP188_SLHA1_convention_xsec_chi00_22(), LEP188_SLHA1_convention_xsec_chi00_23(), LEP188_SLHA1_convention_xsec_chi00_24(), LEP188_SLHA1_convention_xsec_chi00_33(), LEP188_SLHA1_convention_xsec_chi00_34(), LEP188_SLHA1_convention_xsec_chi00_44(), LEP188_SLHA1_convention_xsec_chipm_11(), LEP188_SLHA1_convention_xsec_chipm_12(), LEP188_SLHA1_convention_xsec_chipm_21(), LEP188_SLHA1_convention_xsec_chipm_22(), LEP188_SLHA1_convention_xsec_se1se1bar(), LEP188_SLHA1_convention_xsec_se1se2bar(), LEP188_SLHA1_convention_xsec_se2se1bar(), LEP188_SLHA1_convention_xsec_se2se2bar(), LEP188_SLHA1_convention_xsec_selselbar(), LEP188_SLHA1_convention_xsec_selserbar(), LEP188_SLHA1_convention_xsec_serselbar(), LEP188_SLHA1_convention_xsec_serserbar(), LEP188_SLHA1_convention_xsec_smu1smu1bar(), LEP188_SLHA1_convention_xsec_smu1smu2bar(), LEP188_SLHA1_convention_xsec_smu2smu1bar(), LEP188_SLHA1_convention_xsec_smu2smu2bar(), LEP188_SLHA1_convention_xsec_smulsmulbar(), LEP188_SLHA1_convention_xsec_smulsmurbar(), LEP188_SLHA1_convention_xsec_smursmulbar(), LEP188_SLHA1_convention_xsec_smursmurbar(), LEP188_SLHA1_convention_xsec_stau1stau1bar(), LEP188_SLHA1_convention_xsec_stau1stau2bar(), LEP188_SLHA1_convention_xsec_stau2stau1bar(), LEP188_SLHA1_convention_xsec_stau2stau2bar(), LEP188_SLHA1_convention_xsec_staulstaulbar(), LEP188_SLHA1_convention_xsec_staulstaurbar(), LEP188_SLHA1_convention_xsec_staurstaulbar(), LEP188_SLHA1_convention_xsec_staurstaurbar(), LEP205_SLHA1_convention_xsec_chi00_11(), LEP205_SLHA1_convention_xsec_chi00_12(), LEP205_SLHA1_convention_xsec_chi00_13(), LEP205_SLHA1_convention_xsec_chi00_14(), LEP205_SLHA1_convention_xsec_chi00_22(), LEP205_SLHA1_convention_xsec_chi00_23(), LEP205_SLHA1_convention_xsec_chi00_24(), LEP205_SLHA1_convention_xsec_chi00_33(), LEP205_SLHA1_convention_xsec_chi00_34(), LEP205_SLHA1_convention_xsec_chi00_44(), LEP205_SLHA1_convention_xsec_chipm_11(), LEP205_SLHA1_convention_xsec_chipm_12(), LEP205_SLHA1_convention_xsec_chipm_21(), LEP205_SLHA1_convention_xsec_chipm_22(), LEP205_SLHA1_convention_xsec_se1se1bar(), LEP205_SLHA1_convention_xsec_se1se2bar(), LEP205_SLHA1_convention_xsec_se2se1bar(), LEP205_SLHA1_convention_xsec_se2se2bar(), LEP205_SLHA1_convention_xsec_selselbar(), LEP205_SLHA1_convention_xsec_selserbar(), LEP205_SLHA1_convention_xsec_serselbar(), LEP205_SLHA1_convention_xsec_serserbar(), LEP205_SLHA1_convention_xsec_smu1smu1bar(), LEP205_SLHA1_convention_xsec_smu1smu2bar(), LEP205_SLHA1_convention_xsec_smu2smu1bar(), LEP205_SLHA1_convention_xsec_smu2smu2bar(), LEP205_SLHA1_convention_xsec_smulsmulbar(), LEP205_SLHA1_convention_xsec_smulsmurbar(), LEP205_SLHA1_convention_xsec_smursmulbar(), LEP205_SLHA1_convention_xsec_smursmurbar(), LEP205_SLHA1_convention_xsec_stau1stau1bar(), LEP205_SLHA1_convention_xsec_stau1stau2bar(), LEP205_SLHA1_convention_xsec_stau2stau1bar(), LEP205_SLHA1_convention_xsec_stau2stau2bar(), LEP205_SLHA1_convention_xsec_staulstaulbar(), LEP205_SLHA1_convention_xsec_staulstaurbar(), LEP205_SLHA1_convention_xsec_staurstaulbar(), LEP205_SLHA1_convention_xsec_staurstaurbar(), LEP207_SLHA1_convention_xsec_chi00_11(), LEP208_SLHA1_convention_xsec_chi00_11(), LEP208_SLHA1_convention_xsec_chi00_12(), LEP208_SLHA1_convention_xsec_chi00_13(), LEP208_SLHA1_convention_xsec_chi00_14(), LEP208_SLHA1_convention_xsec_chi00_22(), LEP208_SLHA1_convention_xsec_chi00_23(), LEP208_SLHA1_convention_xsec_chi00_24(), LEP208_SLHA1_convention_xsec_chi00_33(), LEP208_SLHA1_convention_xsec_chi00_34(), LEP208_SLHA1_convention_xsec_chi00_44(), LEP208_SLHA1_convention_xsec_chipm_11(), LEP208_SLHA1_convention_xsec_chipm_12(), LEP208_SLHA1_convention_xsec_chipm_21(), LEP208_SLHA1_convention_xsec_chipm_22(), LEP208_SLHA1_convention_xsec_se1se1bar(), LEP208_SLHA1_convention_xsec_se1se2bar(), LEP208_SLHA1_convention_xsec_se2se1bar(), LEP208_SLHA1_convention_xsec_se2se2bar(), LEP208_SLHA1_convention_xsec_selselbar(), LEP208_SLHA1_convention_xsec_selserbar(), LEP208_SLHA1_convention_xsec_serselbar(), LEP208_SLHA1_convention_xsec_serserbar(), LEP208_SLHA1_convention_xsec_smu1smu1bar(), LEP208_SLHA1_convention_xsec_smu1smu2bar(), LEP208_SLHA1_convention_xsec_smu2smu1bar(), LEP208_SLHA1_convention_xsec_smu2smu2bar(), LEP208_SLHA1_convention_xsec_smulsmulbar(), LEP208_SLHA1_convention_xsec_smulsmurbar(), LEP208_SLHA1_convention_xsec_smursmulbar(), LEP208_SLHA1_convention_xsec_smursmurbar(), LEP208_SLHA1_convention_xsec_stau1stau1bar(), LEP208_SLHA1_convention_xsec_stau1stau2bar(), LEP208_SLHA1_convention_xsec_stau2stau1bar(), LEP208_SLHA1_convention_xsec_stau2stau2bar(), LEP208_SLHA1_convention_xsec_staulstaulbar(), LEP208_SLHA1_convention_xsec_staulstaurbar(), LEP208_SLHA1_convention_xsec_staurstaulbar(), and LEP208_SLHA1_convention_xsec_staurstaurbar().

97  {
98  double xsec = xsecWithError.central;
99  double dxsec_upper = xsecWithError.upper - xsecWithError.central;
100  double dxsec_lower = xsecWithError.central - xsecWithError.lower;
101  if (xsec < 0.0 or dxsec_upper < 0.0 or dxsec_lower < 0.0
102  or Utils::isnan(xsec) or Utils::isnan(dxsec_upper) or Utils::isnan(dxsec_lower))
103  {
104  cout << "xsec: " << xsec << endl;
105  cout << "dxsec_upper: " << dxsec_upper << endl;
106  cout << "dxsec_lower: " << dxsec_lower << endl;
107  return false;
108  }
109  return true;
110  }
Here is the caller graph for this function:

◆ isFinalB()

template<typename EventT >
bool Gambit::ColliderBit::isFinalB ( int  n,
const EventT &  evt 
)
inline

Definition at line 105 of file Py8Utils.hpp.

106  {
107  // Root particle is invalid
108  if (n == 0) return false;
109  // *This* particle must be a b or b-hadron
110  if (!MCUtils::PID::hasBottom(evt[n].id())) return false;
111  // Daughters must *not* include a b or b-hadron
112  for (int m : evt.daughterList(n)) {
113  if (MCUtils::PID::hasBottom(evt[m].id())) return false;
114  }
115  return true;
116  }

◆ isFinalLepton()

template<typename EventT >
bool Gambit::ColliderBit::isFinalLepton ( int  n,
const EventT &  evt 
)
inline

Definition at line 176 of file Py8Utils.hpp.

Referenced by Gambit::ColliderBit::BuckFast< EventT >::convertPartonEvent().

177  {
178  // Root particle is invalid
179  if (n == 0) return false;
180  const auto& p = evt[n];
181  // Check if it's a lepton at all (including taus and neutrinos) & early exit
182  if (!HEPUtils::in_closed_range(abs(p.id()), 11, 16)) return false;
183  // Must have no daughters
184  return evt.daughterList(n).empty();
185  }
Here is the caller graph for this function:

◆ isFinalParton()

template<typename EventT >
bool Gambit::ColliderBit::isFinalParton ( int  n,
const EventT &  evt 
)
inline

Definition at line 147 of file Py8Utils.hpp.

References isParton().

148  {
149  // Root particle is invalid
150  if (n == 0) return false;
151  // Check if it's a parton at all & early exit
152  if (!isParton(n, evt)) return false;
153  // Daughters must *not* be partons
154  for (int m : evt.daughterList(n)) {
155  if (m == 0) continue; // 0 shouldn't be possible here, but just to be sure...
156  if (isParton(m, evt)) return false;
157  }
158  return true;
159  }
bool isParton(int n, const EventT &evt)
Definition: Py8Utils.hpp:135
Here is the call graph for this function:

◆ isFinalPhoton()

template<typename EventT >
bool Gambit::ColliderBit::isFinalPhoton ( int  n,
const EventT &  evt 
)
inline

Definition at line 163 of file Py8Utils.hpp.

Referenced by Gambit::ColliderBit::BuckFast< EventT >::convertPartonEvent().

164  {
165  // Root particle is invalid
166  if (n == 0) return false;
167  const auto& p = evt[n];
168  // Check if it's a photon at all & early exit
169  if (p.id() != 22) return false;
170  // Must have no daughters
171  return evt.daughterList(n).empty();
172  }
Here is the caller graph for this function:

◆ isFinalTau()

template<typename EventT >
bool Gambit::ColliderBit::isFinalTau ( int  n,
const EventT &  evt 
)
inline

Definition at line 120 of file Py8Utils.hpp.

121  {
122  // Root particle is invalid
123  if (n == 0) return false;
124  // *This* particle must be a tau
125  if (abs(evt[n].id()) != 15) return false;
126  // Daughters must *not* include a tau
127  for (int m : evt.daughterList(n)) {
128  if (abs(evt[m].id()) == 15) return false;
129  }
130  return true;
131  }

◆ isParton()

template<typename EventT >
bool Gambit::ColliderBit::isParton ( int  n,
const EventT &  evt 
)
inline

Definition at line 135 of file Py8Utils.hpp.

Referenced by isFinalParton().

136  {
137  // Root particle is invalid
138  if (n == 0) return false;
139  // This particle must be a parton (could use Py8 P::isParton() + apid == 6?)
140  int apid = abs(evt[n].id());
141  if (!HEPUtils::in_closed_range(apid, 1, 6) && apid != 21) return false;
142  return true;
143  }
Here is the caller graph for this function:

◆ L3_Chargino_All_Channels_Conservative_LLike()

void Gambit::ColliderBit::L3_Chargino_All_Channels_Conservative_LLike ( double result)

Definition at line 1870 of file ColliderBit_LEP.cpp.

References Gambit::DecayTable::at(), Gambit::DecayTable::Entry::BF(), Gambit::triplet< TYPE >::central, decay_rates, Gambit::Spectrum::get(), L3_Chargino_All_Channels_Conservative_LLike, LEP188_xsec_chipm_11, LEP188_xsec_chipm_22, limit_LLike(), Gambit::ColliderBit::BaseLimitContainer::limitAverage(), Gambit::triplet< TYPE >::lower, MSSM_spectrum, Gambit::Par::Pole_Mass, Gambit::Scanner::pow(), and Gambit::triplet< TYPE >::upper.

1871  {
1873  using std::pow;
1874  using std::log;
1875 
1876  const Spectrum& spec = *Dep::MSSM_spectrum;
1877 
1878  const DecayTable& decays = *Dep::decay_rates;
1879  const double mass_neut1 = spec.get(Par::Pole_Mass,1000022, 0);
1880  const double mass_char1 = spec.get(Par::Pole_Mass,1000024, 0);
1881  const double mass_char2 = spec.get(Par::Pole_Mass,1000037, 0);
1882  const double mZ = spec.get(Par::Pole_Mass,23, 0);
1883  triplet<double> xsecWithError;
1884  double xsecLimit, totalBR;
1885 
1886  static const L3CharginoAllChannelsLimitAt188pt6GeV limitContainer;
1887  // #ifdef COLLIDERBIT_DEBUG
1888  // static bool dumped=false;
1889  // if(!dumped)
1890  // {
1891  // limitContainer.dumpPlotData(45., 100., 0., 100., mZ, "lepLimitPlanev2/L3CharginoAllChannelsLimitAt188pt6GeV.dump");
1892  // dumped=true;
1893  // }
1894  // #endif
1895 
1896  result = 0;
1897  // Due to the nature of the analysis details of the model independent limit in
1898  // the paper, the best we can do is to try these processes individually:
1899 
1900  // char1, neut1
1901  xsecLimit = limitContainer.limitAverage(mass_char1, mass_neut1, mZ);
1902 
1903  xsecWithError = *Dep::LEP188_xsec_chipm_11;
1904  // Total up all channels which look like W* decays
1905  totalBR = 0;
1906  totalBR += decays.at("~chi+_1").BF("~chi0_1", "W+");
1907  totalBR += decays.at("~chi+_1").BF("~chi0_1", "u", "dbar");
1908  totalBR += decays.at("~chi+_1").BF("~chi0_1", "c", "sbar");
1909  totalBR += decays.at("~chi+_1").BF("~chi0_1", "e+", "nu_e");
1910  totalBR += decays.at("~chi+_1").BF("~chi0_1", "mu+", "nu_mu");
1911  totalBR += decays.at("~chi+_1").BF("~chi0_1", "tau+", "nu_tau");
1912  xsecWithError.upper *= pow(totalBR, 2);
1913  xsecWithError.central *= pow(totalBR, 2);
1914  xsecWithError.lower *= pow(totalBR, 2);
1915 
1916  if (xsecWithError.central < xsecLimit)
1917  {
1918  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1919  }
1920  else
1921  {
1922  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1923  }
1924 
1925  // char2, neut1
1926  xsecLimit = limitContainer.limitAverage(mass_char2, mass_neut1, mZ);
1927 
1928  xsecWithError = *Dep::LEP188_xsec_chipm_22;
1929  // Total up all channels which look like W* decays
1930  totalBR = 0;
1931  totalBR += decays.at("~chi+_2").BF("~chi0_1", "W+");
1932  totalBR += decays.at("~chi+_2").BF("~chi0_1", "u", "dbar");
1933  totalBR += decays.at("~chi+_2").BF("~chi0_1", "c", "sbar");
1934  totalBR += decays.at("~chi+_2").BF("~chi0_1", "e+", "nu_e");
1935  totalBR += decays.at("~chi+_2").BF("~chi0_1", "mu+", "nu_mu");
1936  totalBR += decays.at("~chi+_2").BF("~chi0_1", "tau+", "nu_tau");
1937  xsecWithError.upper *= pow(totalBR, 2);
1938  xsecWithError.central *= pow(totalBR, 2);
1939  xsecWithError.lower *= pow(totalBR, 2);
1940 
1941  if (xsecWithError.central < xsecLimit)
1942  {
1943  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
1944  }
1945  else
1946  {
1947  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
1948  }
1949 
1950  }
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 DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable L3_Chargino_All_Channels_Conservative_LLike
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 LEP188_xsec_chipm_22
double limit_LLike(double x, double x95, double sigma)
LEP limit likelihood function.
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 LEP188_xsec_chipm_11
Spectrum Spectrum
double pow(const double &a)
Outputs a^i.
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 DecayTable::Entry DecayTable::Entry DecayTable::Entry decay_rates
Here is the call graph for this function:

◆ L3_Chargino_Leptonic_Conservative_LLike()

void Gambit::ColliderBit::L3_Chargino_Leptonic_Conservative_LLike ( double result)

Definition at line 1952 of file ColliderBit_LEP.cpp.

References Gambit::DecayTable::at(), Gambit::DecayTable::Entry::BF(), Gambit::triplet< TYPE >::central, decay_rates, Gambit::Spectrum::get(), LEP188_xsec_chipm_11, LEP188_xsec_chipm_22, limit_LLike(), Gambit::ColliderBit::BaseLimitContainer::limitAverage(), Gambit::triplet< TYPE >::lower, MSSM_spectrum, Gambit::Par::Pole_Mass, Gambit::Scanner::pow(), and Gambit::triplet< TYPE >::upper.

1953  {
1955  using std::pow;
1956  using std::log;
1957 
1958  const Spectrum& spec = *Dep::MSSM_spectrum;
1959 
1960  const DecayTable& decays = *Dep::decay_rates;
1961  const double mass_neut1 = spec.get(Par::Pole_Mass,1000022, 0);
1962  const double mass_char1 = spec.get(Par::Pole_Mass,1000024, 0);
1963  const double mass_char2 = spec.get(Par::Pole_Mass,1000037, 0);
1964  const double mZ = spec.get(Par::Pole_Mass,23, 0);
1965  triplet<double> xsecWithError;
1966  double xsecLimit, totalBR;
1967 
1968  static const L3CharginoLeptonicLimitAt188pt6GeV limitContainer;
1969  // #ifdef COLLIDERBIT_DEBUG
1970  // static bool dumped=false;
1971  // if(!dumped)
1972  // {
1973  // limitContainer.dumpPlotData(45., 100., 0., 100., mZ, "lepLimitPlanev2/L3CharginoLeptonicLimitAt188pt6GeV.dump");
1974  // dumped=true;
1975  // }
1976  // #endif
1977 
1978  result = 0;
1979  // Due to the nature of the analysis details of the model independent limit in
1980  // the paper, the best we can do is to try these processes individually:
1981 
1982  // char1, neut1
1983  xsecLimit = limitContainer.limitAverage(mass_char1, mass_neut1, mZ);
1984 
1985  xsecWithError = *Dep::LEP188_xsec_chipm_11;
1986  // Total up all channels which look like leptonic W* decays
1987  // Total up the leptonic W decays first...
1988  totalBR = 0;
1989  totalBR += decays.at("W+").BF("e+", "nu_e");
1990  totalBR += decays.at("W+").BF("mu+", "nu_mu");
1991  totalBR += decays.at("W+").BF("tau+", "nu_tau");
1992  totalBR = decays.at("~chi+_1").BF("~chi0_1", "W+") * totalBR;
1993 
1994  totalBR += decays.at("~chi+_1").BF("~chi0_1", "e+", "nu_e");
1995  totalBR += decays.at("~chi+_1").BF("~chi0_1", "mu+", "nu_mu");
1996  totalBR += decays.at("~chi+_1").BF("~chi0_1", "tau+", "nu_tau");
1997  xsecWithError.upper *= pow(totalBR, 2);
1998  xsecWithError.central *= pow(totalBR, 2);
1999  xsecWithError.lower *= pow(totalBR, 2);
2000 
2001  if (xsecWithError.central < xsecLimit)
2002  {
2003  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.upper - xsecWithError.central);
2004  }
2005  else
2006  {
2007  result += limit_LLike(xsecWithError.central, xsecLimit, xsecWithError.central - xsecWithError.lower);
2008  }
2009 
2010  // char2, neut1
2011  xsecLimit = limitContainer.limitAverage(mass_char2, mass_neut1, mZ);
2012 
2013  xsecWithError = *Dep::LEP188_xsec_chipm_22;
2014  // Total up all channels which look like leptonic W* decays
2015  // Total up the leptonic W decays first...
2016  totalBR = 0;
2017  totalBR += decays.at("W+").BF("e+", "nu_e");
2018  totalBR += decays.at("W+").BF("mu+", "nu_mu");
2019  totalBR += decays.at("W+").BF("tau+", "nu_tau");
2020  totalBR = decays.at("~chi+_2").BF("~chi0_1", "W+") * totalBR;
2021 
2022  totalBR += decays.at("~chi+_2").BF("~chi0_1", "e+", "nu_e");
2023  totalBR += decays.at("~chi+_2").BF("~chi0_1", "mu+", "nu_mu");
2024  totalBR += decays.at("~chi+_2").BF("~chi0_1", "tau+", "nu_tau");
2025  xsecWithError.upper *= pow(totalBR, 2);
2026  xsecWithError.central *= pow(totalBR, 2);
2027  xsecWithError.lower *= pow(totalBR, 2);
2028 
2029  if (xsecWithError.central < xsecLimit)
2030  {