EIC code displayed by LXR master/​include/​Pythia8/​DireBasics.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

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

 // DireBasics.h is a part of the PYTHIA event generator.
 // Copyright (C) 2024 Stefan Prestel, Torbjorn Sjostrand.
 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
 // Please respect the MCnet Guidelines, see GUIDELINES for details.
 
 // Header file for Dire basics.
 
 #ifndef Pythia8_DireBasics_H
 #define Pythia8_DireBasics_H
 
 #define DIRE_BASICS_VERSION "2.002"
 
 #include "Pythia8/Event.h"
 #include <limits>
 #include <unordered_map>
 
 using std::unordered_map;
 
 namespace Pythia8 {
 
 bool checkSIJ(const Event& e, double minSIJ=0.0);
 void printSI(const Event& e);
 void printSIJ(const Event& e);
 
 typedef unsigned long ulong;
 
 //==========================================================================
 
 // Function to hash string into long integer.
 
 ulong shash(const string& str);
 
 //==========================================================================
 
 // Template to make initializing maps simpler, while not relying on C++11.
 // Usage: map(createmap<T,U>(a,b)(c,d)(e,f));
 
 template <typename T, typename U> class createmap {
 
 private:
 
   map<T, U> m_map;
 
 public:
 
   createmap(const T& key, const U& val) { m_map[key] = val; }
   createmap<T, U>& operator()(const T& key, const U& val) {
     m_map[key] = val;
     return *this;
   }
   operator map<T, U>() { return m_map; }
 
 };
 
 //==========================================================================
 
 // Template to make initializing maps simpler, while not relying on C++11.
 // Usage: map(createmap<T,U>(a,b)(c,d)(e,f));
 
 template <typename T, typename U> class create_unordered_map {
 
 private:
 
   unordered_map<T, U> um_map;
 
 public:
 
   create_unordered_map(const T& key, const U& val) { um_map[key] = val; }
   create_unordered_map<T, U>& operator()(const T& key, const U& val) {
     um_map[key] = val;
     return *this;
   }
   operator unordered_map<T, U>() { return um_map; }
 
 };
 
 //==========================================================================
 
 // Template to make initializing maps simpler, while not relying on C++11.
 // Usage: map(createmap<T,U>(a,b)(c,d)(e,f));
 
 template <typename T> class createvector {
 
 private:
 
   vector<T> m_vector;
 
 public:
 
   createvector(const T& val) { m_vector.push_back(val); }
   createvector<T>& operator()(const T& val) {
     m_vector.push_back(val);
     return *this;
   }
   operator vector<T>() { return m_vector; }
 
 };
 
 //==========================================================================
 
 // Helper function to calculate dilogarithm.
 
 double polev(double x,double* coef,int N );
 // Function to calculate dilogarithm.
 double dilog(double x);
 
 //==========================================================================
 
 // Kallen function and derived quantities.
 
 double lABC(double a, double b, double c);
 double bABC(double a, double b, double c);
 double gABC(double a, double b, double c);
 
 //==========================================================================
 
 class DireFunction {
 
 public:
 
   virtual double f(double) { return 0.; }
   virtual double f(double, double) { return 0.; }
   virtual double f(double, vector<double>) { return 0.; }
 
 };
 
 class DireRootFinder {
 
 public:
 
   DireRootFinder() {};
   virtual ~DireRootFinder() {};
 
   double findRootSecant1D( DireFunction* f, double xmin, double xmax,
     double constant, vector<double> xb = vector<double>(), int N=10 ) {
     vector<double> x;
     x.push_back(xmin);
     x.push_back(xmax);
     for ( int i=2; i < N; ++i ) {
       double xn = x[i-1]
       - ( f->f(x[i-1],xb) - constant)
       * ( x[i-1] - x[i-2] )
       / ( f->f(x[i-1],xb) - f->f(x[i-2],xb) );
       x.push_back(xn);
     }
     return x.back();
   }
 
   double findRoot1D( DireFunction* f, double xmin, double xmax,
     double constant, vector<double> xx = vector<double>(), int N=10,
     double tol = 1e-10 ) {
 
     double a(xmin), b(xmax), c(xmax), d(0.), e(0.),
       fa(f->f(a,xx)-constant), fb(f->f(b,xx)-constant), fc(fb),
       p(0.), q(0.), r(0.), s(0.),
       tol1(tol), xm(0.);
     double EPS = numeric_limits<double>::epsilon();
 
     // No root.
     if ( (fa>0. && fb>0.) || (fa<0. && fb<0.) ) {
      cout << "no root " << constant << " " << f->f(a,xx) << " " << f->f(b,xx)
      << endl;
      return numeric_limits<double>::quiet_NaN();
     }
 
     for ( int i=0; i < N; ++i ) {
 
       if ( (fb>0. && fc>0.) || (fb<0. && fc<0.) ) {
         c  = a;
         fc = fa;
         e  = d = b-a;
       }
 
       if ( abs(fc) < abs(fb) ) {
         a = b;
         b = c;
         c = a;
         fa = fb;
         fb = fc;
         fc = fa;
       }
 
       tol1 = 2.*EPS*abs(b) + 0.5*tol;
       xm = 0.5*(c-b);
 
       if (abs(xm) <= tol1 || fb == 0.) return b;
 
       if (abs(e) >= tol1 && abs(fa) > abs(fb) ) {
         s = fb/fa;
         if ( a == c ) {
           p = 2.*xm*s;
           q = 1.-s;
         } else {
           q = fa/fc;
           r = fb/fc;
           p = s*(2.*xm*q*(q-r) - (b-a)*(r-1.));
           q = (q-1.)*(r-1.)*(s-1.);
         }
         if (p>0.) q = -q;
         p = abs(p);
         double min1 = 3.*xm*q - abs(tol1*q);
         double min2 = abs(e*q);
         if (2.*p < ((min1 < min2) ? min1 : min2)) {
           e = d;
           d = p/q;
         } else {
           d = xm;
           e = d;
         }
 
       } else {
         d = xm;
         e = d;
       }
 
       a = b;
       fa = fb;
 
       if (abs(d) > tol1) { b += d; }
       else {
         b += (xm> 0.) ? tol1 : -tol1;
       }
       fb = f->f(b,xx)-constant;
     }
 
     // Failed. Return NaN
     return numeric_limits<double>::quiet_NaN();
 
   }
 
 };
 
 //==========================================================================
 
 class DireEventInfo {
 
   public:
 
   DireEventInfo() {}
 
   // Bookkeeping of soft particles.
   int sizeSoftPos () const { return softPosSave.size(); }
   int getSoftPos(unsigned int i) const {
     return (i > softPosSave.size()-1) ? -1 : softPosSave[i]; }
   bool isSoft(int iPos) {
     vector<int>::iterator it = find( softPosSave.begin(),
       softPosSave.end(), iPos);
     return (it != softPosSave.end());
   }
   void addSoftPos(int iPos) { if (!isSoft(iPos)) softPosSave.push_back(iPos); }
   void removeSoftPos(int iPos) {
     vector<int>::iterator it = find( softPosSave.begin(),
       softPosSave.end(), iPos);
     if (it != softPosSave.end()) softPosSave.erase(it);
   }
   void updateSoftPos(int iPosOld, int iPosNew) {
     if (isSoft(iPosOld)) removeSoftPos(iPosOld);
     if (isSoft(iPosNew)) removeSoftPos(iPosNew);
     addSoftPos(iPosNew);
   }
   void updateSoftPosIfMatch(int iPosOld, int iPosNew) {
     if (isSoft(iPosOld)) {
       vector<int>::iterator it
         = find (softPosSave.begin(), softPosSave.end(), iPosOld);
       *it = iPosNew;
     }
   }
   vector<int> softPos () const { return softPosSave; }
   void clearSoftPos () { softPosSave.clear(); }
   void listSoft() const {
     cout << " 'Soft' particles: ";
     for (int i=0; i < sizeSoftPos(); ++i) cout << setw(5) << getSoftPos(i);
     cout << endl;
   }
 
   // Bookkeeping of resonances.
   void removeResPos(int iPos) {
     vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPos);
     if (it == iPosRes.end()) return;
     iPosRes.erase(it);
     sort (iPosRes.begin(), iPosRes.end());
   }
   void addResPos(int iPos) {
     vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPos);
     if (it != iPosRes.end()) return;
     iPosRes.push_back(iPos);
     sort (iPosRes.begin(), iPosRes.end());
   }
   void updateResPos(int iPosOld, int iPosNew) {
     vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPosOld);
     if (it == iPosRes.end()) iPosRes.push_back(iPosNew);
     else                    *it = iPosNew;
     sort (iPosRes.begin(), iPosRes.end());
   }
   void updateResPosIfMatch(int iPosOld, int iPosNew) {
     vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPosOld);
     if (it != iPosRes.end()) {
       iPosRes.erase(it);
       iPosRes.push_back(iPosNew);
       sort (iPosRes.begin(), iPosRes.end());
     }
   }
   bool isRes(int iPos) {
     vector<int>::iterator it = find (iPosRes.begin(), iPosRes.end(), iPos);
     return (it != iPosRes.end());
   }
   int sizeResPos() const { return iPosRes.size(); }
   int getResPos(unsigned int i) const {
     return (i > iPosRes.size()-1) ? -1 : iPosRes[i]; }
   void clearResPos() { iPosRes.resize(0); }
   void listRes() const {
     cout << " 'Resonant' particles: ";
     for (int i=0; i < sizeResPos(); ++i) cout << setw(5) <<  getResPos(i);
     cout << endl;
   }
 
   // Data members.
   vector<int> softPosSave;
   vector<int> iPosRes;
 
 };
 
 
 //==========================================================================
 
 class DireDebugInfo {
 
   public:
 
   DireDebugInfo() = default;
 
   void clearMessages() {
     messageStream0.str("");
     messageStream1.str("");
     messageStream2.str("");
   }
 
   void printMessages( int verbosity = 0) {
     cout << "\n"
       << "*------------------------------------------------------------*\n"
       << "*----------------- Begin diagnostic output ------------------*\n\n";
     if (verbosity == 0) cout << scientific << setprecision(8)
     << messageStream0.str();
     if (verbosity == 1) cout << scientific << setprecision(8)
     << messageStream1.str();
     if (verbosity == 2) cout << scientific << setprecision(8)
     << messageStream2.str();
     cout << "\n\n"
       << "*----------------- End diagnostic output -------------------*\n"
       << "*-----------------------------------------------------------*"
       << endl;
   }
 
   void printHistograms() {}
 
   void fillHistograms(int type, int nfinal, double mec, double pT, double z) {
     if (false) cout << type*nfinal*mec*pT*z;}
 
   // Add debug messages to message stream.
   ostream & message ( int verbosity = 0) {
     if (verbosity == 0) return messageStream0;
     if (verbosity == 1) return messageStream1;
     if (verbosity == 2) return messageStream2;
     return messageStream0;
   }
 
   // Debug message streams.
   ostringstream messageStream0, messageStream1, messageStream2;
 
 };
 
 //==========================================================================
 
 class DireInfo {
 
   public:
 
   DireInfo() {}
 
   void clearAll() {
     direEventInfo.clearResPos();
     direEventInfo.clearSoftPos();
     direDebugInfo.clearMessages();
   }
 
   // Resonance info forwards.
   void removeResPos(int iPos)   { return direEventInfo.removeResPos(iPos); }
   void addResPos(int iPos)      { return direEventInfo.addResPos(iPos); }
   bool isRes(int iPos)          { return direEventInfo.isRes(iPos); }
   void clearResPos ()           { return direEventInfo.clearResPos(); }
   int sizeResPos() const        { return direEventInfo.sizeResPos(); }
   void listRes() const          { return direEventInfo.listRes(); }
   int getResPos(unsigned int i) const { return direEventInfo.getResPos(i); }
   void updateResPos(int iPosOld, int iPosNew) {
     return direEventInfo.updateResPos(iPosOld,iPosNew); }
   void updateResPosIfMatch(int iPosOld, int iPosNew) {
     return direEventInfo.updateResPosIfMatch(iPosOld,iPosNew); }
 
   // Debug info forwards.
   void printMessages( int verbosity = 0) {
     return direDebugInfo.printMessages(verbosity); }
   ostream & message ( int verbosity = 0) {
     return direDebugInfo.message(verbosity); }
   void clearMessages() { direDebugInfo.clearMessages(); }
 
   void fillHistograms(int type, int nfinal, double mec, double pT, double z) {
     direDebugInfo.fillHistograms(type, nfinal,  mec, pT, z); }
   void printHistograms() { direDebugInfo.printHistograms(); }
 
   // Soft particle info forwards.
   bool isSoft(int iPos)          { return direEventInfo.isSoft(iPos); }
   void addSoftPos(int iPos)      { return direEventInfo.addSoftPos(iPos); }
   void removeSoftPos(int iPos)   { return direEventInfo.removeSoftPos(iPos); }
   vector<int> softPos ()         { return direEventInfo.softPos(); }
   void clearSoftPos ()           { return direEventInfo.clearSoftPos(); }
   int sizeSoftPos () const       { return direEventInfo.sizeSoftPos(); }
   void listSoft() const          { return direEventInfo.listSoft(); }
   int getSoftPos(unsigned int i) const { return direEventInfo.getSoftPos(i); }
   void updateSoftPos(int iPosOld, int iPosNew) {
     return direEventInfo.updateSoftPos(iPosOld, iPosNew);
   }
   void updateSoftPosIfMatch(int iPosOld, int iPosNew) {
     return direEventInfo.updateSoftPosIfMatch(iPosOld, iPosNew);
   }
 
   DireEventInfo direEventInfo;
   DireDebugInfo direDebugInfo;
 
 };
 
 //==========================================================================
 
 } // end namespace Pythia8
 
 #endif // Pythia8_DireBasics_H

This page was automatically generated by the 2.3.7 LXR engine. The LXR team