StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireWeightContainer.h
1 // DireWeightContainer.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Stefan Prestel, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Header file for weight constainers for Dire parton shower.
7 
8 #ifndef Pythia8_DireWeightContainer_H
9 #define Pythia8_DireWeightContainer_H
10 
11 #define DIRE_WEIGHTCONTAINER_VERSION "2.002"
12 
13 #include "Pythia8/PythiaStdlib.h"
14 #include "Pythia8/Settings.h"
15 #include "Pythia8/DireBasics.h"
16 #include "Pythia8/ShowerMEs.h"
17 
18 namespace Pythia8 {
19 
20 //==========================================================================
21 
22 // Container for a single weight with auxiliary information.
23 
24 class DirePSWeight {
25 
26 public:
27 
28  // Constructors.
29  DirePSWeight()
30  : wt(1.0), type(0), iAtt(0), dAtt(0.0), sAtt("") {}
31 
32  DirePSWeight( double w)
33  : wt(w), type(0), iAtt(0), dAtt(0.0), sAtt("") {}
34 
35  DirePSWeight( double w, int typeIn, int iAttIn=0, double dAttIn=0.0,
36  string sAttIn="") : wt(w), type(typeIn), iAtt(iAttIn), dAtt(dAttIn),
37  sAtt(sAttIn) {}
38 
39  DirePSWeight( const DirePSWeight& wgt) : wt(wgt.wt), type(wgt.type),
40  iAtt(wgt.iAtt), dAtt(wgt.dAtt), sAtt(wgt.sAtt), auxwt(wgt.auxwt) {}
41 
42  DirePSWeight( double w, int typeIn, double full, double over, double aux,
43  int iAttIn=0, double dAttIn=0.0, string sAttIn="") : wt(w), type(typeIn),
44  iAtt(iAttIn), dAtt(dAttIn), sAtt(sAttIn) { auxwt.push_back(full);
45  auxwt.push_back(over); auxwt.push_back(aux); }
46 
47  // Return functions.
48  double weight() { return wt; }
49  int iAttribute() { return iAtt; }
50  double dAttribute() { return dAtt; }
51  string sAttribute() { return sAtt; }
52 
53  // Set functions.
54  void setWeight ( double w) { wt = w; }
55  inline DirePSWeight& operator*=(double f) { wt *= f; return *this; }
56 
57 private:
58 
59  // Value of the weight.
60  double wt;
61 
62  // Remember type: 1-> Accept weight, 0-> Not assigned, -1->Reject weight
63  int type;
64 
65  // Auxiliary attributes and identifiers.
66  int iAtt;
67  double dAtt;
68  string sAtt;
69 
70  // Auxiliary weights, to e.g. store kernel, overestimate and auxiliary
71  // estimate separately.
72  // Ordering:
73  // auxwt[0] -> full kernel,
74  // auxwt[1] -> overestimate,
75  // auxwt[2] -> auxiliary overestimate
76  vector<double> auxwt;
77 
78 };
79 
80 // Container for all shower weights, including handling.
81 
83 
84 public:
85 
86  // Constructor.
87  DireWeightContainer() : card(""), hasMEs(false), settingsPtr(nullptr),
88  beamA(nullptr), beamB(nullptr), infoPtr(nullptr), direInfoPtr(nullptr)
89  { init(); }
90 
91  DireWeightContainer(Settings* settingsPtrIn) : card(""), hasMEs(false),
92  settingsPtr(settingsPtrIn), beamA(nullptr), beamB(nullptr),
93  infoPtr(nullptr), direInfoPtr(nullptr)
94  { init(); }
95 
96  // Initialize weights.
97  void init() {
98  reset();
99  for ( unordered_map<string, double>::iterator itw = showerWeight.begin();
100  itw != showerWeight.end(); ++itw ) itw->second = 1.;
101  }
102  void setup();
103 
104  void initPtrs(BeamParticle* beamAIn, BeamParticle* beamBIn,
105  Settings* settingsPtrIn, Info* infoPtrIn, DireInfo* direInfoPtrIn) {
106  beamA = beamAIn;
107  beamB = beamBIn;
108  settingsPtr = settingsPtrIn;
109  infoPtr = infoPtrIn;
110  direInfoPtr = direInfoPtrIn;
111  return;
112  }
113 
114  // Reset current accept/reject probabilities.
115  void reset() {
116  for ( unordered_map<string, map<ulong, DirePSWeight> >::iterator
117  it = rejectWeight.begin(); it != rejectWeight.end(); ++it )
118  it->second.clear();
119  for ( unordered_map<string, map<ulong, DirePSWeight> >::iterator
120  it = acceptWeight.begin(); it != acceptWeight.end(); ++it )
121  it->second.clear();
122  }
123  void clear() {
124  reset();
125  for ( unordered_map<string, double >::iterator it = showerWeight.begin();
126  it != showerWeight.end(); ++it ) it->second = 1.0;
127  }
128 
129  // Function to initialize new maps for a new shower variation.
130  void bookWeightVar(string varKey, bool checkSettings = true);
131 
132  // To avoid rounding problems, maps will be indexed with long keys.
133  // Round double inputs to four decimals, as long will should be >10 digits.
134  ulong key(double a) { return (ulong)(a*1e8+0.5); }
135  double dkey(ulong a) { return (double(a)/1e8); }
136 
137  void setWeight( string varKey, double value);
138  void resetAcceptWeight( double pT2key, double value, string varKey);
139  void resetRejectWeight( double pT2key, double value, string varKey);
140  void eraseAcceptWeight( double pT2key, string varKey);
141  void eraseRejectWeight( double pT2key, string varKey);
142  double getAcceptWeight( double pT2key, string varKey);
143  double getRejectWeight( double pT2key, string varKey);
144 
145  // Attach accept/reject probabilities for a proposed shower step.
146  void insertWeights( map<double,double> aWeight,
147  multimap<double,double> rWeight, string varKey );
148 
149  // Function to calculate the weight of the shower evolution step.
150  void calcWeight(double pT2, bool includeAcceptAtPT2 = true,
151  bool includeRejectAtPT2 = false);
152 
153  pair<double,double> getWeight(double pT2, string valKey = "base");
154 
155  // Function to return weight of the shower evolution.
156  double getShowerWeight(string valKey = "base") {
157  // First try to return an individual shower weight indexed by "valKey".
158  unordered_map<string, double>::iterator it1 = showerWeight.find( valKey );
159  if ( it1 != showerWeight.end() ) return it1->second;
160 
161  // If not possible, return a product of shower weights indexed by "valKey".
162  unordered_map<string, vector<string> >::iterator it2
163  = weightCombineList.find(valKey);
164  if ( it2 != weightCombineList.end() ) {
165  double wtNow = 1.;
166  // Loop through group of weights and combine all weights into one weight.
167  for (int iwgtname=0; iwgtname < int(it2->second.size()); ++iwgtname) {
168  unordered_map<string, double>::iterator it3
169  = showerWeight.find( it2->second[iwgtname] );
170  if ( it3 != showerWeight.end() ) wtNow *= it3->second;
171  }
172  return wtNow;
173  }
174  // Done.
175  return 0.;
176  }
177 
178  unordered_map<string,double>* getShowerWeights() { return &showerWeight; }
179  double sizeWeights() const { return showerWeight.size(); }
180  string weightName (int i) const { return weightNames[i]; }
181 
182  double sizeWeightgroups() const { return weightCombineList.size(); }
183  string weightgroupName (int i) const {
184  return weightCombineListNames[i];
185  }
186 
187  // Returns additional user-supplied enhancements factors.
188  double enhanceOverestimate( string name );
189  double getTrialEnhancement(double pT2key);
190  void clearTrialEnhancements() { trialEnhancements.clear(); }
191  void addTrialEnhancement( double pT2key, double value) {
192  trialEnhancements.insert(make_pair(key(pT2key), value));
193  }
194 
195  // MG5 matrix element access.
196  string card;
197  ShowerMEsPlugin matrixElements;
198  bool hasMEs;
199 
200  bool hasME(vector <int> in_pdgs = vector<int>(),
201  vector<int> out_pdgs = vector<int>() );
202  bool hasME(const Event& event);
203  double getME(const Event& event);
204 
205 private:
206 
207  static const double LARGEWT;
208 
209  Settings* settingsPtr;
210 
211  unordered_map<string, map<ulong, DirePSWeight> > acceptWeight;
212  unordered_map<string, map<ulong, DirePSWeight> > rejectWeight;
213  unordered_map<string, double> showerWeight;
214  vector<string> weightNames;
215  unordered_map<string, vector<string> > weightCombineList;
216  vector<string> weightCombineListNames;
217 
218  // Additonal enhancement factors to boost emission probabilities.
219  unordered_map<string,double> enhanceFactors;
220  map<ulong, double> trialEnhancements;
221 
222  BeamParticle* beamA;
223  BeamParticle* beamB;
224  Info* infoPtr;
225  DireInfo* direInfoPtr;
226 
227 };
228 
229 //==========================================================================
230 
231 } // end namespace Pythia8
232 
233 #endif // end Pythia8_DireWeightContainer_H