StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Info.h
1 // Info.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file contains a class that keep track of generic event info.
7 // Info: contains information on the generation process and errors.
8 
9 #ifndef Pythia8_Info_H
10 #define Pythia8_Info_H
11 
12 #include "Pythia8/PythiaStdlib.h"
13 
14 namespace Pythia8 {
15 
16 //==========================================================================
17 
18 // The Info class contains a mixed bag of information on the event
19 // generation activity, especially on the current subprocess properties,
20 // and on the number of errors encountered. This is used by the
21 // generation machinery, but can also be read by the user.
22 // Note: some methods that maybe should not be accessible to the user
23 // are still public, to work also for user-written FSR/ISR classes.
24 
25 class Info {
26 
27 public:
28 
29  // Constructor.
30  Info() : eCMSave(0.), lowPTmin(false), a0MPISave(0.), weightCKKWLSave(1.),
31  weightFIRSTSave(0.) {
32  for (int i = 0; i < 40; ++i) counters[i] = 0;}
33 
34  // Listing of most available information on current event.
35  void list(ostream& os = cout) const;
36 
37  // Beam particles (in rest frame). CM energy of event.
38  int idA() const {return idASave;}
39  int idB() const {return idBSave;}
40  double pzA() const {return pzASave;}
41  double pzB() const {return pzBSave;}
42  double eA() const {return eASave;}
43  double eB() const {return eBSave;}
44  double mA() const {return mASave;}
45  double mB() const {return mBSave;}
46  double eCM() const {return eCMSave;}
47  double s() const {return sSave;}
48 
49  // Warnings from initialization.
50  bool tooLowPTmin() const {return lowPTmin;}
51 
52  // Process name and code, and the number of final-state particles.
53  string name() const {return nameSave;}
54  int code() const {return codeSave;}
55  int nFinal() const {return nFinalSave;}
56 
57  // Are beam particles resolved, with pdf's? Are they diffractive?
58  bool isResolved() const {return isRes;}
59  bool isDiffractiveA() const {return isDiffA;}
60  bool isDiffractiveB() const {return isDiffB;}
61  bool isDiffractiveC() const {return isDiffC;}
62  bool isNonDiffractive() const {return isND;}
63  // Retained for backwards compatibility.
64  bool isMinBias() const {return isND;}
65 
66  // Information for Les Houches Accord and reading files.
67  bool isLHA() const {return isLH;}
68  bool atEndOfFile() const {return atEOF;}
69 
70  // For nondiffractive and Les Houches Accord identify hardest subprocess.
71  bool hasSub(int i = 0) const {return hasSubSave[i];}
72  string nameSub(int i = 0) const {return nameSubSave[i];}
73  int codeSub(int i = 0) const {return codeSubSave[i];}
74  int nFinalSub(int i = 0) const {return nFinalSubSave[i];}
75 
76  // Incoming parton flavours and x values.
77  int id1(int i = 0) const {return id1Save[i];}
78  int id2(int i = 0) const {return id2Save[i];}
79  double x1(int i = 0) const {return x1Save[i];}
80  double x2(int i = 0) const {return x2Save[i];}
81  double y(int i = 0) const {return 0.5 * log( x1Save[i] / x2Save[i]);}
82  double tau(int i = 0) const {return x1Save[i] * x2Save[i];}
83 
84  // Hard process flavours, x values, parton densities, couplings, Q2 scales.
85  int id1pdf(int i = 0) const {return id1pdfSave[i];}
86  int id2pdf(int i = 0) const {return id2pdfSave[i];}
87  double x1pdf(int i = 0) const {return x1pdfSave[i];}
88  double x2pdf(int i = 0) const {return x2pdfSave[i];}
89  double pdf1(int i = 0) const {return pdf1Save[i];}
90  double pdf2(int i = 0) const {return pdf2Save[i];}
91  double QFac(int i = 0) const {return sqrtpos(Q2FacSave[i]);}
92  double Q2Fac(int i = 0) const {return Q2FacSave[i];}
93  bool isValence1() const {return isVal1;}
94  bool isValence2() const {return isVal2;}
95  double alphaS(int i = 0) const {return alphaSSave[i];}
96  double alphaEM(int i = 0) const {return alphaEMSave[i];}
97  double QRen(int i = 0) const {return sqrtpos(Q2RenSave[i]);}
98  double Q2Ren(int i = 0) const {return Q2RenSave[i];}
99  double scalup(int i = 0) const {return scalupSave[i];}
100 
101  // Mandelstam variables (notation as if subcollision).
102  double mHat(int i = 0) const {return sqrt(sH[i]);}
103  double sHat(int i = 0) const {return sH[i];}
104  double tHat(int i = 0) const {return tH[i];}
105  double uHat(int i = 0) const {return uH[i];}
106  double pTHat(int i = 0) const {return pTH[i];}
107  double pT2Hat(int i = 0) const {return pTH[i] * pTH[i];}
108  double m3Hat(int i = 0) const {return m3H[i];}
109  double m4Hat(int i = 0) const {return m4H[i];}
110  double thetaHat(int i = 0) const {return thetaH[i];}
111  double phiHat(int i = 0) const {return phiH[i];}
112 
113  // Weight of current event; normally 1, but used for Les Houches events
114  // or when reweighting phase space selection. Conversion from mb to pb
115  // for LHA strategy +-4. Also cumulative sum.
116  double weight() const;
117  double weightSum() const;
118  double lhaStrategy() const {return lhaStrategySave;}
119 
120  // Number of times other steps have been carried out.
121  int nISR() const {return nISRSave;}
122  int nFSRinProc() const {return nFSRinProcSave;}
123  int nFSRinRes() const {return nFSRinResSave;}
124 
125  // Maximum pT scales for MPI, ISR and FSR (in hard process).
126  double pTmaxMPI() const {return pTmaxMPISave;}
127  double pTmaxISR() const {return pTmaxISRSave;}
128  double pTmaxFSR() const {return pTmaxFSRSave;}
129 
130  // Current evolution scale (for UserHooks).
131  double pTnow() const {return pTnowSave;}
132 
133  // Impact parameter picture, global information
134  double a0MPI() const {return a0MPISave;}
135 
136  // Impact parameter picture, as set by hardest interaction.
137  double bMPI() const {return (bIsSet) ? bMPISave : 1.;}
138  double enhanceMPI() const {return (bIsSet) ? enhanceMPISave : 1.;}
139  double eMPI(int i) const {return (bIsSet) ? eMPISave[i] : 1.;}
140 
141  // Number of multiparton interactions, with code and pT for them.
142  int nMPI() const {return nMPISave;}
143  int codeMPI(int i) const {return codeMPISave[i];}
144  double pTMPI(int i) const {return pTMPISave[i];}
145  int iAMPI(int i) const {return iAMPISave[i];}
146  int iBMPI(int i) const {return iBMPISave[i];}
147 
148  // Cross section estimate, optionally process by process.
149  vector<int> codesHard();
150  string nameProc(int i = 0) {return (i == 0) ? "sum"
151  : ( (procNameM[i] == "") ? "unknown process" : procNameM[i] );}
152  long nTried(int i = 0) {return (i == 0) ? nTry : nTryM[i];}
153  long nSelected(int i = 0) {return (i == 0) ? nSel : nSelM[i];}
154  long nAccepted(int i = 0) {return (i == 0) ? nAcc : nAccM[i];}
155  double sigmaGen(int i = 0) {return (i == 0) ? sigGen : sigGenM[i];}
156  double sigmaErr(int i = 0) {return (i == 0) ? sigErr : sigErrM[i];}
157 
158  // Counters for number of loops in various places.
159  int getCounter( int i) const {return counters[i];}
160 
161  // Set or increase the value stored in a counter.
162  void setCounter( int i, int value = 0) {counters[i] = value;}
163  void addCounter( int i, int value = 1) {counters[i] += value;}
164 
165  // Reset to empty map of error messages.
166  void errorReset() {messages.clear();}
167 
168  // Print a message the first few times. Insert in database.
169  void errorMsg(string messageIn, string extraIn = " ",
170  bool showAlways = false, ostream& os = cout);
171 
172  // Provide total number of errors/aborts/warnings experienced to date.
173  int errorTotalNumber();
174 
175  // Print statistics on errors/aborts/warnings.
176  void errorStatistics(ostream& os = cout);
177 
178  // Set initialization warning flag when too low pTmin in ISR/FSR/MPI.
179  void setTooLowPTmin(bool lowPTminIn) {lowPTmin = lowPTminIn;}
180 
181  // Set info on valence character of hard collision partons.
182  void setValence( bool isVal1In, bool isVal2In) {isVal1 = isVal1In;
183  isVal2 = isVal2In;}
184 
185  // Set and get some MPI/ISR/FSR properties needed for matching,
186  // i.e. mainly of internal relevance.
187  void hasHistory( bool hasHistoryIn) {hasHistorySave = hasHistoryIn;}
188  bool hasHistory() {return hasHistorySave;}
189  void zNowISR( double zNowIn) {zNowISRSave = zNowIn;}
190  double zNowISR() {return zNowISRSave;}
191  void pT2NowISR( double pT2NowIn) {pT2NowISRSave = pT2NowIn;}
192  double pT2NowISR() {return pT2NowISRSave;}
193 
194  // Return CKKW-L weight.
195  double getWeightCKKWL() const { return weightCKKWLSave;}
196  // Set CKKW-L weight.
197  void setWeightCKKWL( double weightIn) { weightCKKWLSave = weightIn;}
198  // Return merging weight.
199  double mergingWeight() const { return weightCKKWLSave;}
200 
201  // Return the complete NLO weight.
202  double mergingWeightNLO() const {
203  return (weightCKKWLSave - weightFIRSTSave); }
204  // Return the O(\alpha_s)-term of the CKKW-L weight.
205  double getWeightFIRST() const { return weightFIRSTSave;}
206  // Set the O(\alpha_s)-term of the CKKW-L weight.
207  void setWeightFIRST( double weightIn) { weightFIRSTSave = weightIn;}
208 
209  // Return an LHEF header
210  string header(const string &key) {
211  if (headers.find(key) == headers.end()) return "";
212  else return headers[key];
213  }
214 
215  // Return a list of all header key names
216  vector < string > headerKeys() {
217  vector < string > keys;
218  for (map < string, string >::iterator it = headers.begin();
219  it != headers.end(); it++)
220  keys.push_back(it->first);
221  return keys;
222  }
223 
224 private:
225 
226  // Number of times the same error message is repeated, unless overridden.
227  static const int TIMESTOPRINT;
228 
229  // Allow conversion from mb to pb.
230  static const double CONVERTMB2PB;
231 
232  // Store common beam quantities.
233  int idASave, idBSave;
234  double pzASave, eASave,mASave, pzBSave, eBSave, mBSave, eCMSave, sSave;
235 
236  // Store initialization information.
237  bool lowPTmin;
238 
239  // Store common integrated cross section quantities.
240  long nTry, nSel, nAcc;
241  double sigGen, sigErr, wtAccSum;
242  map<int, string> procNameM;
243  map<int, long> nTryM, nSelM, nAccM;
244  map<int, double> sigGenM, sigErrM;
245  int lhaStrategySave;
246 
247  // Store common MPI information.
248  double a0MPISave;
249 
250  // Store current-event quantities.
251  bool isRes, isDiffA, isDiffB, isDiffC, isND, isLH, hasSubSave[4],
252  bIsSet, evolIsSet, atEOF, isVal1, isVal2, hasHistorySave;
253  int codeSave, codeSubSave[4], nFinalSave, nFinalSubSave[4], nTotal,
254  id1Save[4], id2Save[4], id1pdfSave[4], id2pdfSave[4], nMPISave,
255  nISRSave, nFSRinProcSave, nFSRinResSave;
256  double x1Save[4], x2Save[4], x1pdfSave[4], x2pdfSave[4], pdf1Save[4],
257  pdf2Save[4], Q2FacSave[4], alphaEMSave[4], alphaSSave[4],
258  Q2RenSave[4], scalupSave[4], sH[4], tH[4], uH[4], pTH[4], m3H[4],
259  m4H[4], thetaH[4], phiH[4], weightSave, bMPISave, enhanceMPISave,
260  pTmaxMPISave, pTmaxISRSave, pTmaxFSRSave, pTnowSave,
261  zNowISRSave, pT2NowISRSave;
262  string nameSave, nameSubSave[4];
263  vector<int> codeMPISave, iAMPISave, iBMPISave;
264  vector<double> pTMPISave, eMPISave;
265 
266  // Vector of various loop counters.
267  int counters[50];
268 
269  // Map for all error messages.
270  map<string, int> messages;
271 
272  // Map for LHEF headers
273  map<string, string> headers;
274 
275  // Friend classes allowed to set info.
276  friend class Pythia;
277  friend class ProcessLevel;
278  friend class ProcessContainer;
279  friend class PartonLevel;
280  friend class MultipartonInteractions;
281  friend class LHAup;
282 
283  // Set info on the two incoming beams: only from Pythia class.
284  void setBeamA( int idAin, double pzAin, double eAin, double mAin) {
285  idASave = idAin; pzASave = pzAin; eASave = eAin; mASave = mAin;}
286  void setBeamB( int idBin, double pzBin, double eBin, double mBin) {
287  idBSave = idBin; pzBSave = pzBin; eBSave = eBin; mBSave = mBin;}
288  void setECM( double eCMin) {eCMSave = eCMin; sSave = eCMSave * eCMSave;}
289 
290  // Reset info for current event: only from Pythia class.
291  void clear() {
292  isRes = isDiffA = isDiffB = isDiffC = isND = isLH = bIsSet
293  = evolIsSet = atEOF = isVal1 = isVal2 = hasHistorySave = false;
294  codeSave = nFinalSave = nTotal = nMPISave = nISRSave = nFSRinProcSave
295  = nFSRinResSave = 0;
296  weightSave = bMPISave = enhanceMPISave = weightCKKWLSave = 1.;
297  pTmaxMPISave = pTmaxISRSave = pTmaxFSRSave = pTnowSave = zNowISRSave
298  = pT2NowISRSave = weightFIRSTSave = 0.;
299  nameSave = " ";
300  for (int i = 0; i < 4; ++i) {
301  hasSubSave[i] = false;
302  codeSubSave[i] = nFinalSubSave[i] = id1pdfSave[i] = id2pdfSave[i]
303  = id1Save[i] = id2Save[i] = 0;
304  x1pdfSave[i] = x2pdfSave[i] = pdf1Save[i] = pdf2Save[i]
305  = Q2FacSave[i] = alphaEMSave[i] = alphaSSave[i] = Q2RenSave[i]
306  = scalupSave[i] = x1Save[i] = x2Save[i] = sH[i] = tH[i] = uH[i]
307  = pTH[i] = m3H[i] = m4H[i] = thetaH[i] = phiH[i] = 0.;
308  nameSubSave[i] = " ";
309  }
310  codeMPISave.resize(0); iAMPISave.resize(0); iBMPISave.resize(0);
311  pTMPISave.resize(0); eMPISave.resize(0); }
312 
313  // Reset info arrays only for parton and hadron level.
314  int sizeMPIarrays() const { return codeMPISave.size(); }
315  void resizeMPIarrays(int newSize) { codeMPISave.resize(newSize);
316  iAMPISave.resize(newSize); iBMPISave.resize(newSize);
317  pTMPISave.resize(newSize); eMPISave.resize(newSize); }
318 
319  // Set info on the (sub)process: from ProcessLevel, ProcessContainer or
320  // MultipartonInteractions classes.
321  void setType( string nameIn, int codeIn, int nFinalIn,
322  bool isNonDiffIn = false, bool isResolvedIn = true,
323  bool isDiffractiveAin = false, bool isDiffractiveBin = false,
324  bool isDiffractiveCin = false, bool isLHAin = false) {
325  nameSave = nameIn; codeSave = codeIn; nFinalSave = nFinalIn;
326  isND = isNonDiffIn; isRes = isResolvedIn; isDiffA = isDiffractiveAin;
327  isDiffB = isDiffractiveBin; isDiffC = isDiffractiveCin; isLH = isLHAin;
328  nTotal = 2 + nFinalSave; bIsSet = false; hasSubSave[0] = false;
329  nameSubSave[0] = " "; codeSubSave[0] = 0; nFinalSubSave[0] = 0;
330  evolIsSet = false;}
331  void setSubType( int iDS, string nameSubIn, int codeSubIn,
332  int nFinalSubIn) { hasSubSave[iDS] = true; nameSubSave[iDS] = nameSubIn;
333  codeSubSave[iDS] = codeSubIn; nFinalSubSave[iDS] = nFinalSubIn;}
334  void setPDFalpha( int iDS, int id1pdfIn, int id2pdfIn, double x1pdfIn,
335  double x2pdfIn, double pdf1In, double pdf2In, double Q2FacIn,
336  double alphaEMIn, double alphaSIn, double Q2RenIn, double scalupIn) {
337  id1pdfSave[iDS] = id1pdfIn; id2pdfSave[iDS] = id2pdfIn;
338  x1pdfSave[iDS] = x1pdfIn; x2pdfSave[iDS] = x2pdfIn;
339  pdf1Save[iDS] = pdf1In; pdf2Save[iDS] = pdf2In; Q2FacSave[iDS] = Q2FacIn;
340  alphaEMSave[iDS] = alphaEMIn; alphaSSave[iDS] = alphaSIn;
341  Q2RenSave[iDS] = Q2RenIn; scalupSave[iDS] = scalupIn;}
342  void setScalup( int iDS, double scalupIn) {scalupSave[iDS] = scalupIn;}
343  void setKin( int iDS, int id1In, int id2In, double x1In, double x2In,
344  double sHatIn, double tHatIn, double uHatIn, double pTHatIn,
345  double m3HatIn, double m4HatIn, double thetaHatIn, double phiHatIn) {
346  id1Save[iDS] = id1In; id2Save[iDS] = id2In; x1Save[iDS] = x1In;
347  x2Save[iDS] = x2In; sH[iDS] = sHatIn; tH[iDS] = tHatIn;
348  uH[iDS] = uHatIn; pTH[iDS] = pTHatIn; m3H[iDS] = m3HatIn;
349  m4H[iDS] = m4HatIn; thetaH[iDS] = thetaHatIn; phiH[iDS] = phiHatIn;}
350  void setTypeMPI( int codeMPIIn, double pTMPIIn, int iAMPIIn = 0,
351  int iBMPIIn = 0, double eMPIIn = 1.) {codeMPISave.push_back(codeMPIIn);
352  pTMPISave.push_back(pTMPIIn); iAMPISave.push_back(iAMPIIn);
353  iBMPISave.push_back(iBMPIIn); eMPISave.push_back(eMPIIn); }
354 
355  // Set info on cross section: from ProcessLevel (reset in Pythia).
356  void sigmaReset() { nTry = nSel = nAcc = 0; sigGen = sigErr = wtAccSum = 0.;
357  procNameM.clear(); nTryM.clear(); nSelM.clear(); nAccM.clear();
358  sigGenM.clear(); sigErrM.clear();}
359  void setSigma( int i, string procNameIn, long nTryIn, long nSelIn,
360  long nAccIn, double sigGenIn, double sigErrIn, double wtAccSumIn) {
361  if (i == 0) {nTry = nTryIn; nSel = nSelIn; nAcc = nAccIn;
362  sigGen = sigGenIn; sigErr = sigErrIn; wtAccSum = wtAccSumIn;}
363  else {procNameM[i] = procNameIn; nTryM[i] = nTryIn; nSelM[i] = nSelIn;
364  nAccM[i] = nAccIn; sigGenM[i] = sigGenIn; sigErrM[i] = sigErrIn;} }
365  void addSigma( int i, long nTryIn, long nSelIn, long nAccIn, double sigGenIn,
366  double sigErrIn) { nTryM[i] += nTryIn; nSelM[i] += nSelIn;
367  nAccM[i] += nAccIn; sigGenM[i] += sigGenIn;
368  sigErrM[i] = sqrtpos(sigErrM[i]*sigErrM[i] + sigErrIn*sigErrIn); }
369 
370  // Set info on impact parameter: from PartonLevel.
371  void setImpact( double bMPIIn, double enhanceMPIIn, bool bIsSetIn = true) {
372  bMPISave = bMPIIn; enhanceMPISave = eMPISave[0] = enhanceMPIIn,
373  bIsSet = bIsSetIn;}
374 
375  // Set info on pTmax scales and number of evolution steps: from PartonLevel.
376  void setPartEvolved( int nMPIIn, int nISRIn) {
377  nMPISave = nMPIIn; nISRSave = nISRIn;}
378  void setEvolution( double pTmaxMPIIn, double pTmaxISRIn, double pTmaxFSRIn,
379  int nMPIIn, int nISRIn, int nFSRinProcIn, int nFSRinResIn) {
380  pTmaxMPISave = pTmaxMPIIn; pTmaxISRSave = pTmaxISRIn;
381  pTmaxFSRSave = pTmaxFSRIn; nMPISave = nMPIIn; nISRSave = nISRIn;
382  nFSRinProcSave = nFSRinProcIn; nFSRinResSave = nFSRinResIn;
383  evolIsSet = true;}
384 
385  // Set current pT evolution scale for MPI/ISR/FSR; from PartonLevel.
386  void setPTnow( double pTnowIn) {pTnowSave = pTnowIn;}
387 
388  // Set a0 from MultipartonInteractions.
389  void seta0MPI(double a0MPIin) {a0MPISave = a0MPIin;}
390 
391  // Set info whether reading of Les Houches Accord file at end.
392  void setEndOfFile( bool atEOFin) {atEOF = atEOFin;}
393 
394  // Set event weight; currently only for Les Houches description.
395  void setWeight( double weightIn, int lhaStrategyIn) {
396  weightSave = weightIn; lhaStrategySave = lhaStrategyIn; }
397 
398  // Save merging weight (i.e. CKKW-L-type weight, summed O(\alpha_s) weight)
399  double weightCKKWLSave, weightFIRSTSave;
400 
401  // Set LHEF headers
402  void setHeader(const string &key, const string &val)
403  { headers[key] = val; }
404 
405 };
406 
407 //==========================================================================
408 
409 } // end namespace Pythia8
410 
411 #endif // Pythia8_Info_H