include/Pythia8/DireWeightContainer.h

0001 // DireWeightContainer.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Stefan Prestel, 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 // Header file for weight constainers for Dire parton shower.
0007
0008 #ifndef Pythia8_DireWeightContainer_H
0009 #define Pythia8_DireWeightContainer_H
0010
0011 #define DIRE_WEIGHTCONTAINER_VERSION "2.002"
0012
0013 #include "Pythia8/PythiaStdlib.h"
0014 #include "Pythia8/Settings.h"
0015 #include "Pythia8/DireBasics.h"
0016 #include "Pythia8/ExternalMEs.h"
0017
0018 namespace Pythia8 {
0019
0020 //==========================================================================
0021
0022 // Container for a single weight with auxiliary information.
0023
0024 class DirePSWeight {
0025
0026 public:
0027
0028   // Constructors.
0029   DirePSWeight()
0030     : wt(1.0), type(0), iAtt(0), dAtt(0.0), sAtt("") {}
0031
0032   DirePSWeight( double w)
0033     : wt(w), type(0), iAtt(0), dAtt(0.0), sAtt("") {}
0034
0035   DirePSWeight( double w, int typeIn, int iAttIn=0, double dAttIn=0.0,
0036     string sAttIn="") : wt(w), type(typeIn), iAtt(iAttIn), dAtt(dAttIn),
0037     sAtt(sAttIn) {}
0038
0039   DirePSWeight( const DirePSWeight& wgt) : wt(wgt.wt), type(wgt.type),
0040     iAtt(wgt.iAtt), dAtt(wgt.dAtt), sAtt(wgt.sAtt), auxwt(wgt.auxwt) {}
0041
0042   DirePSWeight( double w, int typeIn, double full, double over, double aux,
0043     int iAttIn=0, double dAttIn=0.0, string sAttIn="") : wt(w), type(typeIn),
0044     iAtt(iAttIn), dAtt(dAttIn), sAtt(sAttIn) { auxwt.push_back(full);
0045     auxwt.push_back(over); auxwt.push_back(aux); }
0046
0047   // Return functions.
0048   double weight()     { return wt; }
0049   int    iAttribute() { return iAtt; }
0050   double dAttribute() { return dAtt; }
0051   string sAttribute() { return sAtt; }
0052
0053   // Set functions.
0054   void setWeight  ( double w) { wt  = w; }
0055   inline DirePSWeight& operator*=(double f) { wt *= f; return *this; }
0056
0057 private:
0058
0059   // Value of the weight.
0060   double wt;
0061
0062   // Remember type: 1-> Accept weight, 0-> Not assigned, -1->Reject weight
0063   int type;
0064
0065   // Auxiliary attributes and identifiers.
0066   int iAtt;
0067   double dAtt;
0068   string sAtt;
0069
0070   // Auxiliary weights, to e.g. store kernel, overestimate and auxiliary
0071   // estimate separately.
0072   // Ordering:
0073   // auxwt[0] -> full kernel,
0074   // auxwt[1] -> overestimate,
0075   // auxwt[2] -> auxiliary overestimate
0076   vector<double> auxwt;
0077
0078 };
0079
0080 // Container for all shower weights, including handling.
0081
0082 class DireWeightContainer {
0083
0084 public:
0085
0086   // Constructor.
0087   DireWeightContainer() : card(""), hasMEs(false), settingsPtr(nullptr),
0088     beamA(nullptr), beamB(nullptr), infoPtr(nullptr), direInfoPtr(nullptr)
0089     { init(); }
0090
0091   DireWeightContainer(Settings* settingsPtrIn) : card(""), hasMEs(false),
0092     settingsPtr(settingsPtrIn), beamA(nullptr), beamB(nullptr),
0093     infoPtr(nullptr), direInfoPtr(nullptr)
0094     { init(); }
0095
0096   // Initialize weights.
0097   void init() {
0098     reset();
0099     for ( unordered_map<string, double>::iterator itw = showerWeight.begin();
0100       itw != showerWeight.end(); ++itw ) itw->second = 1.;
0101   }
0102   void setup();
0103
0104   void initPtrs(BeamParticle* beamAIn, BeamParticle* beamBIn,
0105     Settings* settingsPtrIn, Info* infoPtrIn, DireInfo* direInfoPtrIn) {
0106     beamA    = beamAIn;
0107     beamB    = beamBIn;
0108     settingsPtr = settingsPtrIn;
0109     infoPtr = infoPtrIn;
0110     direInfoPtr = direInfoPtrIn;
0111     return;
0112   }
0113
0114   // Reset current accept/reject probabilities.
0115   void reset() {
0116     for ( unordered_map<string, map<ulong, DirePSWeight> >::iterator
0117       it = rejectWeight.begin(); it != rejectWeight.end(); ++it )
0118       it->second.clear();
0119     for ( unordered_map<string, map<ulong, DirePSWeight> >::iterator
0120       it = acceptWeight.begin(); it != acceptWeight.end(); ++it )
0121       it->second.clear();
0122   }
0123   void clear() {
0124     reset();
0125     for ( unordered_map<string, double >::iterator it = showerWeight.begin();
0126       it != showerWeight.end(); ++it ) it->second = 1.0;
0127   }
0128
0129   // Function to initialize new maps for a new shower variation.
0130   void bookWeightVar(string varKey, bool checkSettings = true);
0131
0132   // To avoid rounding problems, maps will be indexed with long keys.
0133   // Round double inputs to four decimals, as long will should be >10 digits.
0134   ulong key(double a) { return (ulong)(a*1e8+0.5); }
0135   double dkey(ulong a) { return (double(a)/1e8); }
0136
0137   void setWeight( string varKey, double value);
0138   void resetAcceptWeight( double pT2key, double value, string varKey);
0139   void resetRejectWeight( double pT2key, double value, string varKey);
0140   void eraseAcceptWeight( double pT2key, string varKey);
0141   void eraseRejectWeight( double pT2key, string varKey);
0142   double getAcceptWeight( double pT2key, string varKey);
0143   double getRejectWeight( double pT2key, string varKey);
0144
0145   // Attach accept/reject probabilities for a proposed shower step.
0146   void insertWeights( map<double,double> aWeight,
0147     multimap<double,double> rWeight, string varKey );
0148
0149   // Function to calculate the weight of the shower evolution step.
0150   void calcWeight(double pT2, bool includeAcceptAtPT2 = true,
0151     bool includeRejectAtPT2 = false);
0152
0153   pair<double,double> getWeight(double pT2, string valKey = "base");
0154
0155   // Function to return weight of the shower evolution.
0156   double getShowerWeight(string valKey = "base") {
0157     // First try to return an individual shower weight indexed by "valKey".
0158     unordered_map<string, double>::iterator it1 = showerWeight.find( valKey );
0159     if ( it1 != showerWeight.end() ) return it1->second;
0160
0161     // If not possible, return a product of shower weights indexed by "valKey".
0162     unordered_map<string, vector<string> >::iterator it2
0163       = weightCombineList.find(valKey);
0164     if ( it2 != weightCombineList.end() ) {
0165       double wtNow = 1.;
0166       // Loop through group of weights and combine all weights into one weight.
0167       for (int iwgtname=0; iwgtname < int(it2->second.size()); ++iwgtname) {
0168         unordered_map<string, double>::iterator it3
0169           = showerWeight.find( it2->second[iwgtname] );
0170        if ( it3 != showerWeight.end() ) wtNow *= it3->second;
0171       }
0172       return wtNow;
0173     }
0174     // Done.
0175     return 0.;
0176   }
0177
0178   unordered_map<string,double>* getShowerWeights() { return &showerWeight; }
0179   double sizeWeights() const { return showerWeight.size(); }
0180   string weightName (int i) const  { return weightNames[i]; }
0181
0182   double sizeWeightgroups() const { return weightCombineList.size(); }
0183   string weightgroupName (int i) const  {
0184     return weightCombineListNames[i];
0185   }
0186
0187   // Returns additional user-supplied enhancements factors.
0188   double enhanceOverestimate( string name );
0189   double getTrialEnhancement(double pT2key);
0190   void   clearTrialEnhancements() { trialEnhancements.clear(); }
0191   void   addTrialEnhancement( double pT2key, double value) {
0192     trialEnhancements.insert(make_pair(key(pT2key), value));
0193   }
0194
0195   // MG5 matrix element access.
0196   string card;
0197   ExternalMEsPtr matrixElements{};
0198   bool hasMEs;
0199
0200   bool hasME(vector <int> in_pdgs = vector<int>(),
0201     vector<int> out_pdgs = vector<int>() );
0202   bool hasME(const Event& event);
0203   double getME(const Event& event);
0204
0205 private:
0206
0207   static const double LARGEWT;
0208
0209   Settings* settingsPtr;
0210
0211   unordered_map<string, map<ulong, DirePSWeight> > acceptWeight;
0212   unordered_map<string, map<ulong, DirePSWeight> > rejectWeight;
0213   unordered_map<string, double> showerWeight;
0214   vector<string> weightNames;
0215   unordered_map<string, vector<string> > weightCombineList;
0216   vector<string> weightCombineListNames;
0217
0218   // Additonal enhancement factors to boost emission probabilities.
0219   unordered_map<string,double> enhanceFactors;
0220   map<ulong, double> trialEnhancements;
0221
0222   BeamParticle* beamA;
0223   BeamParticle* beamB;
0224   Info* infoPtr;
0225   DireInfo* direInfoPtr;
0226
0227 };
0228
0229 //==========================================================================
0230
0231 } // end namespace Pythia8
0232
0233 #endif // Pythia8_DireWeightContainer_H