StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireBasics.h
1 // DireBasics.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 Dire basics.
7 
8 #ifndef Pythia8_DireBasics_H
9 #define Pythia8_DireBasics_H
10 
11 #define DIRE_BASICS_VERSION "2.002"
12 
13 #include "Pythia8/Event.h"
14 #include <limits>
15 #include <unordered_map>
16 
17 using std::unordered_map;
18 
19 namespace Pythia8 {
20 
21 bool checkSIJ(const Event& e, double minSIJ=0.0);
22 void printSI(const Event& e);
23 void printSIJ(const Event& e);
24 
25 typedef unsigned long ulong;
26 
27 //==========================================================================
28 
29 // Function to hash string into long integer.
30 
31 ulong shash(const string& str);
32 
33 //==========================================================================
34 
35 // Template to make initializing maps simpler, while not relying on C++11.
36 // Usage: map(createmap<T,U>(a,b)(c,d)(e,f));
37 
38 template <typename T, typename U> class createmap {
39 
40 private:
41 
42  map<T, U> m_map;
43 
44 public:
45 
46  createmap(const T& key, const U& val) { m_map[key] = val; }
47  createmap<T, U>& operator()(const T& key, const U& val) {
48  m_map[key] = val;
49  return *this;
50  }
51  operator map<T, U>() { return m_map; }
52 
53 };
54 
55 //==========================================================================
56 
57 // Template to make initializing maps simpler, while not relying on C++11.
58 // Usage: map(createmap<T,U>(a,b)(c,d)(e,f));
59 
60 template <typename T, typename U> class create_unordered_map {
61 
62 private:
63 
64  unordered_map<T, U> um_map;
65 
66 public:
67 
68  create_unordered_map(const T& key, const U& val) { um_map[key] = val; }
69  create_unordered_map<T, U>& operator()(const T& key, const U& val) {
70  um_map[key] = val;
71  return *this;
72  }
73  operator unordered_map<T, U>() { return um_map; }
74 
75 };
76 
77 //==========================================================================
78 
79 // Template to make initializing maps simpler, while not relying on C++11.
80 // Usage: map(createmap<T,U>(a,b)(c,d)(e,f));
81 
82 template <typename T> class createvector {
83 
84 private:
85 
86  vector<T> m_vector;
87 
88 public:
89 
90  createvector(const T& val) { m_vector.push_back(val); }
91  createvector<T>& operator()(const T& val) {
92  m_vector.push_back(val);
93  return *this;
94  }
95  operator vector<T>() { return m_vector; }
96 
97 };
98 
99 //==========================================================================
100 
101 // Helper function to calculate dilogarithm.
102 
103 double polev(double x,double* coef,int N );
104 // Function to calculate dilogarithm.
105 double dilog(double x);
106 
107 //==========================================================================
108 
109 // Kallen function and derived quantities.
110 
111 double lABC(double a, double b, double c);
112 double bABC(double a, double b, double c);
113 double gABC(double a, double b, double c);
114 
115 //==========================================================================
116 
118 
119 public:
120 
121  virtual double f(double) { return 0.; }
122  virtual double f(double, double) { return 0.; }
123  virtual double f(double, vector<double>) { return 0.; }
124 
125 };
126 
128 
129 public:
130 
131  DireRootFinder() {};
132  virtual ~DireRootFinder() {};
133 
134  double findRootSecant1D( DireFunction* f, double xmin, double xmax,
135  double constant, vector<double> xb = vector<double>(), int N=10 ) {
136  vector<double> x;
137  x.push_back(xmin);
138  x.push_back(xmax);
139  for ( int i=2; i < N; ++i ) {
140  double xn = x[i-1]
141  - ( f->f(x[i-1],xb) - constant)
142  * ( x[i-1] - x[i-2] )
143  / ( f->f(x[i-1],xb) - f->f(x[i-2],xb) );
144  x.push_back(xn);
145  }
146  return x.back();
147  }
148 
149  double findRoot1D( DireFunction* f, double xmin, double xmax,
150  double constant, vector<double> xx = vector<double>(), int N=10,
151  double tol = 1e-10 ) {
152 
153  double a(xmin), b(xmax), c(xmax), d(0.), e(0.),
154  fa(f->f(a,xx)-constant), fb(f->f(b,xx)-constant), fc(fb),
155  p(0.), q(0.), r(0.), s(0.),
156  tol1(tol), xm(0.);
157  double EPS = std::numeric_limits<double>::epsilon();
158 
159  // No root.
160  if ( (fa>0. && fb>0.) || (fa<0. && fb<0.) ) {
161  cout << "no root " << constant << " " << f->f(a,xx) << " " << f->f(b,xx)
162  << endl;
163  return std::numeric_limits<double>::quiet_NaN();
164  }
165 
166  for ( int i=0; i < N; ++i ) {
167 
168  if ( (fb>0. && fc>0.) || (fb<0. && fc<0.) ) {
169  c = a;
170  fc = fa;
171  e = d = b-a;
172  }
173 
174  if ( abs(fc) < abs(fb) ) {
175  a = b;
176  b = c;
177  c = a;
178  fa = fb;
179  fb = fc;
180  fc = fa;
181  }
182 
183  tol1 = 2.*EPS*abs(b) + 0.5*tol;
184  xm = 0.5*(c-b);
185 
186  if (abs(xm) <= tol1 || fb == 0.) return b;
187 
188  if (abs(e) >= tol1 && abs(fa) > abs(fb) ) {
189  s = fb/fa;
190  if ( a == c ) {
191  p = 2.*xm*s;
192  q = 1.-s;
193  } else {
194  q = fa/fc;
195  r = fb/fc;
196  p = s*(2.*xm*q*(q-r) - (b-a)*(r-1.));
197  q = (q-1.)*(r-1.)*(s-1.);
198  }
199  if (p>0.) q = -q;
200  p = abs(p);
201  double min1 = 3.*xm*q - abs(tol1*q);
202  double min2 = abs(e*q);
203  if (2.*p < ((min1 < min2) ? min1 : min2)) {
204  e = d;
205  d = p/q;
206  } else {
207  d = xm;
208  e = d;
209  }
210 
211  } else {
212  d = xm;
213  e = d;
214  }
215 
216  a = b;
217  fa = fb;
218 
219  if (abs(d) > tol1) { b += d; }
220  else {
221  b += (xm> 0.) ? tol1 : -tol1;
222  }
223  fb = f->f(b,xx)-constant;
224  }
225 
226  // Failed. Return NaN
227  return std::numeric_limits<double>::quiet_NaN();
228 
229  }
230 
231 };
232 
233 //==========================================================================
234 
236 
237  public:
238 
239  DireEventInfo() {}
240 
241  // Bookkeeping of soft particles.
242  int sizeSoftPos () const { return softPosSave.size(); }
243  int getSoftPos(unsigned int i) const {
244  return (i > softPosSave.size()-1) ? -1 : softPosSave[i]; }
245  bool isSoft(int iPos) {
246  vector<int>::iterator it = find( softPosSave.begin(),
247  softPosSave.end(), iPos);
248  return (it != softPosSave.end());
249  }
250  void addSoftPos(int iPos) { if (!isSoft(iPos)) softPosSave.push_back(iPos); }
251  void removeSoftPos(int iPos) {
252  vector<int>::iterator it = find( softPosSave.begin(),
253  softPosSave.end(), iPos);
254  if (it != softPosSave.end()) softPosSave.erase(it);
255  }
256  void updateSoftPos(int iPosOld, int iPosNew) {
257  if (isSoft(iPosOld)) removeSoftPos(iPosOld);
258  if (isSoft(iPosNew)) removeSoftPos(iPosNew);
259  addSoftPos(iPosNew);
260  }
261  void updateSoftPosIfMatch(int iPosOld, int iPosNew) {
262  if (isSoft(iPosOld)) {
263  vector<int>::iterator it
264  = find (softPosSave.begin(), softPosSave.end(), iPosOld);
265  *it = iPosNew;
266  }
267  }
268  vector<int> softPos () const { return softPosSave; }
269  void clearSoftPos () { softPosSave.clear(); }
270  void listSoft() const {
271  cout << " 'Soft' particles: ";
272  for (int i=0; i < sizeSoftPos(); ++i) cout << setw(5) << getSoftPos(i);
273  cout << endl;
274  }
275 
276  // Bookkeeping of resonances.
277  void removeResPos(int iPos) {
278  vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPos);
279  if (it == iPosRes.end()) return;
280  iPosRes.erase(it);
281  sort (iPosRes.begin(), iPosRes.end());
282  }
283  void addResPos(int iPos) {
284  vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPos);
285  if (it != iPosRes.end()) return;
286  iPosRes.push_back(iPos);
287  sort (iPosRes.begin(), iPosRes.end());
288  }
289  void updateResPos(int iPosOld, int iPosNew) {
290  vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPosOld);
291  if (it == iPosRes.end()) iPosRes.push_back(iPosNew);
292  else *it = iPosNew;
293  sort (iPosRes.begin(), iPosRes.end());
294  }
295  void updateResPosIfMatch(int iPosOld, int iPosNew) {
296  vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPosOld);
297  if (it != iPosRes.end()) {
298  iPosRes.erase(it);
299  iPosRes.push_back(iPosNew);
300  sort (iPosRes.begin(), iPosRes.end());
301  }
302  }
303  bool isRes(int iPos) {
304  vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPos);
305  return (it != iPosRes.end());
306  }
307  int sizeResPos() const { return iPosRes.size(); }
308  int getResPos(unsigned int i) const {
309  return (i > iPosRes.size()-1) ? -1 : iPosRes[i]; }
310  void clearResPos() { iPosRes.resize(0); }
311  void listRes() const {
312  cout << " 'Resonant' particles: ";
313  for (int i=0; i < sizeResPos(); ++i) cout << setw(5) << getResPos(i);
314  cout << endl;
315  }
316 
317  // Data members.
318  vector<int> softPosSave;
319  vector<int> iPosRes;
320 
321 };
322 
323 
324 //==========================================================================
325 
327 
328  public:
329 
330  DireDebugInfo() = default;
331 
332  void clearMessages() {
333  messageStream0.str("");
334  messageStream1.str("");
335  messageStream2.str("");
336  }
337 
338  void printMessages( int verbosity = 0) {
339  cout << "\n"
340  << "*------------------------------------------------------------*\n"
341  << "*----------------- Begin diagnostic output ------------------*\n\n";
342  if (verbosity == 0) cout << scientific << setprecision(8)
343  << messageStream0.str();
344  if (verbosity == 1) cout << scientific << setprecision(8)
345  << messageStream1.str();
346  if (verbosity == 2) cout << scientific << setprecision(8)
347  << messageStream2.str();
348  cout << "\n\n"
349  << "*----------------- End diagnostic output -------------------*\n"
350  << "*-----------------------------------------------------------*"
351  << endl;
352  }
353 
354  void printHistograms() {}
355 
356  void fillHistograms(int type, int nfinal, double mec, double pT, double z) {
357  if (false) cout << type*nfinal*mec*pT*z;}
358 
359  // Add debug messages to message stream.
360  ostream & message ( int verbosity = 0) {
361  if (verbosity == 0) return messageStream0;
362  if (verbosity == 1) return messageStream1;
363  if (verbosity == 2) return messageStream2;
364  return messageStream0;
365  }
366 
367  // Debug message streams.
368  ostringstream messageStream0, messageStream1, messageStream2;
369 
370 };
371 
372 //==========================================================================
373 
374 class DireInfo {
375 
376  public:
377 
378  DireInfo() {}
379 
380  void clearAll() {
381  direEventInfo.clearResPos();
382  direEventInfo.clearSoftPos();
383  direDebugInfo.clearMessages();
384  }
385 
386  // Resonance info forwards.
387  void removeResPos(int iPos) { return direEventInfo.removeResPos(iPos); }
388  void addResPos(int iPos) { return direEventInfo.addResPos(iPos); }
389  bool isRes(int iPos) { return direEventInfo.isRes(iPos); }
390  void clearResPos () { return direEventInfo.clearResPos(); }
391  int sizeResPos() const { return direEventInfo.sizeResPos(); }
392  void listRes() const { return direEventInfo.listRes(); }
393  int getResPos(unsigned int i) const { return direEventInfo.getResPos(i); }
394  void updateResPos(int iPosOld, int iPosNew) {
395  return direEventInfo.updateResPos(iPosOld,iPosNew); }
396  void updateResPosIfMatch(int iPosOld, int iPosNew) {
397  return direEventInfo.updateResPosIfMatch(iPosOld,iPosNew); }
398 
399  // Debug info forwards.
400  void printMessages( int verbosity = 0) {
401  return direDebugInfo.printMessages(verbosity); }
402  ostream & message ( int verbosity = 0) {
403  return direDebugInfo.message(verbosity); }
404  void clearMessages() { direDebugInfo.clearMessages(); }
405 
406  void fillHistograms(int type, int nfinal, double mec, double pT, double z) {
407  direDebugInfo.fillHistograms(type, nfinal, mec, pT, z); }
408  void printHistograms() { direDebugInfo.printHistograms(); }
409 
410  // Soft particle info forwards.
411  bool isSoft(int iPos) { return direEventInfo.isSoft(iPos); }
412  void addSoftPos(int iPos) { return direEventInfo.addSoftPos(iPos); }
413  void removeSoftPos(int iPos) { return direEventInfo.removeSoftPos(iPos); }
414  vector<int> softPos () { return direEventInfo.softPos(); }
415  void clearSoftPos () { return direEventInfo.clearSoftPos(); }
416  int sizeSoftPos () const { return direEventInfo.sizeSoftPos(); }
417  void listSoft() const { return direEventInfo.listSoft(); }
418  int getSoftPos(unsigned int i) const { return direEventInfo.getSoftPos(i); }
419  void updateSoftPos(int iPosOld, int iPosNew) {
420  return direEventInfo.updateSoftPos(iPosOld, iPosNew);
421  }
422  void updateSoftPosIfMatch(int iPosOld, int iPosNew) {
423  return direEventInfo.updateSoftPosIfMatch(iPosOld, iPosNew);
424  }
425 
426  DireEventInfo direEventInfo;
427  DireDebugInfo direDebugInfo;
428 
429 };
430 
431 //==========================================================================
432 
433 } // end namespace Pythia8
434 
435 #endif // end Pythia8_DireBasics_H
Definition: AgUStep.h:26