ATOOLS/Phys/Variations.H

0001 #ifndef ATOOLS_Phys_Variations_H
0002 #define ATOOLS_Phys_Variations_H
0003
0004 #include <ostream>
0005 #include <map>
0006 #include <string>
0007 #include <vector>
0008
0009 #include "ATOOLS/Org/Exception.H"
0010 #include "ATOOLS/Org/Shell_Tools.H"
0011 #include "ATOOLS/Org/Enum_Flags.H"
0012 #include "ATOOLS/Org/Scoped_Settings.H"
0013
0014 #define ENABLE_REWEIGHTING_FACTORS_HISTOGRAMS 0
0015
0016 namespace PDF { class PDF_Base; }
0017 namespace MODEL { class Running_AlphaS; }
0018 namespace BEAM { class Beam_Spectra_Handler; }
0019
0020
0021 namespace ATOOLS {
0022
0023   //! Constants to denote various variation sources
0024   enum class Variations_Source {
0025     all,     //!< use for total weight retrieval, i.e. the product of the below
0026     main,    //!< use for everything that is not one of the below
0027     sudakov  //!< use for weight factors from the reweighting of the Sudakov
0028              //!< Veto Algorithm, i.e. the reweighting of parton shower
0029              //!< emissions
0030   };
0031   std::ostream& operator<<(std::ostream&, const Variations_Source&);
0032
0033   enum class Variations_Type { qcd, qcut, custom };
0034   enum class Variations_Mode { all, nominal_only };
0035   std::ostream& operator<<(std::ostream&, const Variations_Type&);
0036
0037   enum class Variations_Name_Type { weight_name_convention, human_readable };
0038
0039   class Variations;
0040   struct QCD_Variation_Params;
0041   struct Qcut_Variation_Params;
0042
0043   extern Variations* s_variations;
0044
0045   struct ScaleFactorExpansions {
0046     enum code {
0047       None            = 0,
0048       MuF             = 1,
0049       MuR             = 1<<1,
0050       SevenPoint      = 1<<2,  // vary independently, but exclude (up,down) and (down,up) variations
0051       QCUT            = 1<<3
0052     };
0053   };
0054   DEFINE_SHERPA_ENUM_FLAG_OPERATORS(ScaleFactorExpansions::code)
0055
0056   struct NonfactorizingCoupling {
0057     enum code {
0058       WithCustomVariationWeight =
0059           98, // coupling it is not factorizing, but
0060               // a custom calculator function for a variation weight is provided
0061       WithoutCustomVariationWeight =
0062           99, // coupling is not factorizing, and there is no calculator
0063               // function for a variation weight provided, i.e. no on-the-fly
0064               // reweighting can not be done
0065     };
0066   };
0067
0068   struct Qcut_Variation_Params {
0069     Qcut_Variation_Params(double scale_factor):
0070       m_scale_factor{scale_factor}
0071     {}
0072     bool IsTrivial() const { return (m_scale_factor == 1.0); }
0073     std::string Name(Variations_Source source=ATOOLS::Variations_Source::all) const;
0074     double m_scale_factor;
0075   };
0076
0077   //! Initialise and hold all information that is needed for varying input parameters
0078   class Variations {
0079   public:
0080     static bool NeedsLHAPDF6Interface();
0081     static void CheckConsistencyWithBeamSpectra(BEAM::Beam_Spectra_Handler *);
0082
0083     Variations(Variations_Mode mode=Variations_Mode::all);
0084     ~Variations();
0085
0086     void EnableVariations() { m_enabled = true; }
0087     void DisableVariations() { m_enabled = false; }
0088
0089     std::vector<Variations_Type> ManagedVariationTypes() const
0090     {
0091       return {Variations_Type::qcd, Variations_Type::qcut};
0092     }
0093
0094     typedef std::vector<QCD_Variation_Params *> Parameters_Vector;
0095     const Parameters_Vector * GetParametersVector() const { return &m_parameters_vector; }
0096     QCD_Variation_Params& Parameters(size_t i) { return *m_parameters_vector[i]; }
0097     Qcut_Variation_Params& Qcut_Parameters(size_t i) { return m_qcut_parameters_vector[i]; }
0098     std::string
0099     GetVariationNameAt(Variations::Parameters_Vector::size_type,
0100                        Variations_Type type = Variations_Type::qcd,
0101                        Variations_Source source = Variations_Source::all,
0102                        Variations_Name_Type name_type =
0103                            Variations_Name_Type::weight_name_convention) const;
0104     size_t Size(Variations_Type t=Variations_Type::qcd) const;
0105     bool HasVariations() const {
0106       for (const auto t : ManagedVariationTypes())
0107         if (Size(t) > 0)
0108           return true;
0109       return false;
0110     }
0111     template <typename Handler>
0112     void ForEach(Handler h);
0113
0114     // register and print warnings
0115     void IncrementOrInitialiseWarningCounter(const std::string name) { m_warnings[name]++; }
0116     void PrintStatistics(std::ostream &str);
0117
0118   private:
0119     void ReadDefaults();
0120     void LoadLHAPDFInterfaceIfNecessary();
0121     void InitialiseParametersVector();
0122     void AddParameters(Scoped_Settings&);
0123     struct PDFs_And_AlphaS {
0124       PDFs_And_AlphaS();
0125       PDFs_And_AlphaS(double alphasmz);
0126       PDFs_And_AlphaS(std::string pdfname, int pdfmember, int beammask,
0127                       int alphasbeam);
0128       std::vector<PDF::PDF_Base *> m_pdfs;
0129       MODEL::Running_AlphaS *p_alphas;
0130       int m_shoulddeletepdfmask {0};
0131       bool m_shoulddeletealphas {false};
0132     };
0133     struct PDFs_And_AlphaS_List {
0134       std::vector<PDFs_And_AlphaS> items;
0135       bool did_expand {false};
0136     };
0137     void AddParameterExpandingScaleFactors(std::vector<std::string> scalestringparams,
0138                                            ScaleFactorExpansions::code,
0139                                            PDFs_And_AlphaS_List);
0140     void AddParameters(double, double,
0141                        std::vector<PDFs_And_AlphaS>::const_iterator,
0142                        int deletepdfmask, bool deletealphas);
0143     PDFs_And_AlphaS_List PDFsAndAlphaSList(const std::string& pdfstringparam) const;
0144     PDFs_And_AlphaS_List PDFsAndAlphaSList(std::string pdfstringparam,
0145                                            bool expandpdf, int beammask,
0146                                            int alphasbeam) const;
0147
0148     bool m_enabled;
0149     Parameters_Vector m_parameters_vector;
0150     std::vector<Qcut_Variation_Params> m_qcut_parameters_vector;
0151
0152     //! (name -> counter value) map used to track warnings
0153     std::map<std::string, unsigned long> m_warnings;
0154
0155     //! whether the central value should be included when expanding VARIATIONS arguments
0156     bool m_includecentralvaluevariation;
0157     //! whether the muR variation factor should be applied to the shower AlphaS scale arguments
0158     bool m_reweightsplittingalphasscales;
0159     //! whether the muF variation factor should be applied to the shower PDF scale arguments
0160     bool m_reweightsplittingpdfsscales;
0161   };
0162
0163   template <typename Handler> void Variations::ForEach(Handler h)
0164   {
0165     const auto size = Size();
0166     for (int i {0}; i < size; ++i) {
0167       h(i, Parameters(i));
0168     }
0169   }
0170
0171   std::ostream& operator<<(std::ostream&, const Variations&);
0172
0173 #if ENABLE_REWEIGHTING_FACTORS_HISTOGRAMS
0174   class ReweightingFactorHistogram {
0175
0176   public:
0177
0178     ReweightingFactorHistogram();
0179
0180     void Fill(std::string, double);
0181     void Write(std::string filenameaffix);
0182
0183   private:
0184
0185     typedef std::string EntryKey;
0186     typedef std::vector<long> EntryNumbers;
0187     typedef std::pair<double, EntryNumbers> Bin;
0188
0189     std::vector<EntryKey> m_keys;
0190     std::vector<Bin> m_bins;
0191   };
0192 #endif
0193
0194   /*!
0195    * Hold a set of input parameters and factors for a single input parameter variation
0196    *
0197    * Used as an argument to a reweighting function.
0198    */
0199   struct QCD_Variation_Params {
0200     QCD_Variation_Params(double muR2fac, double muF2fac,
0201                          bool showermuR2enabled,
0202                          bool showermuF2enabled,
0203                          PDF::PDF_Base *pdf1,
0204                          PDF::PDF_Base *pdf2,
0205                          MODEL::Running_AlphaS *alphas,
0206                          int deletepdfmask,
0207                          bool deletealphas):
0208       m_muR2fac(muR2fac), m_muF2fac(muF2fac),
0209       m_showermuR2enabled(showermuR2enabled),
0210       m_showermuF2enabled(showermuF2enabled),
0211       p_pdf1(pdf1), p_pdf2(pdf2), p_alphas(alphas),
0212       m_deletepdfmask(deletepdfmask),
0213       m_deletealphas(deletealphas)
0214     {}
0215     ~QCD_Variation_Params();
0216
0217     void IncrementOrInitialiseWarningCounter(const std::string name) { m_warnings[name]++; }
0218
0219 #if ENABLE_REWEIGHTING_FACTORS_HISTOGRAMS
0220     void FillReweightingFactorsHisto(std::string name, double value) { m_rewfachisto.Fill(name, value); }
0221 #endif
0222
0223     bool IsTrivial() const;
0224     std::string
0225     Name(Variations_Source source = ATOOLS::Variations_Source::all,
0226          Variations_Name_Type name_type =
0227              ATOOLS::Variations_Name_Type::weight_name_convention) const;
0228
0229     const double m_muR2fac, m_muF2fac;
0230     const double m_showermuR2enabled, m_showermuF2enabled;
0231     //! Pointers to the beam PDFs, which can be NULL for (semi-)leptonic events
0232     PDF::PDF_Base * const p_pdf1;
0233     PDF::PDF_Base * const p_pdf2;
0234     MODEL::Running_AlphaS * const p_alphas;
0235     //! Set whether the pointers to the PDFs and the AlphaS should be deleted after usage
0236     const int m_deletepdfmask;
0237     const bool m_deletealphas;
0238
0239   private:
0240     //! (name -> counter value) map used to track warnings
0241     std::map<std::string, unsigned long> m_warnings;
0242     friend class Variations;
0243 #if ENABLE_REWEIGHTING_FACTORS_HISTOGRAMS
0244     ReweightingFactorHistogram m_rewfachisto;
0245 #endif
0246   };
0247
0248   template <typename U>
0249   std::string GenerateVariationNamePart(std::string tag, U value) {
0250     return tag + ToString(value);
0251   }
0252
0253   bool IsQCDVariationTrivial(
0254       double muR2fac, double muF2fac,
0255       PDF::PDF_Base * const pdf1,
0256       PDF::PDF_Base * const pdf2,
0257       MODEL::Running_AlphaS * const);
0258
0259 };
0260
0261 #endif // #ifndef ATOOLS_Phys_Variations_H