Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:06:36

0001 // Weights.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // This file contains classes that keep track of event weights.
0007 
0008 #ifndef Pythia8_Weights_H
0009 #define Pythia8_Weights_H
0010 
0011 #include "Pythia8/Basics.h"
0012 #include "Pythia8/LHEF3.h"
0013 #include "Pythia8/PythiaStdlib.h"
0014 #include "Pythia8/SharedPointers.h"
0015 
0016 namespace Pythia8 {
0017 
0018 // Forward declare classes which need internal member access.
0019 class Info;
0020 class History;
0021 class PartonLevel;
0022 class Merging;
0023 class WeightContainer;
0024 
0025 //==========================================================================
0026 
0027 // This is a base class to store weight information in a way that allows
0028 // unified access in the structure that contains all event generation weights
0029 // information (WeightContainer below). The main purpuse of this class is to
0030 // supply convenience features to make defining new categories of weights easy.
0031 // All weights should inherit from this base class. The specialized classes
0032 // may then contain specialized functions, or only act as a glorified
0033 // book-keeping structure.
0034 
0035 class WeightsBase {
0036 
0037 public:
0038 
0039   // Friends for internal protected members.
0040   friend class History;
0041   friend class PartonLevel;
0042 
0043   // Initialize the weights.
0044   virtual void init() {weightValues.resize(0); weightNames.resize(0);
0045     bookWeight("Baseline");}
0046   virtual void init(bool) {}
0047 
0048   // Reset all internal values.
0049   virtual void clear() {fill(weightValues.begin(), weightValues.end(), 1.);}
0050 
0051   // Store the current event information (replace whitespace with underscore
0052   // for HepMC).
0053   virtual void bookVectors(vector<double> weights, vector<string> names) {
0054     for (int i = 0; i < (int)weights.size(); ++i) {
0055       replace(names[i].begin(), names[i].end(), ' ', '_');
0056       bookWeight(names[i], weights[i]);
0057     }
0058   }
0059 
0060   // Function to return processed weights to weight container, e.g. if
0061   // weights should be combined before proceeding.
0062   virtual void collectWeightValues(vector<double>& outputWeights,
0063     double norm = 1.);
0064 
0065   // Similar function to return processed weight names.
0066   virtual void collectWeightNames(vector<string>& outputNames);
0067 
0068   // Get the stored information.
0069   // Direcly use storage members here in the base class,
0070   // and access those through non-virtual getters.
0071   // Note: NOT opting for a map data structure, since information will anyway
0072   // have to be serialized in output.
0073   // Change weight names for compatibility with Rivet by replacing colons
0074   // with full stops.
0075   string getWeightsName(int iPos) const  {
0076     string name = iPos >= 0
0077       && iPos < (int)weightNames.size() ? weightNames[iPos] : "";
0078     if (name.find(":") != string::npos)
0079       replace(name.begin(), name.end(), ':', '.');
0080     return name == "" ? to_string(iPos) : name;}
0081   virtual double getWeightsValue(int iPos) const { return weightValues[iPos]; }
0082   int getWeightsSize() const { return weightValues.size(); }
0083 
0084   // Function to create a new, synchronized, pair of weight name and value.
0085   void bookWeight(string name, double defaultValue = 1.)  {
0086     // Check if name already exists.
0087     if (findIndexOfName(name) != -1) setValueByName(name, defaultValue);
0088     else {
0089       weightNames.push_back(name);
0090       weightValues.push_back(defaultValue);
0091     }
0092   }
0093 
0094   // Functions to set values of weights.
0095   void setValueByIndex(int iPos, double val) {
0096     if (iPos < 0 || iPos >= (int)weightValues.size()) return;
0097     weightValues[iPos] = val;
0098   }
0099   void setValueByName(string name, double val) {
0100     setValueByIndex(findIndexOfName(name), val);}
0101 
0102   // Functions to reweight weight values, by index or by name.
0103   virtual void reweightValueByIndex(int iPos, double val) {
0104     if (iPos < 0 || iPos >= (int)weightValues.size()) return;
0105     weightValues[iPos] *= val;
0106   };
0107   virtual void reweightValueByName(string name, double val) {
0108     // Use existing functions: Find index of name, then set by index.
0109     int iPos = findIndexOfName(name);
0110     reweightValueByIndex(iPos, val);
0111   }
0112 
0113   // Function to find the index of a given entry in the weightNames vector,
0114   // e.g., to be able to synchronize with the weightValues vector.
0115   int findIndexOfName(string name) {
0116     vector<string>::iterator it
0117       = find(weightNames.begin(), weightNames.end(), name);
0118     unsigned long int index = distance(weightNames.begin(), it);
0119     if (index == weightNames.size()) return -1;
0120     return distance(weightNames.begin(), it);
0121   }
0122 
0123   // Set the pointers.
0124   void setPtrs(Info* infoPtrIn) { infoPtr = infoPtrIn; }
0125 
0126 protected:
0127 
0128   // Parse a WVec of variations into a dictionary.
0129   void parse(string wvecKey, map<string, map<string, double> > &dct);
0130 
0131   // Weight values and names.
0132   vector<double> weightValues;
0133   vector<string> weightNames;
0134 
0135   // Pointers necessary for variation initialization
0136   Info* infoPtr{};
0137 
0138 };
0139 
0140 //==========================================================================
0141 
0142 // Purely virtual base class for shower weights.
0143 
0144 class WeightsShower : public WeightsBase {
0145 
0146 public:
0147 
0148   // Initialize weights (more can be booked at any time)
0149   void init(bool) override {}
0150 
0151   // Weight "groups" (combinations of one or more unique weights).
0152   virtual void initWeightGroups(bool = false) {}
0153   virtual int nWeightGroups() const {return 0;}
0154   virtual string getGroupName( int /*iGroup*/ ) const {return "";}
0155   virtual double getGroupWeight( int /*iGroup*/ ) const {return 1.;}
0156 
0157 };
0158 
0159 //==========================================================================
0160 
0161 // This shows a WeightsShower example implementation for SimpleShower.
0162 
0163 class WeightsSimpleShower : public WeightsShower {
0164 
0165 public:
0166 
0167   // Initialize weights (more can be booked at any time)
0168   void init( bool doMerging) override;
0169 
0170   // Functions to return processed weights to weight container, e.g. if
0171   // weights should be combined before proceeding.
0172   void collectWeightNames(vector<string>& outputNames) override;
0173   void collectWeightValues(vector<double>& outputWeights,
0174     double norm = 1.) override;
0175 
0176   // Initialize the weight groups for automated variations.
0177   void initWeightGroups(bool = false) override;
0178 
0179   // Return group name (want to integrate this in weightNameVector?)
0180   string getGroupName(int iGN) const override;
0181 
0182   // Return group weight (want to integrate this in weightValueVector?)
0183   double getGroupWeight(int iGW) const override;
0184   int    nWeightGroups() const override { return externalVariations.size(); }
0185 
0186   // Initialise list of atomic weight variations to be performed by shower.
0187   bool initUniqueShowerVars();
0188 
0189   // Memorize enhancement factor and pT of enhanced emission
0190   void setEnhancedTrial( double pTIn, double wtIn) { pTEnhanced = pTIn;
0191     wtEnhanced = wtIn; }
0192   double getEnhancedTrialPT() { return pTEnhanced;}
0193   double getEnhancedTrialWeight() { return wtEnhanced;}
0194 
0195   // Extract full list or subset that matches specific keys (e.g., FSR ones).
0196   vector<string> getUniqueShowerVars() { return uniqueShowerVars; }
0197   vector<string> getUniqueShowerVars(vector<string> keys);
0198 
0199   // Initialise map of enhancement factors for shower trial probabilities.
0200   bool initEnhanceFactors();
0201   unordered_map<string,double> getEnhanceFactors() { return enhanceFactors; }
0202 
0203   string getInitialName(int iG) const { return initialNameSave[iG]; }
0204 
0205   // Variations that must be known by TimeShower and Spaceshower.
0206   map<int,double> varPDFplus, varPDFminus, varPDFmember;
0207 
0208   // Vectors for storing shower variatons and enhancement factors.
0209   vector<string> uniqueShowerVars;
0210   unordered_map<string,double> enhanceFactors;
0211 
0212   // Vectors for weight group handling
0213   vector<string>          externalVariations;
0214   vector<vector<string> > externalVarNames;
0215   vector<string>          externalGroupNames;
0216   vector<string>          initialNameSave;
0217   vector<vector<int> >    externalMap;
0218   int                     externalVariationsSize{};
0219 
0220   // Vector for merging requested weight handling
0221   vector<double> getMuRWeightVector();
0222   vector<vector<string> > mergingVarNames;
0223   double pTEnhanced{}, wtEnhanced{};
0224 
0225 };
0226 
0227 //==========================================================================
0228 
0229 // This class collects information on weights generated in the merging
0230 // framework. The first weight is always required for CKKW-L, UMEPS and
0231 // NLO merging. The additional weights are required for simultaneous
0232 // renormalization scale variation in matrix element generation and parton
0233 // shower.
0234 
0235 class WeightsMerging : public WeightsBase {
0236 
0237 public:
0238 
0239   // Friends for internal protected members.
0240   friend class Merging;
0241   friend class WeightContainer;
0242 
0243   // Initialize weights (more can be booked at any time)
0244   void init() override;
0245 
0246   // Reset all internal values;
0247   void clear() override;
0248 
0249   // Function to create a new, synchronized, pair of weight name and value.
0250   void bookWeight(string name, double value, double valueFirst);
0251 
0252   // Store the current event information.
0253   void bookVectors(vector<double> weights, vector<string> names) override;
0254   void bookVectors(vector<double> weights, vector<double> weightsFirst,
0255             vector<string> names);
0256   // Modified weight getter to include first order weight
0257   double getWeightsValue(int iPos) const override {
0258     return weightValues[iPos] - weightValuesFirst[iPos]; }
0259   // Also add getters for UNLOPS-P and -PC schemes
0260   double getWeightsValueP(int iPos) const {
0261     return weightValuesP[iPos] - weightValuesFirstP[iPos]; }
0262   double getWeightsValuePC(int iPos) const {
0263     return weightValuesPC[iPos] - weightValuesFirstPC[iPos]; }
0264 
0265   // Functions to set values of weights.
0266   void reweightValueByIndex(int iPos, double val) override;
0267   void reweightValueByName(string name, double val) override;
0268 
0269   // Functions to set values of first order weights.
0270   void setValueFirstByIndex(int iPos, double val);
0271   void setValueFirstByName(string name, double val);
0272 
0273   // Functions to set values as whole vector.
0274   void setValueVector(vector<double> ValueVector);
0275   void setValueFirstVector(vector<double> ValueFirstVector);
0276 
0277   // Function telling merging which muR variations to perform
0278   vector<double> getMuRVarFactors();
0279 
0280   // Set up mapping between LHEF variations.
0281   void setLHEFvariationMapping();
0282 
0283   // Function to collect weight names.
0284   void collectWeightNames(vector<string>& outputNames) override;
0285 
0286   // Function collecting weight values.
0287   void collectWeightValues(vector<double>& outputWeights,
0288      double norm = 1.) override;
0289 
0290 protected:
0291 
0292   // Corresponding vector with respective LHEF weight indices.
0293   map<int,int> muRVarLHEFindex;
0294 
0295   // Data member for first order weight.
0296   vector<double> weightValuesFirst;
0297 
0298   // Data members for UNLOPS-P and -PC.
0299   vector<double> weightValuesP, weightValuesPC,
0300     weightValuesFirstP, weightValuesFirstPC;
0301 
0302   // Boolean to memorize if LHEF weight needs to be applied (only for NLO).
0303   bool isNLO;
0304 
0305 };
0306 
0307 //==========================================================================
0308 
0309 // This is a short example class to collect information on Les Houches
0310 // Event weights into a container class that can be part of Weight,
0311 // which in turn is part of Info.
0312 
0313 class WeightsLHEF : public WeightsBase {
0314 
0315 public:
0316 
0317   // Friends for internal protected members.
0318   friend class WeightsMerging;
0319 
0320   // Central weight, needed for normalization, set from ProcessContainer.cc
0321   double centralWeight;
0322 
0323   // Reset all internal values;
0324   void clear() override;
0325 
0326   // Store the current event information.
0327   void bookVectors(vector<double> weights, vector<string> names) override;
0328 
0329   // Function to return processed weights to weight container, e.g. if
0330   // weights should be combined before proceeding.
0331   void collectWeightValues(vector<double>& outputWeights,
0332      double norm = 1.) override;
0333   void collectWeightNames(vector<string>& outputNames) override;
0334 
0335   // Convert weight names in MadGraph5 convention to the convention outlined
0336   // in https://arxiv.org/pdf/1405.1067.pdf, page 162ff.
0337   vector<string> convertNames(vector<string> names);
0338 
0339   void identifyVariationsFromLHAinit(map<string,LHAweight>* weights);
0340 
0341 protected:
0342 
0343   // Renormalization variations.
0344   map<int,double> muRvars;
0345 
0346 };
0347 
0348 //==========================================================================
0349 
0350 // This is class to collect information on fragmentation weighting,
0351 // and is in turn part of Info.
0352 
0353 class WeightsFragmentation : public WeightsBase {
0354 
0355 public:
0356 
0357   // Initialize the weights.
0358   void init() override;
0359 
0360   int nWeightGroups() const {return externalGroupNames.size();}
0361 
0362   string getGroupName(int iGN) const {return iGN < 0 || iGN >= nWeightGroups()
0363       ? "Null" : externalGroupNames[iGN];}
0364 
0365   // Return group weight.
0366   double getGroupWeight(int iGW) const {
0367     if (iGW < 0 || iGW >= nWeightGroups()) return 1.0;
0368     double wgt(1.);
0369     for (const int &iWgt : externalMap[iGW]) wgt *= getWeightsValue(iWgt);
0370     return wgt;}
0371 
0372   // Functions to return processed weights to weight container, e.g. if
0373   // weights should be combined before proceeding.
0374   void collectWeightNames(vector<string>& outputNames) override;
0375   void collectWeightValues(vector<double>& outputWeights,
0376     double norm = 1.) override;
0377 
0378   // Vectors for weight group handling.
0379   vector<map<vector<double>, int> > weightParms{};
0380   vector<string>       externalGroupNames{};
0381   vector<vector<int> > externalMap{};
0382 
0383   // Factorization indices.
0384   enum FactIndex{Z, Flav, PT};
0385 
0386 };
0387 
0388 //==========================================================================
0389 
0390 // This is a container class to collect all event generation weight
0391 // information into a wrapper which is in turn is part of Info. In this
0392 // way, we could avoid cluttering Info.
0393 
0394 class WeightContainer {
0395 
0396 public:
0397 
0398   // Default constructor only ensures that members are initialized with
0399   // sensible default values.
0400   WeightContainer() : weightNominal(1.0), xsecIsInit(false) {}
0401 
0402   // The nominal Pythia weight, in pb for lha strategy 4 and -4
0403   double weightNominal;
0404   void setWeightNominal( double weightNow );
0405   double collectWeightNominal();
0406 
0407   // First example of a weight subcategory.
0408   WeightsLHEF          weightsLHEF{};
0409 
0410   // Shower weights. SimpleShower one hardcoded and made default since
0411   // Merging explicitly depends on it, but allow to point to a
0412   // different weightsShower-derived object, eg for Vincia (with own
0413   // merging).
0414   WeightsShower*       weightsShowerPtr{};
0415   WeightsSimpleShower  weightsSimpleShower{};
0416 
0417   // Merging weights.
0418   WeightsMerging       weightsMerging{};
0419 
0420   // Fragmentation weights.
0421   WeightsFragmentation weightsFragmentation{};
0422 
0423   // Userhooks weights.
0424   WeightsBase          weightsUserHooks{};
0425 
0426   // Functions to retrieve information stored in the subcategory members.
0427   int numberOfWeights();
0428   double weightValueByIndex(int key = 0);
0429   string weightNameByIndex(int key = 0);
0430 
0431   // Function to return the vector of weight values, combining all
0432   // weights from all subcategories. Currently, only the nominal
0433   // weight and LHEF weights are considered. Internal Pythia weights
0434   // should also be included eventually.
0435   vector<double> weightValueVector();
0436 
0437   // Function to return the vector of weight names, combining all names from
0438   // all subcategories, cf. weightValueVector function.
0439   vector<string> weightNameVector();
0440 
0441   // Reset all members to default stage.
0442   void clear();
0443 
0444   // Reset total cross section estimate.
0445   void clearTotal();
0446 
0447   // Init, for those classes that need it.
0448   void init(bool doMerging);
0449 
0450   // Function to set Pointers in weight classes.
0451   void initPtrs(Info* infoPtrIn);
0452 
0453   // Initialize the cross-section vector.
0454   void initXsecVec();
0455 
0456   // Get cross-section vectors.
0457   vector<double> getSampleXsec();
0458   vector<double> getTotalXsec();
0459   vector<double> getSampleXsecErr();
0460   vector<double> getTotalXsecErr();
0461 
0462   // Accumulate cross section for all weights.
0463   void accumulateXsec(double norm = 1.);
0464 
0465 private:
0466 
0467   // Pointers necessary for variation initialization.
0468   Info* infoPtr{};
0469 
0470   // Suppress AUX_ weights.
0471   bool doSuppressAUXweights{};
0472 
0473   // Vectors for cross-section.
0474   vector<double> sigmaTotal, sigmaSample, errorTotal, errorSample;
0475   bool xsecIsInit;
0476 
0477 };
0478 
0479 //==========================================================================
0480 
0481 } // end namespace Pythia8
0482 
0483 #endif // Pythia8_Weights_H