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