[1] | 1 | // HadronScatter.h is a part of the PYTHIA event generator. |
---|
| 2 | // Copyright (C) 2012 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 | #ifndef Pythia8_HadronScatter_H |
---|
| 7 | #define Pythia8_HadronScatter_H |
---|
| 8 | |
---|
| 9 | #include "Event.h" |
---|
| 10 | #include "Info.h" |
---|
| 11 | #include "ParticleData.h" |
---|
| 12 | #include "PythiaComplex.h" |
---|
| 13 | |
---|
| 14 | #include <algorithm> |
---|
| 15 | #include <map> |
---|
| 16 | #include <set> |
---|
| 17 | using std::map; |
---|
| 18 | using std::set; |
---|
| 19 | using std::sort; |
---|
| 20 | |
---|
| 21 | namespace Pythia8 { |
---|
| 22 | |
---|
| 23 | class SigmaPartialWave { |
---|
| 24 | public: |
---|
| 25 | // Initialisation |
---|
| 26 | bool init(int, string, string, Info *, ParticleData *, Rndm *); |
---|
| 27 | |
---|
| 28 | // Read data file |
---|
| 29 | bool readFile(string, string); |
---|
| 30 | |
---|
| 31 | // Set the subprocess/incoming particles |
---|
| 32 | bool setSubprocess(int); |
---|
| 33 | bool setSubprocess(int, int); |
---|
| 34 | |
---|
| 35 | // Return sigma total/elastic, dSigma/dCos(theta) |
---|
| 36 | double sigmaEl(double Wcm) { return sigma(0, Wcm); } |
---|
| 37 | double sigmaTot(double Wcm) { return sigma(1, Wcm); } |
---|
| 38 | double dSigma(double Wcm, double cTheta) { return sigma(2, Wcm, cTheta); } |
---|
| 39 | |
---|
| 40 | // Return a cos(theta) value |
---|
| 41 | double pickCosTheta(double); |
---|
| 42 | |
---|
| 43 | // Return maximum sigma elastic |
---|
| 44 | double getSigmaElMax() { return sigElMax; } |
---|
| 45 | |
---|
| 46 | private: |
---|
| 47 | // Pointers |
---|
| 48 | Info *infoPtr; |
---|
| 49 | ParticleData *particleDataPtr; |
---|
| 50 | Rndm *rndmPtr; |
---|
| 51 | |
---|
| 52 | // Constants |
---|
| 53 | static const int LSHIFT, ISHIFT, SUBBIN, ITER; |
---|
| 54 | static const double CONVERT2MB, WCMBIN, CTBIN, MASSSAFETY, GRIDSAFETY; |
---|
| 55 | static const complex BINEND; |
---|
| 56 | |
---|
| 57 | // Settings |
---|
| 58 | int process, subprocess, subprocessMax, norm; |
---|
| 59 | |
---|
| 60 | // Current incoming and maximum L/I values |
---|
| 61 | int idA, idB, Lmax, Imax; |
---|
| 62 | |
---|
| 63 | // Masses of incoming, last bin, maximum sigma elastic |
---|
| 64 | double mA, mB, binMax, sigElMax; |
---|
| 65 | |
---|
| 66 | // Map subprocess to incoming and vice versa: |
---|
| 67 | // sp2in[subprocess] = idA, idB |
---|
| 68 | // in2sp[idA, idB] = subprocess |
---|
| 69 | map < int, pair < int, int > > sp2in; |
---|
| 70 | map < pair < int, int >, int > in2sp; |
---|
| 71 | |
---|
| 72 | // Isospin coefficients isoCoeff[subprocess][2I] |
---|
| 73 | map < int, map < int, double > > isoCoeff; |
---|
| 74 | |
---|
| 75 | // Storage for partial wave data: |
---|
| 76 | // pwData[L * LSHIFT + 2I * ISHIFT + 2J][eNow] = T |
---|
| 77 | map < int, map < double, complex > > pwData; |
---|
| 78 | |
---|
| 79 | // Values of Pl and Pl' as computed by legendreP |
---|
| 80 | vector < double > PlVec, PlpVec; |
---|
| 81 | |
---|
| 82 | // Integration grid [subprocess][WcmBin][cosThetaBin] |
---|
| 83 | vector < vector < vector < double > > > gridMax; |
---|
| 84 | vector < vector < double > > gridNorm; |
---|
| 85 | |
---|
| 86 | // Setup subprocesses (including isospin coefficients) |
---|
| 87 | void setupSubprocesses(); |
---|
| 88 | |
---|
| 89 | // Setup grids for integration |
---|
| 90 | void setupGrid(); |
---|
| 91 | |
---|
| 92 | // Routine for calculating sigma total/elastic and dSigma/dCos(theta) |
---|
| 93 | double sigma(int, double, double = 0.); |
---|
| 94 | |
---|
| 95 | // Generate Legendre polynomials (and optionally derivatives) |
---|
| 96 | void legendreP(double, bool = false); |
---|
| 97 | |
---|
| 98 | }; |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | //========================================================================== |
---|
| 102 | |
---|
| 103 | // HadronScatterPair class |
---|
| 104 | // Simple class to hold details of a pair of hadrons which will scatter. |
---|
| 105 | // Stores indices in event record and the measure used for ordering |
---|
| 106 | |
---|
| 107 | // Store a pair of indices |
---|
| 108 | typedef pair < int, int > HSIndex; |
---|
| 109 | |
---|
| 110 | class HadronScatterPair { |
---|
| 111 | public: |
---|
| 112 | // Constructor |
---|
| 113 | HadronScatterPair(const HSIndex &i1in, int yt1in, int pt1in, |
---|
| 114 | const HSIndex &i2in, int yt2in, int pt2in, |
---|
| 115 | double measureIn) : |
---|
| 116 | i1(i1in), yt1(yt1in), pt1(pt1in), |
---|
| 117 | i2(i2in), yt2(yt2in), pt2(pt2in), |
---|
| 118 | measure(measureIn) {} |
---|
| 119 | |
---|
| 120 | // Operator for sorting according to ordering measure |
---|
| 121 | bool operator<(const HadronScatterPair& in) const { |
---|
| 122 | return this->measure < in.measure; |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | // Indices into event record of hadrons to scatter |
---|
| 126 | HSIndex i1; |
---|
| 127 | int yt1, pt1; |
---|
| 128 | HSIndex i2; |
---|
| 129 | int yt2, pt2; |
---|
| 130 | // Ordering measure |
---|
| 131 | double measure; |
---|
| 132 | }; |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | //========================================================================== |
---|
| 136 | |
---|
| 137 | // HadronScatter class |
---|
| 138 | |
---|
| 139 | class HadronScatter { |
---|
| 140 | |
---|
| 141 | public: |
---|
| 142 | |
---|
| 143 | // Constructor. |
---|
| 144 | HadronScatter() {} |
---|
| 145 | |
---|
| 146 | // Initialisation |
---|
| 147 | bool init(Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn, |
---|
| 148 | ParticleData *particleDataPtr); |
---|
| 149 | |
---|
| 150 | // Perform all hadron scatterings |
---|
| 151 | void scatter(Event&); |
---|
| 152 | |
---|
| 153 | private: |
---|
| 154 | |
---|
| 155 | // Pointer to various information on the generation. |
---|
| 156 | Info* infoPtr; |
---|
| 157 | Rndm* rndmPtr; |
---|
| 158 | |
---|
| 159 | // Settings |
---|
| 160 | bool doHadronScatter, afterDecay, allowDecayProd, |
---|
| 161 | scatterRepeat, doTile; |
---|
| 162 | int hadronSelect, scatterProb; |
---|
| 163 | double Npar, kPar, pPar, jPar, rMax, rMax2; |
---|
| 164 | double pTsigma, pTsigma2, pT0MPI; |
---|
| 165 | |
---|
| 166 | // Tiling |
---|
| 167 | int ytMax, ptMax; |
---|
| 168 | double yMin, yMax, ytSize, ptSize; |
---|
| 169 | vector < vector < set < HSIndex > > > tile; |
---|
| 170 | |
---|
| 171 | // Keep track of scattered pairs |
---|
| 172 | set < HSIndex > scattered; |
---|
| 173 | |
---|
| 174 | // Partial wave amplitudes |
---|
| 175 | SigmaPartialWave sigmaPW[3]; |
---|
| 176 | |
---|
| 177 | // Maximum sigma elastic |
---|
| 178 | double sigElMax; |
---|
| 179 | |
---|
| 180 | // Decide if a hadron can scatter |
---|
| 181 | bool canScatter(Event &, int); |
---|
| 182 | |
---|
| 183 | // Probability for a pair of hadrons to scatter |
---|
| 184 | bool doesScatter(Event &, const HSIndex &, const HSIndex &); |
---|
| 185 | |
---|
| 186 | // Calculate measure for ordering of scatterings |
---|
| 187 | double measure(Event &, int, int); |
---|
| 188 | |
---|
| 189 | // Perform a single hadron scattering |
---|
| 190 | bool hadronScatter(Event &, HadronScatterPair &); |
---|
| 191 | |
---|
| 192 | // Tiling functions |
---|
| 193 | bool tileIntProb(vector < HadronScatterPair > &, Event &, |
---|
| 194 | const HSIndex &, int, int, bool); |
---|
| 195 | int yTile(Event& event, int idx) { |
---|
| 196 | return (doTile) ? int((event[idx].y() - yMin) / ytSize) : 0; |
---|
| 197 | } |
---|
| 198 | int pTile(Event& event, int idx) { |
---|
| 199 | return (doTile) ? int((event[idx].phi() + M_PI) / ptSize) : 0; |
---|
| 200 | } |
---|
| 201 | |
---|
| 202 | // Debug |
---|
| 203 | void debugOutput(); |
---|
| 204 | }; |
---|
| 205 | |
---|
| 206 | //========================================================================== |
---|
| 207 | |
---|
| 208 | |
---|
| 209 | } // end namespace Pythia8 |
---|
| 210 | |
---|
| 211 | #endif // Pythia8_HadronScatter_H |
---|
| 212 | |
---|