1 | // PhaseSpace.cc 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 | // Function definitions (not found in the header) for the |
---|
7 | // PhaseSpace and PhaseSpace2to2tauyz classes. |
---|
8 | |
---|
9 | #include "PhaseSpace.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // The PhaseSpace class. |
---|
16 | // Base class for phase space generators. |
---|
17 | |
---|
18 | //-------------------------------------------------------------------------- |
---|
19 | |
---|
20 | // Constants: could be changed here if desired, but normally should not. |
---|
21 | // These are of technical nature, as described for each. |
---|
22 | |
---|
23 | // Number of trial maxima around which maximum search is performed. |
---|
24 | const int PhaseSpace::NMAXTRY = 2; |
---|
25 | |
---|
26 | // Number of three-body trials in phase space optimization. |
---|
27 | const int PhaseSpace::NTRY3BODY = 20; |
---|
28 | |
---|
29 | // Maximum cross section increase, just in case true maximum not found. |
---|
30 | const double PhaseSpace::SAFETYMARGIN = 1.05; |
---|
31 | |
---|
32 | // Small number to avoid division by zero. |
---|
33 | const double PhaseSpace::TINY = 1e-20; |
---|
34 | |
---|
35 | // Fraction of total weight that is shared evenly between all shapes. |
---|
36 | const double PhaseSpace::EVENFRAC = 0.4; |
---|
37 | |
---|
38 | // Two cross sections with a small relative error are assumed same. |
---|
39 | const double PhaseSpace::SAMESIGMA = 1e-6; |
---|
40 | |
---|
41 | // Do not include resonances peaked too far outside allowed mass region. |
---|
42 | const double PhaseSpace::WIDTHMARGIN = 20.; |
---|
43 | |
---|
44 | // Special optimization treatment when two resonances at almost same mass. |
---|
45 | const double PhaseSpace::SAMEMASS = 0.01; |
---|
46 | |
---|
47 | // Minimum phase space left when kinematics constraints are combined. |
---|
48 | const double PhaseSpace::MASSMARGIN = 0.01; |
---|
49 | |
---|
50 | // When using Breit-Wigners in 2 -> 2 raise maximum weight estimate. |
---|
51 | const double PhaseSpace::EXTRABWWTMAX = 1.25; |
---|
52 | |
---|
53 | // Size of Breit-Wigner threshold region, for mass selection biasing. |
---|
54 | const double PhaseSpace::THRESHOLDSIZE = 3.; |
---|
55 | |
---|
56 | // Step size in optimal-mass search, for mass selection biasing. |
---|
57 | const double PhaseSpace::THRESHOLDSTEP = 0.2; |
---|
58 | |
---|
59 | // Minimal rapidity range for allowed open range (in 2 -> 3). |
---|
60 | const double PhaseSpace::YRANGEMARGIN = 1e-6; |
---|
61 | |
---|
62 | // Cutoff for f_e^e at x < 1 - 10^{-10} to be used in phase space selection. |
---|
63 | // Note: the ...MIN quantities come from 1 - x_max or 1 - tau_max. |
---|
64 | const double PhaseSpace::LEPTONXMIN = 1e-10; |
---|
65 | const double PhaseSpace::LEPTONXMAX = 1. - 1e-10; |
---|
66 | const double PhaseSpace::LEPTONXLOGMIN = log(1e-10); |
---|
67 | const double PhaseSpace::LEPTONXLOGMAX = log(1. - 1e-10); |
---|
68 | const double PhaseSpace::LEPTONTAUMIN = 2e-10; |
---|
69 | |
---|
70 | // Safety to avoid division with unreasonably small value for z selection. |
---|
71 | const double PhaseSpace::SHATMINZ = 1.; |
---|
72 | |
---|
73 | // Regularization for small pT2min in z = cos(theta) selection. |
---|
74 | const double PhaseSpace::PT2RATMINZ = 0.0001; |
---|
75 | |
---|
76 | // These numbers are hardwired empirical parameters, |
---|
77 | // intended to speed up the M-generator. |
---|
78 | const double PhaseSpace::WTCORRECTION[11] = { 1., 1., 1., |
---|
79 | 2., 5., 15., 60., 250., 1250., 7000., 50000. }; |
---|
80 | |
---|
81 | // MBR: form factor appoximation with two exponents, [FFB1,FFB2] = GeV-2. |
---|
82 | const double PhaseSpace::FFA1 = 0.9; |
---|
83 | const double PhaseSpace::FFA2 = 0.1; |
---|
84 | const double PhaseSpace::FFB1 = 4.6; |
---|
85 | const double PhaseSpace::FFB2 = 0.6; |
---|
86 | |
---|
87 | //-------------------------------------------------------------------------- |
---|
88 | |
---|
89 | // Perform simple initialization and store pointers. |
---|
90 | |
---|
91 | void PhaseSpace::init(bool isFirst, SigmaProcess* sigmaProcessPtrIn, |
---|
92 | Info* infoPtrIn, Settings* settingsPtrIn, ParticleData* particleDataPtrIn, |
---|
93 | Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, |
---|
94 | Couplings* couplingsPtrIn, SigmaTotal* sigmaTotPtrIn, |
---|
95 | UserHooks* userHooksPtrIn) { |
---|
96 | |
---|
97 | // Store input pointers for future use. |
---|
98 | sigmaProcessPtr = sigmaProcessPtrIn; |
---|
99 | infoPtr = infoPtrIn; |
---|
100 | settingsPtr = settingsPtrIn; |
---|
101 | particleDataPtr = particleDataPtrIn; |
---|
102 | rndmPtr = rndmPtrIn; |
---|
103 | beamAPtr = beamAPtrIn; |
---|
104 | beamBPtr = beamBPtrIn; |
---|
105 | couplingsPtr = couplingsPtrIn; |
---|
106 | sigmaTotPtr = sigmaTotPtrIn; |
---|
107 | userHooksPtr = userHooksPtrIn; |
---|
108 | |
---|
109 | // Some commonly used beam information. |
---|
110 | idA = beamAPtr->id(); |
---|
111 | idB = beamBPtr->id(); |
---|
112 | mA = beamAPtr->m(); |
---|
113 | mB = beamBPtr->m(); |
---|
114 | eCM = infoPtr->eCM(); |
---|
115 | s = eCM * eCM; |
---|
116 | |
---|
117 | // Flag if lepton beams, and if non-resolved ones. |
---|
118 | hasLeptonBeams = ( beamAPtr->isLepton() || beamBPtr->isLepton() ); |
---|
119 | hasPointLeptons = ( hasLeptonBeams |
---|
120 | && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) ); |
---|
121 | |
---|
122 | // Standard phase space cuts. |
---|
123 | if (isFirst || settingsPtr->flag("PhaseSpace:sameForSecond")) { |
---|
124 | mHatGlobalMin = settingsPtr->parm("PhaseSpace:mHatMin"); |
---|
125 | mHatGlobalMax = settingsPtr->parm("PhaseSpace:mHatMax"); |
---|
126 | pTHatGlobalMin = settingsPtr->parm("PhaseSpace:pTHatMin"); |
---|
127 | pTHatGlobalMax = settingsPtr->parm("PhaseSpace:pTHatMax"); |
---|
128 | |
---|
129 | // Optionally separate phase space cuts for second hard process. |
---|
130 | } else { |
---|
131 | mHatGlobalMin = settingsPtr->parm("PhaseSpace:mHatMinSecond"); |
---|
132 | mHatGlobalMax = settingsPtr->parm("PhaseSpace:mHatMaxSecond"); |
---|
133 | pTHatGlobalMin = settingsPtr->parm("PhaseSpace:pTHatMinSecond"); |
---|
134 | pTHatGlobalMax = settingsPtr->parm("PhaseSpace:pTHatMaxSecond"); |
---|
135 | } |
---|
136 | |
---|
137 | // Cutoff against divergences at pT -> 0. |
---|
138 | pTHatMinDiverge = settingsPtr->parm("PhaseSpace:pTHatMinDiverge"); |
---|
139 | |
---|
140 | // When to use Breit-Wigners. |
---|
141 | useBreitWigners = settingsPtr->flag("PhaseSpace:useBreitWigners"); |
---|
142 | minWidthBreitWigners = settingsPtr->parm("PhaseSpace:minWidthBreitWigners"); |
---|
143 | |
---|
144 | // Whether generation is with variable energy. |
---|
145 | doEnergySpread = settingsPtr->flag("Beams:allowMomentumSpread"); |
---|
146 | |
---|
147 | // Flags for maximization information and violation handling. |
---|
148 | showSearch = settingsPtr->flag("PhaseSpace:showSearch"); |
---|
149 | showViolation = settingsPtr->flag("PhaseSpace:showViolation"); |
---|
150 | increaseMaximum = settingsPtr->flag("PhaseSpace:increaseMaximum"); |
---|
151 | |
---|
152 | // Know whether a Z0 is pure Z0 or admixed with gamma*. |
---|
153 | gmZmodeGlobal = settingsPtr->mode("WeakZ0:gmZmode"); |
---|
154 | |
---|
155 | // Flags if user should be allowed to reweight cross section. |
---|
156 | canModifySigma = (userHooksPtr != 0) |
---|
157 | ? userHooksPtr->canModifySigma() : false; |
---|
158 | canBiasSelection = (userHooksPtr != 0) |
---|
159 | ? userHooksPtr->canBiasSelection() : false; |
---|
160 | |
---|
161 | // Parameters for simplified reweighting of 2 -> 2 processes. |
---|
162 | canBias2Sel = settingsPtr->flag("PhaseSpace:bias2Selection"); |
---|
163 | bias2SelPow = settingsPtr->parm("PhaseSpace:bias2SelectionPow"); |
---|
164 | bias2SelRef = settingsPtr->parm("PhaseSpace:bias2SelectionRef"); |
---|
165 | |
---|
166 | // Default event-specific kinematics properties. |
---|
167 | x1H = 1.; |
---|
168 | x2H = 1.; |
---|
169 | m3 = 0.; |
---|
170 | m4 = 0.; |
---|
171 | m5 = 0.; |
---|
172 | s3 = m3 * m3; |
---|
173 | s4 = m4 * m4; |
---|
174 | s5 = m5 * m5; |
---|
175 | mHat = eCM; |
---|
176 | sH = s; |
---|
177 | tH = 0.; |
---|
178 | uH = 0.; |
---|
179 | pTH = 0.; |
---|
180 | theta = 0.; |
---|
181 | phi = 0.; |
---|
182 | runBW3H = 1.; |
---|
183 | runBW4H = 1.; |
---|
184 | runBW5H = 1.; |
---|
185 | |
---|
186 | // Default cross section information. |
---|
187 | sigmaNw = 0.; |
---|
188 | sigmaMx = 0.; |
---|
189 | sigmaPos = 0.; |
---|
190 | sigmaNeg = 0.; |
---|
191 | newSigmaMx = false; |
---|
192 | biasWt = 1.; |
---|
193 | |
---|
194 | } |
---|
195 | |
---|
196 | //-------------------------------------------------------------------------- |
---|
197 | |
---|
198 | // Allow for nonisotropic decays when ME's available. |
---|
199 | |
---|
200 | void PhaseSpace::decayKinematics( Event& process) { |
---|
201 | |
---|
202 | // Identify sets of sister partons. |
---|
203 | int iResEnd = 4; |
---|
204 | for (int iResBeg = 5; iResBeg < process.size(); ++iResBeg) { |
---|
205 | if (iResBeg <= iResEnd) continue; |
---|
206 | iResEnd = iResBeg; |
---|
207 | while ( iResEnd < process.size() - 1 |
---|
208 | && process[iResEnd + 1].mother1() == process[iResBeg].mother1() |
---|
209 | && process[iResEnd + 1].mother2() == process[iResBeg].mother2() ) |
---|
210 | ++iResEnd; |
---|
211 | |
---|
212 | // Check that at least one of them is a resonance. |
---|
213 | bool hasRes = false; |
---|
214 | for (int iRes = iResBeg; iRes <= iResEnd; ++iRes) |
---|
215 | if ( !process[iRes].isFinal() ) hasRes = true; |
---|
216 | if ( !hasRes ) continue; |
---|
217 | |
---|
218 | // Evaluate matrix element and decide whether to keep kinematics. |
---|
219 | double decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd); |
---|
220 | if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay" |
---|
221 | "Kinematics: negative angular weight"); |
---|
222 | if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay" |
---|
223 | "Kinematics: angular weight above unity"); |
---|
224 | while (decWt < rndmPtr->flat() ) { |
---|
225 | |
---|
226 | // Find resonances for which to redo decay angles. |
---|
227 | for (int iRes = iResBeg; iRes < process.size(); ++iRes) { |
---|
228 | if ( process[iRes].isFinal() ) continue; |
---|
229 | int iResMother = iRes; |
---|
230 | while (iResMother > iResEnd) |
---|
231 | iResMother = process[iResMother].mother1(); |
---|
232 | if (iResMother < iResBeg) continue; |
---|
233 | |
---|
234 | // Do decay of this mother isotropically in phase space. |
---|
235 | decayKinematicsStep( process, iRes); |
---|
236 | |
---|
237 | // End loop over resonance decay chains. |
---|
238 | } |
---|
239 | |
---|
240 | // Ready to allow new test of matrix element. |
---|
241 | decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd); |
---|
242 | if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay" |
---|
243 | "Kinematics: negative angular weight"); |
---|
244 | if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay" |
---|
245 | "Kinematics: angular weight above unity"); |
---|
246 | } |
---|
247 | |
---|
248 | // End loop over sets of sister resonances/partons. |
---|
249 | } |
---|
250 | |
---|
251 | } |
---|
252 | |
---|
253 | //-------------------------------------------------------------------------- |
---|
254 | |
---|
255 | // Reselect decay products momenta isotropically in phase space. |
---|
256 | // Does not redo secondary vertex position! |
---|
257 | |
---|
258 | void PhaseSpace::decayKinematicsStep( Event& process, int iRes) { |
---|
259 | |
---|
260 | // Multiplicity and mother mass and four-momentum. |
---|
261 | int i1 = process[iRes].daughter1(); |
---|
262 | int mult = process[iRes].daughter2() + 1 - i1; |
---|
263 | double m0 = process[iRes].m(); |
---|
264 | Vec4 pRes = process[iRes].p(); |
---|
265 | |
---|
266 | // Description of two-body decays as simple special case. |
---|
267 | if (mult == 2) { |
---|
268 | |
---|
269 | // Products and product masses. |
---|
270 | int i2 = i1 + 1; |
---|
271 | double m1t = process[i1].m(); |
---|
272 | double m2t = process[i2].m(); |
---|
273 | |
---|
274 | // Energies and absolute momentum in the rest frame. |
---|
275 | double e1 = 0.5 * (m0*m0 + m1t*m1t - m2t*m2t) / m0; |
---|
276 | double e2 = 0.5 * (m0*m0 + m2t*m2t - m1t*m1t) / m0; |
---|
277 | double p12 = 0.5 * sqrtpos( (m0 - m1t - m2t) * (m0 + m1t + m2t) |
---|
278 | * (m0 + m1t - m2t) * (m0 - m1t + m2t) ) / m0; |
---|
279 | |
---|
280 | // Pick isotropic angles to give three-momentum. |
---|
281 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
282 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
283 | double phi12 = 2. * M_PI * rndmPtr->flat(); |
---|
284 | double pX = p12 * sinTheta * cos(phi12); |
---|
285 | double pY = p12 * sinTheta * sin(phi12); |
---|
286 | double pZ = p12 * cosTheta; |
---|
287 | |
---|
288 | // Fill four-momenta in mother rest frame and then boost to lab frame. |
---|
289 | Vec4 p1( pX, pY, pZ, e1); |
---|
290 | Vec4 p2( -pX, -pY, -pZ, e2); |
---|
291 | p1.bst( pRes ); |
---|
292 | p2.bst( pRes ); |
---|
293 | |
---|
294 | // Done for two-body decay. |
---|
295 | process[i1].p( p1 ); |
---|
296 | process[i2].p( p2 ); |
---|
297 | return; |
---|
298 | } |
---|
299 | |
---|
300 | // Description of three-body decays as semi-simple special case. |
---|
301 | if (mult == 3) { |
---|
302 | |
---|
303 | // Products and product masses. |
---|
304 | int i2 = i1 + 1; |
---|
305 | int i3 = i2 + 1; |
---|
306 | double m1t = process[i1].m(); |
---|
307 | double m2t = process[i2].m(); |
---|
308 | double m3t = process[i3].m(); |
---|
309 | double mDiff = m0 - (m1t + m2t + m3t); |
---|
310 | |
---|
311 | // Kinematical limits for 2+3 mass. Maximum phase-space weight. |
---|
312 | double m23Min = m2t + m3t; |
---|
313 | double m23Max = m0 - m1t; |
---|
314 | double p1Max = 0.5 * sqrtpos( (m0 - m1t - m23Min) |
---|
315 | * (m0 + m1t + m23Min) * (m0 + m1t - m23Min) |
---|
316 | * (m0 - m1t + m23Min) ) / m0; |
---|
317 | double p23Max = 0.5 * sqrtpos( (m23Max - m2t - m3t) |
---|
318 | * (m23Max + m2t + m3t) * (m23Max + m2t - m3t) |
---|
319 | * (m23Max - m2t + m3t) ) / m23Max; |
---|
320 | double wtPSmax = 0.5 * p1Max * p23Max; |
---|
321 | |
---|
322 | // Pick an intermediate mass m23 flat in the allowed range. |
---|
323 | double wtPS, m23, p1Abs, p23Abs; |
---|
324 | do { |
---|
325 | m23 = m23Min + rndmPtr->flat() * mDiff; |
---|
326 | |
---|
327 | // Translate into relative momenta and find phase-space weight. |
---|
328 | p1Abs = 0.5 * sqrtpos( (m0 - m1t - m23) * (m0 + m1t + m23) |
---|
329 | * (m0 + m1t - m23) * (m0 - m1t + m23) ) / m0; |
---|
330 | p23Abs = 0.5 * sqrtpos( (m23 - m2t - m3t) * (m23 + m2t + m3t) |
---|
331 | * (m23 + m2t - m3t) * (m23 - m2t + m3t) ) / m23; |
---|
332 | wtPS = p1Abs * p23Abs; |
---|
333 | |
---|
334 | // If rejected, try again with new invariant masses. |
---|
335 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); |
---|
336 | |
---|
337 | // Set up m23 -> m2 + m3 isotropic in its rest frame. |
---|
338 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
339 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
340 | double phi23 = 2. * M_PI * rndmPtr->flat(); |
---|
341 | double pX = p23Abs * sinTheta * cos(phi23); |
---|
342 | double pY = p23Abs * sinTheta * sin(phi23); |
---|
343 | double pZ = p23Abs * cosTheta; |
---|
344 | double e2 = sqrt( m2t*m2t + p23Abs*p23Abs); |
---|
345 | double e3 = sqrt( m3t*m3t + p23Abs*p23Abs); |
---|
346 | Vec4 p2( pX, pY, pZ, e2); |
---|
347 | Vec4 p3( -pX, -pY, -pZ, e3); |
---|
348 | |
---|
349 | // Set up 0 -> 1 + 23 isotropic in its rest frame. |
---|
350 | cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
351 | sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
352 | phi23 = 2. * M_PI * rndmPtr->flat(); |
---|
353 | pX = p1Abs * sinTheta * cos(phi23); |
---|
354 | pY = p1Abs * sinTheta * sin(phi23); |
---|
355 | pZ = p1Abs * cosTheta; |
---|
356 | double e1 = sqrt( m1t*m1t + p1Abs*p1Abs); |
---|
357 | double e23 = sqrt( m23*m23 + p1Abs*p1Abs); |
---|
358 | Vec4 p1( pX, pY, pZ, e1); |
---|
359 | |
---|
360 | // Boost 2 + 3 to the 0 rest frame and then boost to lab frame. |
---|
361 | Vec4 p23( -pX, -pY, -pZ, e23); |
---|
362 | p2.bst( p23 ); |
---|
363 | p3.bst( p23 ); |
---|
364 | p1.bst( pRes ); |
---|
365 | p2.bst( pRes ); |
---|
366 | p3.bst( pRes ); |
---|
367 | |
---|
368 | // Done for three-body decay. |
---|
369 | process[i1].p( p1 ); |
---|
370 | process[i2].p( p2 ); |
---|
371 | process[i3].p( p3 ); |
---|
372 | return; |
---|
373 | } |
---|
374 | |
---|
375 | // Do a multibody decay using the M-generator algorithm. |
---|
376 | |
---|
377 | // Set up masses and four-momenta in a vector, with mother in slot 0. |
---|
378 | vector<double> mProd; |
---|
379 | mProd.push_back( m0); |
---|
380 | for (int i = i1; i <= process[iRes].daughter2(); ++i) |
---|
381 | mProd.push_back( process[i].m() ); |
---|
382 | vector<Vec4> pProd; |
---|
383 | pProd.push_back( pRes); |
---|
384 | |
---|
385 | // Sum of daughter masses. |
---|
386 | double mSum = mProd[1]; |
---|
387 | for (int i = 2; i <= mult; ++i) mSum += mProd[i]; |
---|
388 | double mDiff = m0 - mSum; |
---|
389 | |
---|
390 | // Begin setup of intermediate invariant masses. |
---|
391 | vector<double> mInv; |
---|
392 | for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]); |
---|
393 | |
---|
394 | // Calculate the maximum weight in the decay. |
---|
395 | double wtPSmax = 1. / WTCORRECTION[mult]; |
---|
396 | double mMaxWT = mDiff + mProd[mult]; |
---|
397 | double mMinWT = 0.; |
---|
398 | for (int i = mult - 1; i > 0; --i) { |
---|
399 | mMaxWT += mProd[i]; |
---|
400 | mMinWT += mProd[i+1]; |
---|
401 | double mNow = mProd[i]; |
---|
402 | wtPSmax *= 0.5 * sqrtpos( (mMaxWT - mMinWT - mNow) |
---|
403 | * (mMaxWT + mMinWT + mNow) * (mMaxWT + mMinWT - mNow) |
---|
404 | * (mMaxWT - mMinWT + mNow) ) / mMaxWT; |
---|
405 | } |
---|
406 | |
---|
407 | // Begin loop to find the set of intermediate invariant masses. |
---|
408 | vector<double> rndmOrd; |
---|
409 | double wtPS; |
---|
410 | do { |
---|
411 | wtPS = 1.; |
---|
412 | |
---|
413 | // Find and order random numbers in descending order. |
---|
414 | rndmOrd.resize(0); |
---|
415 | rndmOrd.push_back(1.); |
---|
416 | for (int i = 1; i < mult - 1; ++i) { |
---|
417 | double rndm = rndmPtr->flat(); |
---|
418 | rndmOrd.push_back(rndm); |
---|
419 | for (int j = i - 1; j > 0; --j) { |
---|
420 | if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] ); |
---|
421 | else break; |
---|
422 | } |
---|
423 | } |
---|
424 | rndmOrd.push_back(0.); |
---|
425 | |
---|
426 | // Translate into intermediate masses and find weight. |
---|
427 | for (int i = mult - 1; i > 0; --i) { |
---|
428 | mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff; |
---|
429 | wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i]) |
---|
430 | * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i]) |
---|
431 | * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i]; |
---|
432 | } |
---|
433 | |
---|
434 | // If rejected, try again with new invariant masses. |
---|
435 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); |
---|
436 | |
---|
437 | // Perform two-particle decays in the respective rest frame. |
---|
438 | vector<Vec4> pInv; |
---|
439 | pInv.resize(mult + 1); |
---|
440 | for (int i = 1; i < mult; ++i) { |
---|
441 | double p12 = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i]) |
---|
442 | * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i]) |
---|
443 | * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i]; |
---|
444 | |
---|
445 | // Isotropic angles give three-momentum. |
---|
446 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
447 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
448 | double phiLoc = 2. * M_PI * rndmPtr->flat(); |
---|
449 | double pX = p12 * sinTheta * cos(phiLoc); |
---|
450 | double pY = p12 * sinTheta * sin(phiLoc); |
---|
451 | double pZ = p12 * cosTheta; |
---|
452 | |
---|
453 | // Calculate energies, fill four-momenta. |
---|
454 | double eHad = sqrt( mProd[i]*mProd[i] + p12*p12); |
---|
455 | double eInv = sqrt( mInv[i+1]*mInv[i+1] + p12*p12); |
---|
456 | pProd.push_back( Vec4( pX, pY, pZ, eHad) ); |
---|
457 | pInv[i+1].p( -pX, -pY, -pZ, eInv); |
---|
458 | } |
---|
459 | pProd.push_back( pInv[mult] ); |
---|
460 | |
---|
461 | // Boost decay products to the mother rest frame and on to lab frame. |
---|
462 | pInv[1] = pProd[0]; |
---|
463 | for (int iFrame = mult - 1; iFrame > 0; --iFrame) |
---|
464 | for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]); |
---|
465 | |
---|
466 | // Done for multibody decay. |
---|
467 | for (int i = 1; i < mult; ++i) |
---|
468 | process[i1 + i - 1].p( pProd[i] ); |
---|
469 | return; |
---|
470 | |
---|
471 | } |
---|
472 | |
---|
473 | //-------------------------------------------------------------------------- |
---|
474 | |
---|
475 | // Determine how 3-body phase space should be sampled. |
---|
476 | |
---|
477 | void PhaseSpace::setup3Body() { |
---|
478 | |
---|
479 | // Check for massive t-channel propagator particles. |
---|
480 | int idTchan1 = abs( sigmaProcessPtr->idTchan1() ); |
---|
481 | int idTchan2 = abs( sigmaProcessPtr->idTchan2() ); |
---|
482 | mTchan1 = (idTchan1 == 0) ? pTHatMinDiverge |
---|
483 | : particleDataPtr->m0(idTchan1); |
---|
484 | mTchan2 = (idTchan2 == 0) ? pTHatMinDiverge |
---|
485 | : particleDataPtr->m0(idTchan2); |
---|
486 | sTchan1 = mTchan1 * mTchan1; |
---|
487 | sTchan2 = mTchan2 * mTchan2; |
---|
488 | |
---|
489 | // Find coefficients of different pT2 selection terms. Mirror choice. |
---|
490 | frac3Pow1 = sigmaProcessPtr->tChanFracPow1(); |
---|
491 | frac3Pow2 = sigmaProcessPtr->tChanFracPow2(); |
---|
492 | frac3Flat = 1. - frac3Pow1 - frac3Pow2; |
---|
493 | useMirrorWeight = sigmaProcessPtr->useMirrorWeight(); |
---|
494 | |
---|
495 | } |
---|
496 | |
---|
497 | //-------------------------------------------------------------------------- |
---|
498 | |
---|
499 | // Determine how phase space should be sampled. |
---|
500 | |
---|
501 | bool PhaseSpace::setupSampling123(bool is2, bool is3, ostream& os) { |
---|
502 | |
---|
503 | // Optional printout. |
---|
504 | if (showSearch) os << "\n PYTHIA Optimization printout for " |
---|
505 | << sigmaProcessPtr->name() << "\n \n" << scientific << setprecision(3); |
---|
506 | |
---|
507 | // Check that open range in tau (+ set tauMin, tauMax). |
---|
508 | if (!limitTau(is2, is3)) return false; |
---|
509 | |
---|
510 | // Reset coefficients and matrices of equation system to solve. |
---|
511 | int binTau[8], binY[8], binZ[8]; |
---|
512 | double vecTau[8], matTau[8][8], vecY[8], matY[8][8], vecZ[8], matZ[8][8]; |
---|
513 | for (int i = 0; i < 8; ++i) { |
---|
514 | tauCoef[i] = 0.; |
---|
515 | yCoef[i] = 0.; |
---|
516 | zCoef[i] = 0.; |
---|
517 | binTau[i] = 0; |
---|
518 | binY[i] = 0; |
---|
519 | binZ[i] = 0; |
---|
520 | vecTau[i] = 0.; |
---|
521 | vecY[i] = 0.; |
---|
522 | vecZ[i] = 0.; |
---|
523 | for (int j = 0; j < 8; ++j) { |
---|
524 | matTau[i][j] = 0.; |
---|
525 | matY[i][j] = 0.; |
---|
526 | matZ[i][j] = 0.; |
---|
527 | } |
---|
528 | } |
---|
529 | sigmaMx = 0.; |
---|
530 | sigmaNeg = 0.; |
---|
531 | |
---|
532 | // Number of used coefficients/points for each dimension: tau, y, c. |
---|
533 | nTau = (hasPointLeptons) ? 1 : 2; |
---|
534 | nY = (hasPointLeptons) ? 1 : 5; |
---|
535 | nZ = (is2) ? 5 : 1; |
---|
536 | |
---|
537 | // Identify if any resonances contribute in s-channel. |
---|
538 | idResA = sigmaProcessPtr->resonanceA(); |
---|
539 | if (idResA != 0) { |
---|
540 | mResA = particleDataPtr->m0(idResA); |
---|
541 | GammaResA = particleDataPtr->mWidth(idResA); |
---|
542 | if (mHatMin > mResA + WIDTHMARGIN * GammaResA || (mHatMax > 0. |
---|
543 | && mHatMax < mResA - WIDTHMARGIN * GammaResA) ) idResA = 0; |
---|
544 | } |
---|
545 | idResB = sigmaProcessPtr->resonanceB(); |
---|
546 | if (idResB != 0) { |
---|
547 | mResB = particleDataPtr->m0(idResB); |
---|
548 | GammaResB = particleDataPtr->mWidth(idResB); |
---|
549 | if (mHatMin > mResB + WIDTHMARGIN * GammaResB || (mHatMax > 0. |
---|
550 | && mHatMax < mResB - WIDTHMARGIN * GammaResB) ) idResB = 0; |
---|
551 | } |
---|
552 | if (idResA == 0 && idResB != 0) { |
---|
553 | idResA = idResB; |
---|
554 | mResA = mResB; |
---|
555 | GammaResA = GammaResB; |
---|
556 | idResB = 0; |
---|
557 | } |
---|
558 | |
---|
559 | // More sampling in tau if resonances in s-channel. |
---|
560 | if (idResA !=0 && !hasPointLeptons) { |
---|
561 | nTau += 2; |
---|
562 | tauResA = mResA * mResA / s; |
---|
563 | widResA = mResA * GammaResA / s; |
---|
564 | } |
---|
565 | if (idResB != 0 && !hasPointLeptons) { |
---|
566 | nTau += 2; |
---|
567 | tauResB = mResB * mResB / s; |
---|
568 | widResB = mResB * GammaResB / s; |
---|
569 | } |
---|
570 | |
---|
571 | // More sampling in tau (and different in y) if incoming lepton beams. |
---|
572 | if (hasLeptonBeams && !hasPointLeptons) ++nTau; |
---|
573 | |
---|
574 | // Special case when both resonances have same mass. |
---|
575 | sameResMass = false; |
---|
576 | if (idResB != 0 && abs(mResA - mResB) < SAMEMASS * (GammaResA + GammaResB)) |
---|
577 | sameResMass = true; |
---|
578 | |
---|
579 | // Default z value and weight required for 2 -> 1. Number of dimensions. |
---|
580 | z = 0.; |
---|
581 | wtZ = 1.; |
---|
582 | int nVar = (is2) ? 3 : 2; |
---|
583 | |
---|
584 | // Initial values, to be modified later. |
---|
585 | tauCoef[0] = 1.; |
---|
586 | yCoef[1] = 0.5; |
---|
587 | yCoef[2] = 0.5; |
---|
588 | zCoef[0] = 1.; |
---|
589 | |
---|
590 | // Step through grid in tau. Set limits on y and z generation. |
---|
591 | for (int iTau = 0; iTau < nTau; ++iTau) { |
---|
592 | double posTau = 0.5; |
---|
593 | if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6; |
---|
594 | selectTau( iTau, posTau, is2); |
---|
595 | if (!limitY()) continue; |
---|
596 | if (is2 && !limitZ()) continue; |
---|
597 | |
---|
598 | // Step through grids in y and z. |
---|
599 | for (int iY = 0; iY < nY; ++iY) { |
---|
600 | selectY( iY, 0.5); |
---|
601 | for (int iZ = 0; iZ < nZ; ++iZ) { |
---|
602 | if (is2) selectZ( iZ, 0.5); |
---|
603 | double sigmaTmp = 0.; |
---|
604 | |
---|
605 | // 2 -> 1: calculate cross section, weighted by phase-space volume. |
---|
606 | if (!is2 && !is3) { |
---|
607 | sigmaProcessPtr->set1Kin( x1H, x2H, sH); |
---|
608 | sigmaTmp = sigmaProcessPtr->sigmaPDF(); |
---|
609 | sigmaTmp *= wtTau * wtY; |
---|
610 | |
---|
611 | // 2 -> 2: calculate cross section, weighted by phase-space volume |
---|
612 | // and Breit-Wigners for masses |
---|
613 | } else if (is2) { |
---|
614 | sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, |
---|
615 | runBW3H, runBW4H); |
---|
616 | sigmaTmp = sigmaProcessPtr->sigmaPDF(); |
---|
617 | sigmaTmp *= wtTau * wtY * wtZ * wtBW; |
---|
618 | |
---|
619 | // 2 -> 3: repeat internal 3-body phase space several times and |
---|
620 | // keep maximal cross section, weighted by phase-space volume |
---|
621 | // and Breit-Wigners for masses |
---|
622 | } else if (is3) { |
---|
623 | for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) { |
---|
624 | if (!select3Body()) continue; |
---|
625 | sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm, |
---|
626 | m3, m4, m5, runBW3H, runBW4H, runBW5H); |
---|
627 | double sigmaTry = sigmaProcessPtr->sigmaPDF(); |
---|
628 | sigmaTry *= wtTau * wtY * wt3Body * wtBW; |
---|
629 | if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry; |
---|
630 | } |
---|
631 | } |
---|
632 | |
---|
633 | // Allow possibility for user to modify cross section. (3body??) |
---|
634 | if (canModifySigma) sigmaTmp |
---|
635 | *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false); |
---|
636 | if (canBiasSelection) sigmaTmp |
---|
637 | *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, false); |
---|
638 | if (canBias2Sel) sigmaTmp *= pow( pTH / bias2SelRef, bias2SelPow); |
---|
639 | |
---|
640 | // Check if current maximum exceeded. |
---|
641 | if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp; |
---|
642 | |
---|
643 | // Optional printout. Protect against negative cross sections. |
---|
644 | if (showSearch) os << " tau =" << setw(11) << tau << " y =" |
---|
645 | << setw(11) << y << " z =" << setw(11) << z |
---|
646 | << " sigma =" << setw(11) << sigmaTmp << "\n"; |
---|
647 | if (sigmaTmp < 0.) sigmaTmp = 0.; |
---|
648 | |
---|
649 | // Sum up tau cross-section pieces in points used. |
---|
650 | if (!hasPointLeptons) { |
---|
651 | binTau[iTau] += 1; |
---|
652 | vecTau[iTau] += sigmaTmp; |
---|
653 | matTau[iTau][0] += 1. / intTau0; |
---|
654 | matTau[iTau][1] += (1. / intTau1) / tau; |
---|
655 | if (idResA != 0) { |
---|
656 | matTau[iTau][2] += (1. / intTau2) / (tau + tauResA); |
---|
657 | matTau[iTau][3] += (1. / intTau3) |
---|
658 | * tau / ( pow2(tau - tauResA) + pow2(widResA) ); |
---|
659 | } |
---|
660 | if (idResB != 0) { |
---|
661 | matTau[iTau][4] += (1. / intTau4) / (tau + tauResB); |
---|
662 | matTau[iTau][5] += (1. / intTau5) |
---|
663 | * tau / ( pow2(tau - tauResB) + pow2(widResB) ); |
---|
664 | } |
---|
665 | if (hasLeptonBeams) matTau[iTau][nTau - 1] += (1. / intTau6) |
---|
666 | * tau / max( LEPTONTAUMIN, 1. - tau); |
---|
667 | } |
---|
668 | |
---|
669 | // Sum up y cross-section pieces in points used. |
---|
670 | if (!hasPointLeptons) { |
---|
671 | binY[iY] += 1; |
---|
672 | vecY[iY] += sigmaTmp; |
---|
673 | matY[iY][0] += (yMax / intY0) / cosh(y); |
---|
674 | matY[iY][1] += (yMax / intY12) * (y + yMax); |
---|
675 | matY[iY][2] += (yMax / intY12) * (yMax - y); |
---|
676 | if (!hasLeptonBeams) { |
---|
677 | matY[iY][3] += (yMax / intY34) * exp(y); |
---|
678 | matY[iY][4] += (yMax / intY34) * exp(-y); |
---|
679 | } else { |
---|
680 | matY[iY][3] += (yMax / intY56) |
---|
681 | / max( LEPTONXMIN, 1. - exp( y - yMax) ); |
---|
682 | matY[iY][4] += (yMax / intY56) |
---|
683 | / max( LEPTONXMIN, 1. - exp(-y - yMax) ); |
---|
684 | } |
---|
685 | } |
---|
686 | |
---|
687 | // Integrals over z expressions at tauMax, to be used below. |
---|
688 | if (is2) { |
---|
689 | double p2AbsMax = 0.25 * (pow2(tauMax * s - s3 - s4) |
---|
690 | - 4. * s3 * s4) / (tauMax * s); |
---|
691 | double zMaxMax = sqrtpos( 1. - pT2HatMin / p2AbsMax ); |
---|
692 | double zPosMaxMax = max(ratio34, unity34 + zMaxMax); |
---|
693 | double zNegMaxMax = max(ratio34, unity34 - zMaxMax); |
---|
694 | double intZ0Max = 2. * zMaxMax; |
---|
695 | double intZ12Max = log( zPosMaxMax / zNegMaxMax); |
---|
696 | double intZ34Max = 1. / zNegMaxMax - 1. / zPosMaxMax; |
---|
697 | |
---|
698 | // Sum up z cross-section pieces in points used. |
---|
699 | binZ[iZ] += 1; |
---|
700 | vecZ[iZ] += sigmaTmp; |
---|
701 | matZ[iZ][0] += 1.; |
---|
702 | matZ[iZ][1] += (intZ0Max / intZ12Max) / zNeg; |
---|
703 | matZ[iZ][2] += (intZ0Max / intZ12Max) / zPos; |
---|
704 | matZ[iZ][3] += (intZ0Max / intZ34Max) / pow2(zNeg); |
---|
705 | matZ[iZ][4] += (intZ0Max / intZ34Max) / pow2(zPos); |
---|
706 | } |
---|
707 | |
---|
708 | // End of loops over phase space points. |
---|
709 | } |
---|
710 | } |
---|
711 | } |
---|
712 | |
---|
713 | // Fail if no non-vanishing cross sections. |
---|
714 | if (sigmaMx <= 0.) { |
---|
715 | sigmaMx = 0.; |
---|
716 | return false; |
---|
717 | } |
---|
718 | |
---|
719 | // Solve respective equation system for better phase space coefficients. |
---|
720 | if (!hasPointLeptons) solveSys( nTau, binTau, vecTau, matTau, tauCoef); |
---|
721 | if (!hasPointLeptons) solveSys( nY, binY, vecY, matY, yCoef); |
---|
722 | if (is2) solveSys( nZ, binZ, vecZ, matZ, zCoef); |
---|
723 | if (showSearch) os << "\n"; |
---|
724 | |
---|
725 | // Provide cumulative sum of coefficients. |
---|
726 | tauCoefSum[0] = tauCoef[0]; |
---|
727 | yCoefSum[0] = yCoef[0]; |
---|
728 | zCoefSum[0] = zCoef[0]; |
---|
729 | for (int i = 1; i < 8; ++ i) { |
---|
730 | tauCoefSum[i] = tauCoefSum[i - 1] + tauCoef[i]; |
---|
731 | yCoefSum[i] = yCoefSum[i - 1] + yCoef[i]; |
---|
732 | zCoefSum[i] = zCoefSum[i - 1] + zCoef[i]; |
---|
733 | } |
---|
734 | // The last element should be > 1 to be on safe side in selection below. |
---|
735 | tauCoefSum[nTau - 1] = 2.; |
---|
736 | yCoefSum[nY - 1] = 2.; |
---|
737 | zCoefSum[nZ - 1] = 2.; |
---|
738 | |
---|
739 | |
---|
740 | // Begin find two most promising maxima among same points as before. |
---|
741 | int iMaxTau[NMAXTRY + 2], iMaxY[NMAXTRY + 2], iMaxZ[NMAXTRY + 2]; |
---|
742 | double sigMax[NMAXTRY + 2]; |
---|
743 | int nMax = 0; |
---|
744 | |
---|
745 | // Scan same grid as before in tau, y, z. |
---|
746 | for (int iTau = 0; iTau < nTau; ++iTau) { |
---|
747 | double posTau = 0.5; |
---|
748 | if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6; |
---|
749 | selectTau( iTau, posTau, is2); |
---|
750 | if (!limitY()) continue; |
---|
751 | if (is2 && !limitZ()) continue; |
---|
752 | for (int iY = 0; iY < nY; ++iY) { |
---|
753 | selectY( iY, 0.5); |
---|
754 | for (int iZ = 0; iZ < nZ; ++iZ) { |
---|
755 | if (is2) selectZ( iZ, 0.5); |
---|
756 | double sigmaTmp = 0.; |
---|
757 | |
---|
758 | // 2 -> 1: calculate cross section, weighted by phase-space volume. |
---|
759 | if (!is2 && !is3) { |
---|
760 | sigmaProcessPtr->set1Kin( x1H, x2H, sH); |
---|
761 | sigmaTmp = sigmaProcessPtr->sigmaPDF(); |
---|
762 | sigmaTmp *= wtTau * wtY; |
---|
763 | |
---|
764 | // 2 -> 2: calculate cross section, weighted by phase-space volume |
---|
765 | // and Breit-Wigners for masses |
---|
766 | } else if (is2) { |
---|
767 | sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, |
---|
768 | runBW3H, runBW4H); |
---|
769 | sigmaTmp = sigmaProcessPtr->sigmaPDF(); |
---|
770 | sigmaTmp *= wtTau * wtY * wtZ * wtBW; |
---|
771 | |
---|
772 | // 2 -> 3: repeat internal 3-body phase space several times and |
---|
773 | // keep maximal cross section, weighted by phase-space volume |
---|
774 | // and Breit-Wigners for masses |
---|
775 | } else if (is3) { |
---|
776 | for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) { |
---|
777 | if (!select3Body()) continue; |
---|
778 | sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm, |
---|
779 | m3, m4, m5, runBW3H, runBW4H, runBW5H); |
---|
780 | double sigmaTry = sigmaProcessPtr->sigmaPDF(); |
---|
781 | sigmaTry *= wtTau * wtY * wt3Body * wtBW; |
---|
782 | if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry; |
---|
783 | } |
---|
784 | } |
---|
785 | |
---|
786 | // Allow possibility for user to modify cross section. (3body??) |
---|
787 | if (canModifySigma) sigmaTmp |
---|
788 | *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false); |
---|
789 | if (canBiasSelection) sigmaTmp |
---|
790 | *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, false); |
---|
791 | if (canBias2Sel) sigmaTmp *= pow( pTH / bias2SelRef, bias2SelPow); |
---|
792 | |
---|
793 | // Optional printout. Protect against negative cross section. |
---|
794 | if (showSearch) os << " tau =" << setw(11) << tau << " y =" |
---|
795 | << setw(11) << y << " z =" << setw(11) << z |
---|
796 | << " sigma =" << setw(11) << sigmaTmp << "\n"; |
---|
797 | if (sigmaTmp < 0.) sigmaTmp = 0.; |
---|
798 | |
---|
799 | // Check that point is not simply mirror of already found one. |
---|
800 | bool mirrorPoint = false; |
---|
801 | for (int iMove = 0; iMove < nMax; ++iMove) |
---|
802 | if (abs(sigmaTmp - sigMax[iMove]) < SAMESIGMA |
---|
803 | * (sigmaTmp + sigMax[iMove])) mirrorPoint = true; |
---|
804 | |
---|
805 | // Add to or insert in maximum list. Only first two count. |
---|
806 | if (!mirrorPoint) { |
---|
807 | int iInsert = 0; |
---|
808 | for (int iMove = nMax - 1; iMove >= -1; --iMove) { |
---|
809 | iInsert = iMove + 1; |
---|
810 | if (iInsert == 0 || sigmaTmp < sigMax[iMove]) break; |
---|
811 | iMaxTau[iMove + 1] = iMaxTau[iMove]; |
---|
812 | iMaxY[iMove + 1] = iMaxY[iMove]; |
---|
813 | iMaxZ[iMove + 1] = iMaxZ[iMove]; |
---|
814 | sigMax[iMove + 1] = sigMax[iMove]; |
---|
815 | } |
---|
816 | iMaxTau[iInsert] = iTau; |
---|
817 | iMaxY[iInsert] = iY; |
---|
818 | iMaxZ[iInsert] = iZ; |
---|
819 | sigMax[iInsert] = sigmaTmp; |
---|
820 | if (nMax < NMAXTRY) ++nMax; |
---|
821 | } |
---|
822 | |
---|
823 | // Found two most promising maxima. |
---|
824 | } |
---|
825 | } |
---|
826 | } |
---|
827 | if (showSearch) os << "\n"; |
---|
828 | |
---|
829 | // Read out starting position for search. |
---|
830 | sigmaMx = sigMax[0]; |
---|
831 | int beginVar = (hasPointLeptons) ? 2 : 0; |
---|
832 | for (int iMax = 0; iMax < nMax; ++iMax) { |
---|
833 | int iTau = iMaxTau[iMax]; |
---|
834 | int iY = iMaxY[iMax]; |
---|
835 | int iZ = iMaxZ[iMax]; |
---|
836 | double tauVal = 0.5; |
---|
837 | double yVal = 0.5; |
---|
838 | double zVal = 0.5; |
---|
839 | int iGrid; |
---|
840 | double varVal, varNew, deltaVar, marginVar, sigGrid[3]; |
---|
841 | |
---|
842 | // Starting point and step size in parameter space. |
---|
843 | for (int iRepeat = 0; iRepeat < 2; ++iRepeat) { |
---|
844 | // Run through (possibly a subset of) tau, y and z. |
---|
845 | for (int iVar = beginVar; iVar < nVar; ++iVar) { |
---|
846 | if (iVar == 0) varVal = tauVal; |
---|
847 | else if (iVar == 1) varVal = yVal; |
---|
848 | else varVal = zVal; |
---|
849 | deltaVar = (iRepeat == 0) ? 0.1 |
---|
850 | : max( 0.01, min( 0.05, min( varVal - 0.02, 0.98 - varVal) ) ); |
---|
851 | marginVar = (iRepeat == 0) ? 0.02 : 0.002; |
---|
852 | int moveStart = (iRepeat == 0 && iVar == 0) ? 0 : 1; |
---|
853 | for (int move = moveStart; move < 9; ++move) { |
---|
854 | |
---|
855 | // Define new parameter-space point by step in one dimension. |
---|
856 | if (move == 0) { |
---|
857 | iGrid = 1; |
---|
858 | varNew = varVal; |
---|
859 | } else if (move == 1) { |
---|
860 | iGrid = 2; |
---|
861 | varNew = varVal + deltaVar; |
---|
862 | } else if (move == 2) { |
---|
863 | iGrid = 0; |
---|
864 | varNew = varVal - deltaVar; |
---|
865 | } else if (sigGrid[2] >= max( sigGrid[0], sigGrid[1]) |
---|
866 | && varVal + 2. * deltaVar < 1. - marginVar) { |
---|
867 | varVal += deltaVar; |
---|
868 | sigGrid[0] = sigGrid[1]; |
---|
869 | sigGrid[1] = sigGrid[2]; |
---|
870 | iGrid = 2; |
---|
871 | varNew = varVal + deltaVar; |
---|
872 | } else if (sigGrid[0] >= max( sigGrid[1], sigGrid[2]) |
---|
873 | && varVal - 2. * deltaVar > marginVar) { |
---|
874 | varVal -= deltaVar; |
---|
875 | sigGrid[2] = sigGrid[1]; |
---|
876 | sigGrid[1] = sigGrid[0]; |
---|
877 | iGrid = 0; |
---|
878 | varNew = varVal - deltaVar; |
---|
879 | } else if (sigGrid[2] >= sigGrid[0]) { |
---|
880 | deltaVar *= 0.5; |
---|
881 | varVal += deltaVar; |
---|
882 | sigGrid[0] = sigGrid[1]; |
---|
883 | iGrid = 1; |
---|
884 | varNew = varVal; |
---|
885 | } else { |
---|
886 | deltaVar *= 0.5; |
---|
887 | varVal -= deltaVar; |
---|
888 | sigGrid[2] = sigGrid[1]; |
---|
889 | iGrid = 1; |
---|
890 | varNew = varVal; |
---|
891 | } |
---|
892 | |
---|
893 | // Convert to relevant variables and find derived new limits. |
---|
894 | bool insideLimits = true; |
---|
895 | if (iVar == 0) { |
---|
896 | tauVal = varNew; |
---|
897 | selectTau( iTau, tauVal, is2); |
---|
898 | if (!limitY()) insideLimits = false; |
---|
899 | if (is2 && !limitZ()) insideLimits = false; |
---|
900 | if (insideLimits) { |
---|
901 | selectY( iY, yVal); |
---|
902 | if (is2) selectZ( iZ, zVal); |
---|
903 | } |
---|
904 | } else if (iVar == 1) { |
---|
905 | yVal = varNew; |
---|
906 | selectY( iY, yVal); |
---|
907 | } else if (iVar == 2) { |
---|
908 | zVal = varNew; |
---|
909 | selectZ( iZ, zVal); |
---|
910 | } |
---|
911 | |
---|
912 | // Evaluate cross-section. |
---|
913 | double sigmaTmp = 0.; |
---|
914 | if (insideLimits) { |
---|
915 | |
---|
916 | // 2 -> 1: calculate cross section, weighted by phase-space volume. |
---|
917 | if (!is2 && !is3) { |
---|
918 | sigmaProcessPtr->set1Kin( x1H, x2H, sH); |
---|
919 | sigmaTmp = sigmaProcessPtr->sigmaPDF(); |
---|
920 | sigmaTmp *= wtTau * wtY; |
---|
921 | |
---|
922 | // 2 -> 2: calculate cross section, weighted by phase-space volume |
---|
923 | // and Breit-Wigners for masses |
---|
924 | } else if (is2) { |
---|
925 | sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, |
---|
926 | runBW3H, runBW4H); |
---|
927 | sigmaTmp = sigmaProcessPtr->sigmaPDF(); |
---|
928 | sigmaTmp *= wtTau * wtY * wtZ * wtBW; |
---|
929 | |
---|
930 | // 2 -> 3: repeat internal 3-body phase space several times and |
---|
931 | // keep maximal cross section, weighted by phase-space volume |
---|
932 | // and Breit-Wigners for masses |
---|
933 | } else if (is3) { |
---|
934 | for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) { |
---|
935 | if (!select3Body()) continue; |
---|
936 | sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm, |
---|
937 | m3, m4, m5, runBW3H, runBW4H, runBW5H); |
---|
938 | double sigmaTry = sigmaProcessPtr->sigmaPDF(); |
---|
939 | sigmaTry *= wtTau * wtY * wt3Body * wtBW; |
---|
940 | if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry; |
---|
941 | } |
---|
942 | } |
---|
943 | |
---|
944 | // Allow possibility for user to modify cross section. |
---|
945 | if (canModifySigma) sigmaTmp |
---|
946 | *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false); |
---|
947 | if (canBiasSelection) sigmaTmp |
---|
948 | *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, false); |
---|
949 | if (canBias2Sel) sigmaTmp *= pow( pTH / bias2SelRef, bias2SelPow); |
---|
950 | |
---|
951 | // Optional printout. Protect against negative cross section. |
---|
952 | if (showSearch) os << " tau =" << setw(11) << tau << " y =" |
---|
953 | << setw(11) << y << " z =" << setw(11) << z |
---|
954 | << " sigma =" << setw(11) << sigmaTmp << "\n"; |
---|
955 | if (sigmaTmp < 0.) sigmaTmp = 0.; |
---|
956 | } |
---|
957 | |
---|
958 | // Save new maximum. Final maximum. |
---|
959 | sigGrid[iGrid] = sigmaTmp; |
---|
960 | if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp; |
---|
961 | } |
---|
962 | } |
---|
963 | } |
---|
964 | } |
---|
965 | sigmaMx *= SAFETYMARGIN; |
---|
966 | sigmaPos = sigmaMx; |
---|
967 | |
---|
968 | // Optional printout. |
---|
969 | if (showSearch) os << "\n Final maximum = " << setw(11) << sigmaMx << endl; |
---|
970 | |
---|
971 | // Done. |
---|
972 | return true; |
---|
973 | } |
---|
974 | |
---|
975 | //-------------------------------------------------------------------------- |
---|
976 | |
---|
977 | // Select a trial kinematics phase space point. |
---|
978 | // Note: by In is meant the integral over the quantity multiplying |
---|
979 | // coefficient cn. The sum of cn is normalized to unity. |
---|
980 | |
---|
981 | bool PhaseSpace::trialKin123(bool is2, bool is3, bool inEvent, ostream& os) { |
---|
982 | |
---|
983 | // Allow for possibility that energy varies from event to event. |
---|
984 | if (doEnergySpread) { |
---|
985 | eCM = infoPtr->eCM(); |
---|
986 | s = eCM * eCM; |
---|
987 | |
---|
988 | // Find shifted tauRes values. |
---|
989 | if (idResA !=0 && !hasPointLeptons) { |
---|
990 | tauResA = mResA * mResA / s; |
---|
991 | widResA = mResA * GammaResA / s; |
---|
992 | } |
---|
993 | if (idResB != 0 && !hasPointLeptons) { |
---|
994 | tauResB = mResB * mResB / s; |
---|
995 | widResB = mResB * GammaResB / s; |
---|
996 | } |
---|
997 | } |
---|
998 | |
---|
999 | // Choose tau according to h1(tau)/tau, where |
---|
1000 | // h1(tau) = c0/I0 + (c1/I1) * 1/tau |
---|
1001 | // + (c2/I2) / (tau + tauResA) |
---|
1002 | // + (c3/I3) * tau / ((tau - tauResA)^2 + widResA^2) |
---|
1003 | // + (c4/I4) / (tau + tauResB) |
---|
1004 | // + (c5/I5) * tau / ((tau - tauResB)^2 + widResB^2) |
---|
1005 | // + (c6/I6) * tau / (1 - tau). |
---|
1006 | if (!limitTau(is2, is3)) return false; |
---|
1007 | int iTau = 0; |
---|
1008 | if (!hasPointLeptons) { |
---|
1009 | double rTau = rndmPtr->flat(); |
---|
1010 | while (rTau > tauCoefSum[iTau]) ++iTau; |
---|
1011 | } |
---|
1012 | selectTau( iTau, rndmPtr->flat(), is2); |
---|
1013 | |
---|
1014 | // Choose y according to h2(y), where |
---|
1015 | // h2(y) = (c0/I0) * 1/cosh(y) |
---|
1016 | // + (c1/I1) * (y-ymin) + (c2/I2) * (ymax-y) |
---|
1017 | // + (c3/I3) * exp(y) + (c4/i4) * exp(-y) (for hadron; for lepton instead) |
---|
1018 | // + (c5/I5) * 1 / (1 - exp(y-ymax)) + (c6/I6) * 1 / (1 - exp(ymin-y)). |
---|
1019 | if (!limitY()) return false; |
---|
1020 | int iY = 0; |
---|
1021 | if (!hasPointLeptons) { |
---|
1022 | double rY = rndmPtr->flat(); |
---|
1023 | while (rY > yCoefSum[iY]) ++iY; |
---|
1024 | } |
---|
1025 | selectY( iY, rndmPtr->flat()); |
---|
1026 | |
---|
1027 | // Choose z = cos(thetaHat) according to h3(z), where |
---|
1028 | // h3(z) = c0/I0 + (c1/I1) * 1/(A - z) + (c2/I2) * 1/(A + z) |
---|
1029 | // + (c3/I3) * 1/(A - z)^2 + (c4/I4) * 1/(A + z)^2, |
---|
1030 | // where A = 1 + 2*(m3*m4/sH)^2 (= 1 for massless products). |
---|
1031 | if (is2) { |
---|
1032 | if (!limitZ()) return false; |
---|
1033 | int iZ = 0; |
---|
1034 | double rZ = rndmPtr->flat(); |
---|
1035 | while (rZ > zCoefSum[iZ]) ++iZ; |
---|
1036 | selectZ( iZ, rndmPtr->flat()); |
---|
1037 | } |
---|
1038 | |
---|
1039 | // 2 -> 1: calculate cross section, weighted by phase-space volume. |
---|
1040 | if (!is2 && !is3) { |
---|
1041 | sigmaProcessPtr->set1Kin( x1H, x2H, sH); |
---|
1042 | sigmaNw = sigmaProcessPtr->sigmaPDF(); |
---|
1043 | sigmaNw *= wtTau * wtY; |
---|
1044 | |
---|
1045 | // 2 -> 2: calculate cross section, weighted by phase-space volume |
---|
1046 | // and Breit-Wigners for masses |
---|
1047 | } else if (is2) { |
---|
1048 | sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, runBW3H, runBW4H); |
---|
1049 | sigmaNw = sigmaProcessPtr->sigmaPDF(); |
---|
1050 | sigmaNw *= wtTau * wtY * wtZ * wtBW; |
---|
1051 | |
---|
1052 | // 2 -> 3: also sample internal 3-body phase, weighted by |
---|
1053 | // 2 -> 1 phase-space volume and Breit-Wigners for masses |
---|
1054 | } else if (is3) { |
---|
1055 | if (!select3Body()) sigmaNw = 0.; |
---|
1056 | else { |
---|
1057 | sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm, |
---|
1058 | m3, m4, m5, runBW3H, runBW4H, runBW5H); |
---|
1059 | sigmaNw = sigmaProcessPtr->sigmaPDF(); |
---|
1060 | sigmaNw *= wtTau * wtY * wt3Body * wtBW; |
---|
1061 | } |
---|
1062 | } |
---|
1063 | |
---|
1064 | // Allow possibility for user to modify cross section. |
---|
1065 | if (canModifySigma) sigmaNw |
---|
1066 | *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, inEvent); |
---|
1067 | if (canBiasSelection) sigmaNw |
---|
1068 | *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, inEvent); |
---|
1069 | if (canBias2Sel) sigmaNw *= pow( pTH / bias2SelRef, bias2SelPow); |
---|
1070 | |
---|
1071 | // Check if maximum violated. |
---|
1072 | newSigmaMx = false; |
---|
1073 | if (sigmaNw > sigmaMx) { |
---|
1074 | infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin: " |
---|
1075 | "maximum for cross section violated"); |
---|
1076 | |
---|
1077 | // Violation strategy 1: increase maximum (always during initialization). |
---|
1078 | if (increaseMaximum || !inEvent) { |
---|
1079 | double violFact = SAFETYMARGIN * sigmaNw / sigmaMx; |
---|
1080 | sigmaMx = SAFETYMARGIN * sigmaNw; |
---|
1081 | newSigmaMx = true; |
---|
1082 | if (showViolation) { |
---|
1083 | if (violFact < 9.99) os << fixed; |
---|
1084 | else os << scientific; |
---|
1085 | os << " PYTHIA Maximum for " << sigmaProcessPtr->name() |
---|
1086 | << " increased by factor " << setprecision(3) << violFact |
---|
1087 | << " to " << scientific << sigmaMx << endl; |
---|
1088 | } |
---|
1089 | |
---|
1090 | // Violation strategy 2: weight event (done in ProcessContainer). |
---|
1091 | } else if (showViolation && sigmaNw > sigmaPos) { |
---|
1092 | double violFact = sigmaNw / sigmaMx; |
---|
1093 | if (violFact < 9.99) os << fixed; |
---|
1094 | else os << scientific; |
---|
1095 | os << " PYTHIA Maximum for " << sigmaProcessPtr->name() |
---|
1096 | << " exceeded by factor " << setprecision(3) << violFact << endl; |
---|
1097 | sigmaPos = sigmaNw; |
---|
1098 | } |
---|
1099 | } |
---|
1100 | |
---|
1101 | // Check if negative cross section. |
---|
1102 | if (sigmaNw < sigmaNeg) { |
---|
1103 | infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin:" |
---|
1104 | " negative cross section set 0", "for " + sigmaProcessPtr->name() ); |
---|
1105 | sigmaNeg = sigmaNw; |
---|
1106 | |
---|
1107 | // Optional printout of (all) violations. |
---|
1108 | if (showViolation) os << " PYTHIA Negative minimum for " |
---|
1109 | << sigmaProcessPtr->name() << " changed to " << scientific |
---|
1110 | << setprecision(3) << sigmaNeg << endl; |
---|
1111 | } |
---|
1112 | if (sigmaNw < 0.) sigmaNw = 0.; |
---|
1113 | |
---|
1114 | // Set event weight, where relevant. |
---|
1115 | biasWt = (canBiasSelection) ? userHooksPtr->biasedSelectionWeight() : 1.; |
---|
1116 | if (canBias2Sel) biasWt /= pow( pTH / bias2SelRef, bias2SelPow); |
---|
1117 | |
---|
1118 | // Done. |
---|
1119 | return true; |
---|
1120 | } |
---|
1121 | |
---|
1122 | //-------------------------------------------------------------------------- |
---|
1123 | |
---|
1124 | // Find range of allowed tau values. |
---|
1125 | |
---|
1126 | bool PhaseSpace::limitTau(bool is2, bool is3) { |
---|
1127 | |
---|
1128 | // Trivial reply for unresolved lepton beams. |
---|
1129 | if (hasPointLeptons) { |
---|
1130 | tauMin = 1.; |
---|
1131 | tauMax = 1.; |
---|
1132 | return true; |
---|
1133 | } |
---|
1134 | |
---|
1135 | // Requirements from allowed mHat range. |
---|
1136 | tauMin = sHatMin / s; |
---|
1137 | tauMax = (mHatMax < mHatMin) ? 1. : min( 1., sHatMax / s); |
---|
1138 | |
---|
1139 | // Requirements from allowed pT range and masses. |
---|
1140 | if (is2 || is3) { |
---|
1141 | double mT3Min = sqrt(s3 + pT2HatMin); |
---|
1142 | double mT4Min = sqrt(s4 + pT2HatMin); |
---|
1143 | double mT5Min = (is3) ? sqrt(s5 + pT2HatMin) : 0.; |
---|
1144 | tauMin = max( tauMin, pow2(mT3Min + mT4Min + mT5Min) / s); |
---|
1145 | } |
---|
1146 | |
---|
1147 | // Check that there is an open range. |
---|
1148 | return (tauMax > tauMin); |
---|
1149 | } |
---|
1150 | |
---|
1151 | //-------------------------------------------------------------------------- |
---|
1152 | |
---|
1153 | // Find range of allowed y values. |
---|
1154 | |
---|
1155 | bool PhaseSpace::limitY() { |
---|
1156 | |
---|
1157 | // Trivial reply for unresolved lepton beams. |
---|
1158 | if (hasPointLeptons) { |
---|
1159 | yMax = 1.; |
---|
1160 | return true; |
---|
1161 | } |
---|
1162 | |
---|
1163 | // Requirements from selected tau value. |
---|
1164 | yMax = -0.5 * log(tau); |
---|
1165 | |
---|
1166 | // For lepton beams requirements from cutoff for f_e^e. |
---|
1167 | double yMaxMargin = (hasLeptonBeams) ? yMax + LEPTONXLOGMAX : yMax; |
---|
1168 | |
---|
1169 | // Check that there is an open range. |
---|
1170 | return (yMaxMargin > 0.); |
---|
1171 | } |
---|
1172 | |
---|
1173 | //-------------------------------------------------------------------------- |
---|
1174 | |
---|
1175 | // Find range of allowed z = cos(theta) values. |
---|
1176 | |
---|
1177 | bool PhaseSpace::limitZ() { |
---|
1178 | |
---|
1179 | // Default limits. |
---|
1180 | zMin = 0.; |
---|
1181 | zMax = 1.; |
---|
1182 | |
---|
1183 | // Requirements from pTHat limits. |
---|
1184 | zMax = sqrtpos( 1. - pT2HatMin / p2Abs ); |
---|
1185 | if (pTHatMax > pTHatMin) zMin = sqrtpos( 1. - pT2HatMax / p2Abs ); |
---|
1186 | |
---|
1187 | // Check that there is an open range. |
---|
1188 | return (zMax > zMin); |
---|
1189 | } |
---|
1190 | |
---|
1191 | //-------------------------------------------------------------------------- |
---|
1192 | |
---|
1193 | // Select tau according to a choice of shapes. |
---|
1194 | |
---|
1195 | void PhaseSpace::selectTau(int iTau, double tauVal, bool is2) { |
---|
1196 | |
---|
1197 | // Trivial reply for unresolved lepton beams. |
---|
1198 | if (hasPointLeptons) { |
---|
1199 | tau = 1.; |
---|
1200 | wtTau = 1.; |
---|
1201 | sH = s; |
---|
1202 | mHat = sqrt(sH); |
---|
1203 | if (is2) { |
---|
1204 | p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH; |
---|
1205 | pAbs = sqrtpos( p2Abs ); |
---|
1206 | } |
---|
1207 | return; |
---|
1208 | } |
---|
1209 | |
---|
1210 | // Contributions from s-channel resonances. |
---|
1211 | double tRatA = 0.; |
---|
1212 | double aLowA = 0.; |
---|
1213 | double aUppA = 0.; |
---|
1214 | if (idResA !=0) { |
---|
1215 | tRatA = ((tauResA + tauMax) / (tauResA + tauMin)) * (tauMin / tauMax); |
---|
1216 | aLowA = atan( (tauMin - tauResA) / widResA); |
---|
1217 | aUppA = atan( (tauMax - tauResA) / widResA); |
---|
1218 | } |
---|
1219 | double tRatB = 0.; |
---|
1220 | double aLowB = 0.; |
---|
1221 | double aUppB = 0.; |
---|
1222 | if (idResB != 0) { |
---|
1223 | tRatB = ((tauResB + tauMax) / (tauResB + tauMin)) * (tauMin / tauMax); |
---|
1224 | aLowB = atan( (tauMin - tauResB) / widResB); |
---|
1225 | aUppB = atan( (tauMax - tauResB) / widResB); |
---|
1226 | } |
---|
1227 | |
---|
1228 | // Contributions from 1 / (1 - tau) for lepton beams. |
---|
1229 | double aLowT = 0.; |
---|
1230 | double aUppT = 0.; |
---|
1231 | if (hasLeptonBeams) { |
---|
1232 | aLowT = log( max( LEPTONTAUMIN, 1. - tauMin) ); |
---|
1233 | aUppT = log( max( LEPTONTAUMIN, 1. - tauMax) ); |
---|
1234 | intTau6 = aLowT - aUppT; |
---|
1235 | } |
---|
1236 | |
---|
1237 | // Select according to 1/tau or 1/tau^2. |
---|
1238 | if (iTau == 0) tau = tauMin * pow( tauMax / tauMin, tauVal); |
---|
1239 | else if (iTau == 1) tau = tauMax * tauMin |
---|
1240 | / (tauMin + (tauMax - tauMin) * tauVal); |
---|
1241 | |
---|
1242 | // Select according to 1 / (1 - tau) for lepton beams. |
---|
1243 | else if (hasLeptonBeams && iTau == nTau - 1) |
---|
1244 | tau = 1. - exp( aUppT + intTau6 * tauVal ); |
---|
1245 | |
---|
1246 | // Select according to 1 / (tau * (tau + tauRes)) or |
---|
1247 | // 1 / ((tau - tauRes)^2 + widRes^2) for resonances A and B. |
---|
1248 | else if (iTau == 2) tau = tauResA * tauMin |
---|
1249 | / ((tauResA + tauMin) * pow( tRatA, tauVal) - tauMin); |
---|
1250 | else if (iTau == 3) tau = tauResA + widResA |
---|
1251 | * tan( aLowA + (aUppA - aLowA) * tauVal); |
---|
1252 | else if (iTau == 4) tau = tauResB * tauMin |
---|
1253 | / ((tauResB + tauMin) * pow( tRatB, tauVal) - tauMin); |
---|
1254 | else if (iTau == 5) tau = tauResB + widResB |
---|
1255 | * tan( aLowB + (aUppB - aLowB) * tauVal); |
---|
1256 | |
---|
1257 | // Phase space weight in tau. |
---|
1258 | intTau0 = log( tauMax / tauMin); |
---|
1259 | intTau1 = (tauMax - tauMin) / (tauMax * tauMin); |
---|
1260 | double invWtTau = (tauCoef[0] / intTau0) + (tauCoef[1] / intTau1) / tau; |
---|
1261 | if (idResA != 0) { |
---|
1262 | intTau2 = -log(tRatA) / tauResA; |
---|
1263 | intTau3 = (aUppA - aLowA) / widResA; |
---|
1264 | invWtTau += (tauCoef[2] / intTau2) / (tau + tauResA) |
---|
1265 | + (tauCoef[3] / intTau3) * tau / ( pow2(tau - tauResA) + pow2(widResA) ); |
---|
1266 | } |
---|
1267 | if (idResB != 0) { |
---|
1268 | intTau4 = -log(tRatB) / tauResB; |
---|
1269 | intTau5 = (aUppB - aLowB) / widResB; |
---|
1270 | invWtTau += (tauCoef[4] / intTau4) / (tau + tauResB) |
---|
1271 | + (tauCoef[5] / intTau5) * tau / ( pow2(tau - tauResB) + pow2(widResB) ); |
---|
1272 | } |
---|
1273 | if (hasLeptonBeams) |
---|
1274 | invWtTau += (tauCoef[nTau - 1] / intTau6) |
---|
1275 | * tau / max( LEPTONTAUMIN, 1. - tau); |
---|
1276 | wtTau = 1. / invWtTau; |
---|
1277 | |
---|
1278 | // Calculate sHat and absolute momentum of outgoing partons. |
---|
1279 | sH = tau * s; |
---|
1280 | mHat = sqrt(sH); |
---|
1281 | if (is2) { |
---|
1282 | p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH; |
---|
1283 | pAbs = sqrtpos( p2Abs ); |
---|
1284 | } |
---|
1285 | |
---|
1286 | } |
---|
1287 | |
---|
1288 | //-------------------------------------------------------------------------- |
---|
1289 | |
---|
1290 | // Select y according to a choice of shapes. |
---|
1291 | |
---|
1292 | void PhaseSpace::selectY(int iY, double yVal) { |
---|
1293 | |
---|
1294 | // Trivial reply for unresolved lepton beams. |
---|
1295 | if (hasPointLeptons) { |
---|
1296 | y = 0.; |
---|
1297 | wtY = 1.; |
---|
1298 | x1H = 1.; |
---|
1299 | x2H = 1.; |
---|
1300 | return; |
---|
1301 | } |
---|
1302 | |
---|
1303 | // For lepton beams skip options 3&4 and go straight to 5&6. |
---|
1304 | if (hasLeptonBeams && iY > 2) iY += 2; |
---|
1305 | |
---|
1306 | // Standard expressions used below. |
---|
1307 | double expYMax = exp( yMax ); |
---|
1308 | double expYMin = exp(-yMax ); |
---|
1309 | double atanMax = atan( expYMax ); |
---|
1310 | double atanMin = atan( expYMin ); |
---|
1311 | double aUppY = (hasLeptonBeams) |
---|
1312 | ? log( max( LEPTONXMIN, LEPTONXMAX / tau - 1. ) ) : 0.; |
---|
1313 | double aLowY = LEPTONXLOGMIN; |
---|
1314 | |
---|
1315 | // 1 / cosh(y). |
---|
1316 | if (iY == 0) y = log( tan( atanMin + (atanMax - atanMin) * yVal ) ); |
---|
1317 | |
---|
1318 | // y - y_min or mirrored y_max - y. |
---|
1319 | else if (iY <= 2) y = yMax * (2. * sqrt(yVal) - 1.); |
---|
1320 | |
---|
1321 | // exp(y) or mirrored exp(-y). |
---|
1322 | else if (iY <= 4) y = log( expYMin + (expYMax - expYMin) * yVal ); |
---|
1323 | |
---|
1324 | // 1 / (1 - exp(y - y_max)) or mirrored 1 / (1 - exp(y_min - y)). |
---|
1325 | else y = yMax - log( 1. + exp(aLowY + (aUppY - aLowY) * yVal) ); |
---|
1326 | |
---|
1327 | // Mirror two cases. |
---|
1328 | if (iY == 2 || iY == 4 || iY == 6) y = -y; |
---|
1329 | |
---|
1330 | // Phase space integral in y. |
---|
1331 | intY0 = 2. * (atanMax - atanMin); |
---|
1332 | intY12 = 0.5 * pow2(2. * yMax); |
---|
1333 | intY34 = expYMax - expYMin; |
---|
1334 | intY56 = aUppY - aLowY; |
---|
1335 | double invWtY = (yCoef[0] / intY0) / cosh(y) |
---|
1336 | + (yCoef[1] / intY12) * (y + yMax) + (yCoef[2] / intY12) * (yMax - y); |
---|
1337 | if (!hasLeptonBeams) invWtY |
---|
1338 | += (yCoef[3] / intY34) * exp(y) + (yCoef[4] / intY34) * exp(-y); |
---|
1339 | else invWtY |
---|
1340 | += (yCoef[3] / intY56) / max( LEPTONXMIN, 1. - exp( y - yMax) ) |
---|
1341 | + (yCoef[4] / intY56) / max( LEPTONXMIN, 1. - exp(-y - yMax) ); |
---|
1342 | wtY = 1. / invWtY; |
---|
1343 | |
---|
1344 | // Calculate x1 and x2. |
---|
1345 | x1H = sqrt(tau) * exp(y); |
---|
1346 | x2H = sqrt(tau) * exp(-y); |
---|
1347 | } |
---|
1348 | |
---|
1349 | //-------------------------------------------------------------------------- |
---|
1350 | |
---|
1351 | // Select z = cos(theta) according to a choice of shapes. |
---|
1352 | // The selection is split in the positive- and negative-z regions, |
---|
1353 | // since a pTmax cut can remove the region around z = 0. |
---|
1354 | |
---|
1355 | void PhaseSpace::selectZ(int iZ, double zVal) { |
---|
1356 | |
---|
1357 | // Mass-dependent dampening of pT -> 0 limit. |
---|
1358 | ratio34 = max(TINY, 2. * s3 * s4 / pow2(sH)); |
---|
1359 | unity34 = 1. + ratio34; |
---|
1360 | double ratiopT2 = 2. * pT2HatMin / max( SHATMINZ, sH); |
---|
1361 | if (ratiopT2 < PT2RATMINZ) ratio34 = max( ratio34, ratiopT2); |
---|
1362 | |
---|
1363 | // Common expressions in z limits. |
---|
1364 | double zPosMax = max(ratio34, unity34 + zMax); |
---|
1365 | double zNegMax = max(ratio34, unity34 - zMax); |
---|
1366 | double zPosMin = max(ratio34, unity34 + zMin); |
---|
1367 | double zNegMin = max(ratio34, unity34 - zMin); |
---|
1368 | |
---|
1369 | // Flat in z. |
---|
1370 | if (iZ == 0) { |
---|
1371 | if (zVal < 0.5) z = -(zMax + (zMin - zMax) * 2. * zVal); |
---|
1372 | else z = zMin + (zMax - zMin) * (2. * zVal - 1.); |
---|
1373 | |
---|
1374 | // 1 / (unity34 - z). |
---|
1375 | } else if (iZ == 1) { |
---|
1376 | double areaNeg = log(zPosMax / zPosMin); |
---|
1377 | double areaPos = log(zNegMin / zNegMax); |
---|
1378 | double area = areaNeg + areaPos; |
---|
1379 | if (zVal * area < areaNeg) { |
---|
1380 | double zValMod = zVal * area / areaNeg; |
---|
1381 | z = unity34 - zPosMax * pow(zPosMin / zPosMax, zValMod); |
---|
1382 | } else { |
---|
1383 | double zValMod = (zVal * area - areaNeg)/ areaPos; |
---|
1384 | z = unity34 - zNegMin * pow(zNegMax / zNegMin, zValMod); |
---|
1385 | } |
---|
1386 | |
---|
1387 | // 1 / (unity34 + z). |
---|
1388 | } else if (iZ == 2) { |
---|
1389 | double areaNeg = log(zNegMin / zNegMax); |
---|
1390 | double areaPos = log(zPosMax / zPosMin); |
---|
1391 | double area = areaNeg + areaPos; |
---|
1392 | if (zVal * area < areaNeg) { |
---|
1393 | double zValMod = zVal * area / areaNeg; |
---|
1394 | z = zNegMax * pow(zNegMin / zNegMax, zValMod) - unity34; |
---|
1395 | } else { |
---|
1396 | double zValMod = (zVal * area - areaNeg)/ areaPos; |
---|
1397 | z = zPosMin * pow(zPosMax / zPosMin, zValMod) - unity34; |
---|
1398 | } |
---|
1399 | |
---|
1400 | // 1 / (unity34 - z)^2. |
---|
1401 | } else if (iZ == 3) { |
---|
1402 | double areaNeg = 1. / zPosMin - 1. / zPosMax; |
---|
1403 | double areaPos = 1. / zNegMax - 1. / zNegMin; |
---|
1404 | double area = areaNeg + areaPos; |
---|
1405 | if (zVal * area < areaNeg) { |
---|
1406 | double zValMod = zVal * area / areaNeg; |
---|
1407 | z = unity34 - 1. / (1./zPosMax + areaNeg * zValMod); |
---|
1408 | } else { |
---|
1409 | double zValMod = (zVal * area - areaNeg)/ areaPos; |
---|
1410 | z = unity34 - 1. / (1./zNegMin + areaPos * zValMod); |
---|
1411 | } |
---|
1412 | |
---|
1413 | // 1 / (unity34 + z)^2. |
---|
1414 | } else if (iZ == 4) { |
---|
1415 | double areaNeg = 1. / zNegMax - 1. / zNegMin; |
---|
1416 | double areaPos = 1. / zPosMin - 1. / zPosMax; |
---|
1417 | double area = areaNeg + areaPos; |
---|
1418 | if (zVal * area < areaNeg) { |
---|
1419 | double zValMod = zVal * area / areaNeg; |
---|
1420 | z = 1. / (1./zNegMax - areaNeg * zValMod) - unity34; |
---|
1421 | } else { |
---|
1422 | double zValMod = (zVal * area - areaNeg)/ areaPos; |
---|
1423 | z = 1. / (1./zPosMin - areaPos * zValMod) - unity34; |
---|
1424 | } |
---|
1425 | } |
---|
1426 | |
---|
1427 | // Safety check for roundoff errors. Combinations with z. |
---|
1428 | if (z < 0.) z = min(-zMin, max(-zMax, z)); |
---|
1429 | else z = min(zMax, max(zMin, z)); |
---|
1430 | zNeg = max(ratio34, unity34 - z); |
---|
1431 | zPos = max(ratio34, unity34 + z); |
---|
1432 | |
---|
1433 | // Phase space integral in z. |
---|
1434 | double intZ0 = 2. * (zMax - zMin); |
---|
1435 | double intZ12 = log( (zPosMax * zNegMin) / (zPosMin * zNegMax) ); |
---|
1436 | double intZ34 = 1. / zPosMin - 1. / zPosMax + 1. / zNegMax |
---|
1437 | - 1. / zNegMin; |
---|
1438 | wtZ = mHat * pAbs / ( (zCoef[0] / intZ0) + (zCoef[1] / intZ12) / zNeg |
---|
1439 | + (zCoef[2] / intZ12) / zPos + (zCoef[3] / intZ34) / pow2(zNeg) |
---|
1440 | + (zCoef[4] / intZ34) / pow2(zPos) ); |
---|
1441 | |
---|
1442 | // Calculate tHat and uHat. Also gives pTHat. |
---|
1443 | double sH34 = -0.5 * (sH - s3 - s4); |
---|
1444 | tH = sH34 + mHat * pAbs * z; |
---|
1445 | uH = sH34 - mHat * pAbs * z; |
---|
1446 | pTH = sqrtpos( (tH * uH - s3 * s4) / sH); |
---|
1447 | |
---|
1448 | } |
---|
1449 | |
---|
1450 | //-------------------------------------------------------------------------- |
---|
1451 | |
---|
1452 | // Select three-body phase space according to a cylindrically based form |
---|
1453 | // that can be chosen to favour low pT based on the form of propagators. |
---|
1454 | |
---|
1455 | bool PhaseSpace::select3Body() { |
---|
1456 | |
---|
1457 | // Upper and lower limits of pT choice for 4 and 5. |
---|
1458 | double m35S = pow2(m3 + m5); |
---|
1459 | double pT4Smax = 0.25 * ( pow2(sH - s4 - m35S) - 4. * s4 * m35S ) / sH; |
---|
1460 | if (pTHatMax > pTHatMin) pT4Smax = min( pT2HatMax, pT4Smax); |
---|
1461 | double pT4Smin = pT2HatMin; |
---|
1462 | double m34S = pow2(m3 + m4); |
---|
1463 | double pT5Smax = 0.25 * ( pow2(sH - s5 - m34S) - 4. * s5 * m34S ) / sH; |
---|
1464 | if (pTHatMax > pTHatMin) pT5Smax = min( pT2HatMax, pT5Smax); |
---|
1465 | double pT5Smin = pT2HatMin; |
---|
1466 | |
---|
1467 | // Check that pT ranges not closed. |
---|
1468 | if ( pT4Smax < pow2(pTHatMin + MASSMARGIN) ) return false; |
---|
1469 | if ( pT5Smax < pow2(pTHatMin + MASSMARGIN) ) return false; |
---|
1470 | |
---|
1471 | // Select pT4S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2. |
---|
1472 | double pTSmaxProp = pT4Smax + sTchan1; |
---|
1473 | double pTSminProp = pT4Smin + sTchan1; |
---|
1474 | double pTSratProp = pTSmaxProp / pTSminProp; |
---|
1475 | double pTSdiff = pT4Smax - pT4Smin; |
---|
1476 | double rShape = rndmPtr->flat(); |
---|
1477 | double pT4S = 0.; |
---|
1478 | if (rShape < frac3Flat) pT4S = pT4Smin + rndmPtr->flat() * pTSdiff; |
---|
1479 | else if (rShape < frac3Flat + frac3Pow1) pT4S = max( pT2HatMin, |
---|
1480 | pTSminProp * pow( pTSratProp, rndmPtr->flat() ) - sTchan1 ); |
---|
1481 | else pT4S = max( pT2HatMin, pTSminProp * pTSmaxProp |
---|
1482 | / (pTSminProp + rndmPtr->flat()* pTSdiff) - sTchan1 ); |
---|
1483 | double wt4 = pTSdiff / ( frac3Flat |
---|
1484 | + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT4S + sTchan1)) |
---|
1485 | + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT4S + sTchan1) ); |
---|
1486 | |
---|
1487 | // Select pT5S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2. |
---|
1488 | pTSmaxProp = pT5Smax + sTchan2; |
---|
1489 | pTSminProp = pT5Smin + sTchan2; |
---|
1490 | pTSratProp = pTSmaxProp / pTSminProp; |
---|
1491 | pTSdiff = pT5Smax - pT5Smin; |
---|
1492 | rShape = rndmPtr->flat(); |
---|
1493 | double pT5S = 0.; |
---|
1494 | if (rShape < frac3Flat) pT5S = pT5Smin + rndmPtr->flat() * pTSdiff; |
---|
1495 | else if (rShape < frac3Flat + frac3Pow1) pT5S = max( pT2HatMin, |
---|
1496 | pTSminProp * pow( pTSratProp, rndmPtr->flat() ) - sTchan2 ); |
---|
1497 | else pT5S = max( pT2HatMin, pTSminProp * pTSmaxProp |
---|
1498 | / (pTSminProp + rndmPtr->flat()* pTSdiff) - sTchan2 ); |
---|
1499 | double wt5 = pTSdiff / ( frac3Flat |
---|
1500 | + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT5S + sTchan2)) |
---|
1501 | + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT5S + sTchan2) ); |
---|
1502 | |
---|
1503 | // Select azimuthal angles and check that third pT in range. |
---|
1504 | double phi4 = 2. * M_PI * rndmPtr->flat(); |
---|
1505 | double phi5 = 2. * M_PI * rndmPtr->flat(); |
---|
1506 | double pT3S = max( 0., pT4S + pT5S + 2. * sqrt(pT4S * pT5S) |
---|
1507 | * cos(phi4 - phi5) ); |
---|
1508 | if ( pT3S < pT2HatMin || (pTHatMax > pTHatMin && pT3S > pT2HatMax) ) |
---|
1509 | return false; |
---|
1510 | |
---|
1511 | // Calculate transverse masses and check that phase space not closed. |
---|
1512 | double sT3 = s3 + pT3S; |
---|
1513 | double sT4 = s4 + pT4S; |
---|
1514 | double sT5 = s5 + pT5S; |
---|
1515 | double mT3 = sqrt(sT3); |
---|
1516 | double mT4 = sqrt(sT4); |
---|
1517 | double mT5 = sqrt(sT5); |
---|
1518 | if ( mT3 + mT4 + mT5 + MASSMARGIN > mHat ) return false; |
---|
1519 | |
---|
1520 | // Select rapidity for particle 3 and check that phase space not closed. |
---|
1521 | double m45S = pow2(mT4 + mT5); |
---|
1522 | double y3max = log( ( sH + sT3 - m45S + sqrtpos( pow2(sH - sT3 - m45S) |
---|
1523 | - 4 * sT3 * m45S ) ) / (2. * mHat * mT3) ); |
---|
1524 | if (y3max < YRANGEMARGIN) return false; |
---|
1525 | double y3 = (2. * rndmPtr->flat() - 1.) * (1. - YRANGEMARGIN) * y3max; |
---|
1526 | double pz3 = mT3 * sinh(y3); |
---|
1527 | double e3 = mT3 * cosh(y3); |
---|
1528 | |
---|
1529 | // Find momentum transfers in the two mirror solutions (in 4-5 frame). |
---|
1530 | double pz45 = -pz3; |
---|
1531 | double e45 = mHat - e3; |
---|
1532 | double sT45 = e45 * e45 - pz45 * pz45; |
---|
1533 | double lam45 = sqrtpos( pow2(sT45 - sT4 - sT5) - 4. * sT4 * sT5 ); |
---|
1534 | if (lam45 < YRANGEMARGIN * sH) return false; |
---|
1535 | double lam4e = sT45 + sT4 - sT5; |
---|
1536 | double lam5e = sT45 + sT5 - sT4; |
---|
1537 | double tFac = -0.5 * mHat / sT45; |
---|
1538 | double t1Pos = tFac * (e45 - pz45) * (lam4e - lam45); |
---|
1539 | double t1Neg = tFac * (e45 - pz45) * (lam4e + lam45); |
---|
1540 | double t2Pos = tFac * (e45 + pz45) * (lam5e - lam45); |
---|
1541 | double t2Neg = tFac * (e45 + pz45) * (lam5e + lam45); |
---|
1542 | |
---|
1543 | // Construct relative mirror weights and make choice. |
---|
1544 | double wtPosUnnorm = 1.; |
---|
1545 | double wtNegUnnorm = 1.; |
---|
1546 | if (useMirrorWeight) { |
---|
1547 | wtPosUnnorm = 1./ pow2( (t1Pos - sTchan1) * (t2Pos - sTchan2) ); |
---|
1548 | wtNegUnnorm = 1./ pow2( (t1Neg - sTchan1) * (t2Neg - sTchan2) ); |
---|
1549 | } |
---|
1550 | double wtPos = wtPosUnnorm / (wtPosUnnorm + wtNegUnnorm); |
---|
1551 | double wtNeg = wtNegUnnorm / (wtPosUnnorm + wtNegUnnorm); |
---|
1552 | double epsilon = (rndmPtr->flat() < wtPos) ? 1. : -1.; |
---|
1553 | |
---|
1554 | // Construct four-vectors in rest frame of subprocess. |
---|
1555 | double px4 = sqrt(pT4S) * cos(phi4); |
---|
1556 | double py4 = sqrt(pT4S) * sin(phi4); |
---|
1557 | double px5 = sqrt(pT5S) * cos(phi5); |
---|
1558 | double py5 = sqrt(pT5S) * sin(phi5); |
---|
1559 | double pz4 = 0.5 * (pz45 * lam4e + epsilon * e45 * lam45) / sT45; |
---|
1560 | double pz5 = pz45 - pz4; |
---|
1561 | double e4 = sqrt(sT4 + pz4 * pz4); |
---|
1562 | double e5 = sqrt(sT5 + pz5 * pz5); |
---|
1563 | p3cm = Vec4( -(px4 + px5), -(py4 + py5), pz3, e3); |
---|
1564 | p4cm = Vec4( px4, py4, pz4, e4); |
---|
1565 | p5cm = Vec4( px5, py5, pz5, e5); |
---|
1566 | |
---|
1567 | // Total weight to associate with kinematics choice. |
---|
1568 | wt3Body = wt4 * wt5 * (2. * y3max) / (128. * pow3(M_PI) * lam45); |
---|
1569 | wt3Body *= (epsilon > 0.) ? 1. / wtPos : 1. / wtNeg; |
---|
1570 | |
---|
1571 | // Cross section of form |M|^2/(2 sHat) dPS_3 so need 1/(2 sHat). |
---|
1572 | wt3Body /= (2. * sH); |
---|
1573 | |
---|
1574 | // Done. |
---|
1575 | return true; |
---|
1576 | |
---|
1577 | } |
---|
1578 | |
---|
1579 | //-------------------------------------------------------------------------- |
---|
1580 | |
---|
1581 | // Solve linear equation system for better phase space coefficients. |
---|
1582 | |
---|
1583 | void PhaseSpace::solveSys( int n, int bin[8], double vec[8], |
---|
1584 | double mat[8][8], double coef[8], ostream& os) { |
---|
1585 | |
---|
1586 | // Optional printout. |
---|
1587 | if (showSearch) { |
---|
1588 | os << "\n Equation system: " << setw(5) << bin[0]; |
---|
1589 | for (int j = 0; j < n; ++j) os << setw(12) << mat[0][j]; |
---|
1590 | os << setw(12) << vec[0] << "\n"; |
---|
1591 | for (int i = 1; i < n; ++i) { |
---|
1592 | os << " " << setw(5) << bin[i]; |
---|
1593 | for (int j = 0; j < n; ++j) os << setw(12) << mat[i][j]; |
---|
1594 | os << setw(12) << vec[i] << "\n"; |
---|
1595 | } |
---|
1596 | } |
---|
1597 | |
---|
1598 | // Local variables. |
---|
1599 | double vecNor[8], coefTmp[8]; |
---|
1600 | for (int i = 0; i < n; ++i) coefTmp[i] = 0.; |
---|
1601 | |
---|
1602 | // Check if equation system solvable. |
---|
1603 | bool canSolve = true; |
---|
1604 | for (int i = 0; i < n; ++i) if (bin[i] == 0) canSolve = false; |
---|
1605 | double vecSum = 0.; |
---|
1606 | for (int i = 0; i < n; ++i) vecSum += vec[i]; |
---|
1607 | if (abs(vecSum) < TINY) canSolve = false; |
---|
1608 | |
---|
1609 | // Solve to find relative importance of cross-section pieces. |
---|
1610 | if (canSolve) { |
---|
1611 | for (int i = 0; i < n; ++i) vecNor[i] = max( 0.1, vec[i] / vecSum); |
---|
1612 | for (int k = 0; k < n - 1; ++k) { |
---|
1613 | for (int i = k + 1; i < n; ++i) { |
---|
1614 | if (abs(mat[k][k]) < TINY) {canSolve = false; break;} |
---|
1615 | double ratio = mat[i][k] / mat[k][k]; |
---|
1616 | vec[i] -= ratio * vec[k]; |
---|
1617 | for (int j = k; j < n; ++j) mat[i][j] -= ratio * mat[k][j]; |
---|
1618 | } |
---|
1619 | if (!canSolve) break; |
---|
1620 | } |
---|
1621 | if (canSolve) { |
---|
1622 | for (int k = n - 1; k >= 0; --k) { |
---|
1623 | for (int j = k + 1; j < n; ++j) vec[k] -= mat[k][j] * coefTmp[j]; |
---|
1624 | coefTmp[k] = vec[k] / mat[k][k]; |
---|
1625 | } |
---|
1626 | } |
---|
1627 | } |
---|
1628 | |
---|
1629 | // Share evenly if failure. |
---|
1630 | if (!canSolve) for (int i = 0; i < n; ++i) { |
---|
1631 | coefTmp[i] = 1.; |
---|
1632 | vecNor[i] = 0.1; |
---|
1633 | if (vecSum > TINY) vecNor[i] = max(0.1, vec[i] / vecSum); |
---|
1634 | } |
---|
1635 | |
---|
1636 | // Normalize coefficients, with piece shared democratically. |
---|
1637 | double coefSum = 0.; |
---|
1638 | vecSum = 0.; |
---|
1639 | for (int i = 0; i < n; ++i) { |
---|
1640 | coefTmp[i] = max( 0., coefTmp[i]); |
---|
1641 | coefSum += coefTmp[i]; |
---|
1642 | vecSum += vecNor[i]; |
---|
1643 | } |
---|
1644 | if (coefSum > 0.) for (int i = 0; i < n; ++i) coef[i] = EVENFRAC / n |
---|
1645 | + (1. - EVENFRAC) * 0.5 * (coefTmp[i] / coefSum + vecNor[i] / vecSum); |
---|
1646 | else for (int i = 0; i < n; ++i) coef[i] = 1. / n; |
---|
1647 | |
---|
1648 | // Optional printout. |
---|
1649 | if (showSearch) { |
---|
1650 | os << " Solution: "; |
---|
1651 | for (int i = 0; i < n; ++i) os << setw(12) << coef[i]; |
---|
1652 | os << "\n"; |
---|
1653 | } |
---|
1654 | } |
---|
1655 | |
---|
1656 | //-------------------------------------------------------------------------- |
---|
1657 | |
---|
1658 | // Setup mass selection for one resonance at a time - part 1. |
---|
1659 | |
---|
1660 | void PhaseSpace::setupMass1(int iM) { |
---|
1661 | |
---|
1662 | // Identity for mass seletion; is 0 also for light quarks (not yet selected). |
---|
1663 | if (iM == 3) idMass[iM] = abs(sigmaProcessPtr->id3Mass()); |
---|
1664 | if (iM == 4) idMass[iM] = abs(sigmaProcessPtr->id4Mass()); |
---|
1665 | if (iM == 5) idMass[iM] = abs(sigmaProcessPtr->id5Mass()); |
---|
1666 | |
---|
1667 | // Masses and widths of resonances. |
---|
1668 | if (idMass[iM] == 0) { |
---|
1669 | mPeak[iM] = 0.; |
---|
1670 | mWidth[iM] = 0.; |
---|
1671 | mMin[iM] = 0.; |
---|
1672 | mMax[iM] = 0.; |
---|
1673 | } else { |
---|
1674 | mPeak[iM] = particleDataPtr->m0(idMass[iM]); |
---|
1675 | mWidth[iM] = particleDataPtr->mWidth(idMass[iM]); |
---|
1676 | mMin[iM] = particleDataPtr->mMin(idMass[iM]); |
---|
1677 | mMax[iM] = particleDataPtr->mMax(idMass[iM]); |
---|
1678 | // gmZmode == 1 means pure photon propagator; set at lower mass limit. |
---|
1679 | if (idMass[iM] == 23 && gmZmode == 1) mPeak[iM] = mMin[iM]; |
---|
1680 | } |
---|
1681 | |
---|
1682 | // Mass and width combinations for Breit-Wigners. |
---|
1683 | sPeak[iM] = mPeak[iM] * mPeak[iM]; |
---|
1684 | useBW[iM] = useBreitWigners && (mWidth[iM] > minWidthBreitWigners); |
---|
1685 | if (!useBW[iM]) mWidth[iM] = 0.; |
---|
1686 | mw[iM] = mPeak[iM] * mWidth[iM]; |
---|
1687 | wmRat[iM] = (idMass[iM] == 0 || mPeak[iM] == 0.) |
---|
1688 | ? 0. : mWidth[iM] / mPeak[iM]; |
---|
1689 | |
---|
1690 | // Simple Breit-Wigner range, upper edge to be corrected subsequently. |
---|
1691 | if (useBW[iM]) { |
---|
1692 | mLower[iM] = mMin[iM]; |
---|
1693 | mUpper[iM] = mHatMax; |
---|
1694 | } |
---|
1695 | |
---|
1696 | } |
---|
1697 | |
---|
1698 | //-------------------------------------------------------------------------- |
---|
1699 | |
---|
1700 | // Setup mass selection for one resonance at a time - part 2. |
---|
1701 | |
---|
1702 | void PhaseSpace::setupMass2(int iM, double distToThresh) { |
---|
1703 | |
---|
1704 | // Store reduced Breit-Wigner range. |
---|
1705 | if (mMax[iM] > mMin[iM]) mUpper[iM] = min( mUpper[iM], mMax[iM]); |
---|
1706 | sLower[iM] = mLower[iM] * mLower[iM]; |
---|
1707 | sUpper[iM] = mUpper[iM] * mUpper[iM]; |
---|
1708 | |
---|
1709 | // Prepare to select m3 by BW + flat + 1/s_3. |
---|
1710 | // Determine relative coefficients by allowed mass range. |
---|
1711 | if (distToThresh > THRESHOLDSIZE) { |
---|
1712 | fracFlat[iM] = 0.1; |
---|
1713 | fracInv[iM] = 0.1; |
---|
1714 | } else if (distToThresh > - THRESHOLDSIZE) { |
---|
1715 | fracFlat[iM] = 0.25 - 0.15 * distToThresh / THRESHOLDSIZE; |
---|
1716 | fracInv [iM] = 0.15 - 0.05 * distToThresh / THRESHOLDSIZE; |
---|
1717 | } else { |
---|
1718 | fracFlat[iM] = 0.4; |
---|
1719 | fracInv[iM] = 0.2; |
---|
1720 | } |
---|
1721 | |
---|
1722 | // For gamma*/Z0: increase 1/s_i part and introduce 1/s_i^2 part. |
---|
1723 | fracInv2[iM] = 0.; |
---|
1724 | if (idMass[iM] == 23 && gmZmode == 0) { |
---|
1725 | fracFlat[iM] *= 0.5; |
---|
1726 | fracInv[iM] = 0.5 * fracInv[iM] + 0.25; |
---|
1727 | fracInv2[iM] = 0.25; |
---|
1728 | } else if (idMass[iM] == 23 && gmZmode == 1) { |
---|
1729 | fracFlat[iM] = 0.1; |
---|
1730 | fracInv[iM] = 0.4; |
---|
1731 | fracInv2[iM] = 0.4; |
---|
1732 | } |
---|
1733 | |
---|
1734 | // Normalization integrals for the respective contribution. |
---|
1735 | atanLower[iM] = atan( (sLower[iM] - sPeak[iM])/ mw[iM] ); |
---|
1736 | atanUpper[iM] = atan( (sUpper[iM] - sPeak[iM])/ mw[iM] ); |
---|
1737 | intBW[iM] = atanUpper[iM] - atanLower[iM]; |
---|
1738 | intFlat[iM] = sUpper[iM] - sLower[iM]; |
---|
1739 | intInv[iM] = log( sUpper[iM] / sLower[iM] ); |
---|
1740 | intInv2[iM] = 1./sLower[iM] - 1./sUpper[iM]; |
---|
1741 | |
---|
1742 | } |
---|
1743 | |
---|
1744 | //-------------------------------------------------------------------------- |
---|
1745 | |
---|
1746 | // Select Breit-Wigner-distributed or fixed masses. |
---|
1747 | |
---|
1748 | void PhaseSpace::trialMass(int iM) { |
---|
1749 | |
---|
1750 | // References to masses to be set. |
---|
1751 | double& mSet = (iM == 3) ? m3 : ( (iM == 4) ? m4 : m5 ); |
---|
1752 | double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 ); |
---|
1753 | |
---|
1754 | // Distribution for m_i is BW + flat + 1/s_i + 1/s_i^2. |
---|
1755 | if (useBW[iM]) { |
---|
1756 | double pickForm = rndmPtr->flat(); |
---|
1757 | if (pickForm > fracFlat[iM] + fracInv[iM] + fracInv2[iM]) |
---|
1758 | sSet = sPeak[iM] + mw[iM] * tan( atanLower[iM] |
---|
1759 | + rndmPtr->flat() * intBW[iM] ); |
---|
1760 | else if (pickForm > fracInv[iM] + fracInv2[iM]) |
---|
1761 | sSet = sLower[iM] + rndmPtr->flat() * (sUpper[iM] - sLower[iM]); |
---|
1762 | else if (pickForm > fracInv2[iM]) |
---|
1763 | sSet = sLower[iM] * pow( sUpper[iM] / sLower[iM], rndmPtr->flat() ); |
---|
1764 | else sSet = sLower[iM] * sUpper[iM] |
---|
1765 | / (sLower[iM] + rndmPtr->flat() * (sUpper[iM] - sLower[iM])); |
---|
1766 | mSet = sqrt(sSet); |
---|
1767 | |
---|
1768 | // Else m_i is fixed at peak value. |
---|
1769 | } else { |
---|
1770 | mSet = mPeak[iM]; |
---|
1771 | sSet = sPeak[iM]; |
---|
1772 | } |
---|
1773 | |
---|
1774 | } |
---|
1775 | |
---|
1776 | //-------------------------------------------------------------------------- |
---|
1777 | |
---|
1778 | // Naively a fixed-width Breit-Wigner is used to pick the mass. |
---|
1779 | // Here come the correction factors for |
---|
1780 | // (i) preselection according to BW + flat in s_i + 1/s_i + 1/s_i^2, |
---|
1781 | // (ii) reduced allowed mass range, |
---|
1782 | // (iii) running width, i.e. m0*Gamma0 -> s*Gamma0/m0. |
---|
1783 | // In the end, the weighted distribution is a running-width BW. |
---|
1784 | |
---|
1785 | double PhaseSpace::weightMass(int iM) { |
---|
1786 | |
---|
1787 | // Reference to mass and to Breit-Wigner weight to be set. |
---|
1788 | double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 ); |
---|
1789 | double& runBWH = (iM == 3) ? runBW3H : ( (iM == 4) ? runBW4H : runBW5H ); |
---|
1790 | |
---|
1791 | // Default weight if no Breit-Wigner. |
---|
1792 | runBWH = 1.; |
---|
1793 | if (!useBW[iM]) return 1.; |
---|
1794 | |
---|
1795 | // Weight of generated distribution. |
---|
1796 | double genBW = (1. - fracFlat[iM] - fracInv[iM] - fracInv2[iM]) |
---|
1797 | * mw[iM] / ( (pow2(sSet - sPeak[iM]) + pow2(mw[iM])) * intBW[iM]) |
---|
1798 | + fracFlat[iM] / intFlat[iM] + fracInv[iM] / (sSet * intInv[iM]) |
---|
1799 | + fracInv2[iM] / (sSet*sSet * intInv2[iM]); |
---|
1800 | |
---|
1801 | // Weight of distribution with running width in Breit-Wigner. |
---|
1802 | double mwRun = sSet * wmRat[iM]; |
---|
1803 | runBWH = mwRun / (pow2(sSet - sPeak[iM]) + pow2(mwRun)) / M_PI; |
---|
1804 | |
---|
1805 | // Done. |
---|
1806 | return (runBWH / genBW); |
---|
1807 | |
---|
1808 | } |
---|
1809 | |
---|
1810 | //========================================================================== |
---|
1811 | |
---|
1812 | // PhaseSpace2to1tauy class. |
---|
1813 | // 2 -> 1 kinematics for normal subprocesses. |
---|
1814 | |
---|
1815 | //-------------------------------------------------------------------------- |
---|
1816 | |
---|
1817 | // Set limits for resonance mass selection. |
---|
1818 | |
---|
1819 | bool PhaseSpace2to1tauy::setupMass() { |
---|
1820 | |
---|
1821 | // Treat Z0 as such or as gamma*/Z0 |
---|
1822 | gmZmode = gmZmodeGlobal; |
---|
1823 | int gmZmodeProc = sigmaProcessPtr->gmZmode(); |
---|
1824 | if (gmZmodeProc >= 0) gmZmode = gmZmodeProc; |
---|
1825 | |
---|
1826 | // Mass limits for current resonance. |
---|
1827 | int idRes = abs(sigmaProcessPtr->resonanceA()); |
---|
1828 | int idTmp = abs(sigmaProcessPtr->resonanceB()); |
---|
1829 | if (idTmp > 0) idRes = idTmp; |
---|
1830 | double mResMin = (idRes == 0) ? 0. : particleDataPtr->mMin(idRes); |
---|
1831 | double mResMax = (idRes == 0) ? 0. : particleDataPtr->mMax(idRes); |
---|
1832 | |
---|
1833 | // Compare with global mass limits and pick tighter of them. |
---|
1834 | mHatMin = max( mResMin, mHatGlobalMin); |
---|
1835 | sHatMin = mHatMin*mHatMin; |
---|
1836 | mHatMax = eCM; |
---|
1837 | if (mResMax > mResMin) mHatMax = min( mHatMax, mResMax); |
---|
1838 | if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( mHatMax, mHatGlobalMax); |
---|
1839 | sHatMax = mHatMax*mHatMax; |
---|
1840 | |
---|
1841 | // Default Breit-Wigner weight. |
---|
1842 | wtBW = 1.; |
---|
1843 | |
---|
1844 | // Fail if mass window (almost) closed. |
---|
1845 | return (mHatMax > mHatMin + MASSMARGIN); |
---|
1846 | |
---|
1847 | } |
---|
1848 | |
---|
1849 | //-------------------------------------------------------------------------- |
---|
1850 | |
---|
1851 | // Construct the four-vector kinematics from the trial values. |
---|
1852 | |
---|
1853 | bool PhaseSpace2to1tauy::finalKin() { |
---|
1854 | |
---|
1855 | // Particle masses; incoming always on mass shell. |
---|
1856 | mH[1] = 0.; |
---|
1857 | mH[2] = 0.; |
---|
1858 | mH[3] = mHat; |
---|
1859 | |
---|
1860 | // Incoming partons along beam axes. Outgoing has sum of momenta. |
---|
1861 | pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H); |
---|
1862 | pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H); |
---|
1863 | pH[3] = pH[1] + pH[2]; |
---|
1864 | |
---|
1865 | // Done. |
---|
1866 | return true; |
---|
1867 | } |
---|
1868 | |
---|
1869 | //========================================================================== |
---|
1870 | |
---|
1871 | // PhaseSpace2to2tauyz class. |
---|
1872 | // 2 -> 2 kinematics for normal subprocesses. |
---|
1873 | |
---|
1874 | //-------------------------------------------------------------------------- |
---|
1875 | |
---|
1876 | // Set up for fixed or Breit-Wigner mass selection. |
---|
1877 | |
---|
1878 | bool PhaseSpace2to2tauyz::setupMasses() { |
---|
1879 | |
---|
1880 | // Treat Z0 as such or as gamma*/Z0 |
---|
1881 | gmZmode = gmZmodeGlobal; |
---|
1882 | int gmZmodeProc = sigmaProcessPtr->gmZmode(); |
---|
1883 | if (gmZmodeProc >= 0) gmZmode = gmZmodeProc; |
---|
1884 | |
---|
1885 | // Set sHat limits - based on global limits only. |
---|
1886 | mHatMin = mHatGlobalMin; |
---|
1887 | sHatMin = mHatMin*mHatMin; |
---|
1888 | mHatMax = eCM; |
---|
1889 | if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax); |
---|
1890 | sHatMax = mHatMax*mHatMax; |
---|
1891 | |
---|
1892 | // Masses and widths of resonances. |
---|
1893 | setupMass1(3); |
---|
1894 | setupMass1(4); |
---|
1895 | |
---|
1896 | // Reduced mass range when two massive particles. |
---|
1897 | if (useBW[3]) mUpper[3] -= (useBW[4]) ? mMin[4] : mPeak[4]; |
---|
1898 | if (useBW[4]) mUpper[4] -= (useBW[3]) ? mMin[3] : mPeak[3]; |
---|
1899 | |
---|
1900 | // If closed phase space then unallowed process. |
---|
1901 | bool physical = true; |
---|
1902 | if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false; |
---|
1903 | if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false; |
---|
1904 | if (!useBW[3] && !useBW[4] && mHatMax < mPeak[3] + mPeak[4] + MASSMARGIN) |
---|
1905 | physical = false; |
---|
1906 | if (!physical) return false; |
---|
1907 | |
---|
1908 | // If either particle is massless then need extra pTHat cut. |
---|
1909 | pTHatMin = pTHatGlobalMin; |
---|
1910 | if (mPeak[3] < pTHatMinDiverge || mPeak[4] < pTHatMinDiverge) |
---|
1911 | pTHatMin = max( pTHatMin, pTHatMinDiverge); |
---|
1912 | pT2HatMin = pTHatMin * pTHatMin; |
---|
1913 | pTHatMax = pTHatGlobalMax; |
---|
1914 | pT2HatMax = pTHatMax * pTHatMax; |
---|
1915 | |
---|
1916 | // Prepare to select m3 by BW + flat + 1/s_3. |
---|
1917 | if (useBW[3]) { |
---|
1918 | double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[3] |
---|
1919 | / (pow2(mWidth[3]) + pow2(mWidth[4])); |
---|
1920 | double distToThreshB = (mHatMax - mPeak[3] - mMin[4]) / mWidth[3]; |
---|
1921 | double distToThresh = min( distToThreshA, distToThreshB); |
---|
1922 | setupMass2(3, distToThresh); |
---|
1923 | } |
---|
1924 | |
---|
1925 | // Prepare to select m4 by BW + flat + 1/s_4. |
---|
1926 | if (useBW[4]) { |
---|
1927 | double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[4] |
---|
1928 | / (pow2(mWidth[3]) + pow2(mWidth[4])); |
---|
1929 | double distToThreshB = (mHatMax - mMin[3] - mPeak[4]) / mWidth[4]; |
---|
1930 | double distToThresh = min( distToThreshA, distToThreshB); |
---|
1931 | setupMass2(4, distToThresh); |
---|
1932 | } |
---|
1933 | |
---|
1934 | // Initialization masses. Special cases when constrained phase space. |
---|
1935 | m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3]; |
---|
1936 | m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4]; |
---|
1937 | if (m3 + m4 + THRESHOLDSIZE * (mWidth[3] + mWidth[4]) + MASSMARGIN |
---|
1938 | > mHatMax) { |
---|
1939 | if (useBW[3] && useBW[4]) physical = constrainedM3M4(); |
---|
1940 | else if (useBW[3]) physical = constrainedM3(); |
---|
1941 | else if (useBW[4]) physical = constrainedM4(); |
---|
1942 | } |
---|
1943 | s3 = m3*m3; |
---|
1944 | s4 = m4*m4; |
---|
1945 | |
---|
1946 | // Correct selected mass-spectrum to running-width Breit-Wigner. |
---|
1947 | // Extra safety margin for maximum search. |
---|
1948 | wtBW = 1.; |
---|
1949 | if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX; |
---|
1950 | if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX; |
---|
1951 | |
---|
1952 | // Done. |
---|
1953 | return physical; |
---|
1954 | |
---|
1955 | } |
---|
1956 | |
---|
1957 | |
---|
1958 | //-------------------------------------------------------------------------- |
---|
1959 | |
---|
1960 | // Select Breit-Wigner-distributed or fixed masses. |
---|
1961 | |
---|
1962 | bool PhaseSpace2to2tauyz::trialMasses() { |
---|
1963 | |
---|
1964 | // By default vanishing cross section. |
---|
1965 | sigmaNw = 0.; |
---|
1966 | wtBW = 1.; |
---|
1967 | |
---|
1968 | // Pick m3 and m4 independently. |
---|
1969 | trialMass(3); |
---|
1970 | trialMass(4); |
---|
1971 | |
---|
1972 | // If outside phase space then reject event. |
---|
1973 | if (m3 + m4 + MASSMARGIN > mHatMax) return false; |
---|
1974 | |
---|
1975 | // Correct selected mass-spectrum to running-width Breit-Wigner. |
---|
1976 | if (useBW[3]) wtBW *= weightMass(3); |
---|
1977 | if (useBW[4]) wtBW *= weightMass(4); |
---|
1978 | |
---|
1979 | // Done. |
---|
1980 | return true; |
---|
1981 | } |
---|
1982 | |
---|
1983 | //-------------------------------------------------------------------------- |
---|
1984 | |
---|
1985 | // Construct the four-vector kinematics from the trial values. |
---|
1986 | |
---|
1987 | bool PhaseSpace2to2tauyz::finalKin() { |
---|
1988 | |
---|
1989 | // Assign masses to particles assumed massless in matrix elements. |
---|
1990 | int id3 = sigmaProcessPtr->id(3); |
---|
1991 | int id4 = sigmaProcessPtr->id(4); |
---|
1992 | if (idMass[3] == 0) { m3 = particleDataPtr->m0(id3); s3 = m3*m3; } |
---|
1993 | if (idMass[4] == 0) { m4 = particleDataPtr->m0(id4); s4 = m4*m4; } |
---|
1994 | |
---|
1995 | // Sometimes swap tHat <-> uHat to reflect chosen final-state order. |
---|
1996 | if (sigmaProcessPtr->swappedTU()) { |
---|
1997 | swap(tH, uH); |
---|
1998 | z = -z; |
---|
1999 | } |
---|
2000 | |
---|
2001 | // Check that phase space still open after new mass assignment. |
---|
2002 | if (m3 + m4 + MASSMARGIN > mHat) { |
---|
2003 | infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::finalKin: " |
---|
2004 | "failed after mass assignment"); |
---|
2005 | return false; |
---|
2006 | } |
---|
2007 | p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH; |
---|
2008 | pAbs = sqrtpos( p2Abs ); |
---|
2009 | |
---|
2010 | // Particle masses; incoming always on mass shell. |
---|
2011 | mH[1] = 0.; |
---|
2012 | mH[2] = 0.; |
---|
2013 | mH[3] = m3; |
---|
2014 | mH[4] = m4; |
---|
2015 | |
---|
2016 | // Incoming partons along beam axes. |
---|
2017 | pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H); |
---|
2018 | pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H); |
---|
2019 | |
---|
2020 | // Outgoing partons initially in collision CM frame along beam axes. |
---|
2021 | pH[3] = Vec4( 0., 0., pAbs, 0.5 * (sH + s3 - s4) / mHat); |
---|
2022 | pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (sH + s4 - s3) / mHat); |
---|
2023 | |
---|
2024 | // Then rotate and boost them to overall CM frame. |
---|
2025 | theta = acos(z); |
---|
2026 | phi = 2. * M_PI * rndmPtr->flat(); |
---|
2027 | betaZ = (x1H - x2H)/(x1H + x2H); |
---|
2028 | pH[3].rot( theta, phi); |
---|
2029 | pH[4].rot( theta, phi); |
---|
2030 | pH[3].bst( 0., 0., betaZ); |
---|
2031 | pH[4].bst( 0., 0., betaZ); |
---|
2032 | pTH = pAbs * sin(theta); |
---|
2033 | |
---|
2034 | // Done. |
---|
2035 | return true; |
---|
2036 | } |
---|
2037 | |
---|
2038 | //-------------------------------------------------------------------------- |
---|
2039 | |
---|
2040 | // Special choice of m3 and m4 when mHatMax push them off mass shell. |
---|
2041 | // Vary x in expression m3 + m4 = mHatMax - x * (Gamma3 + Gamma4). |
---|
2042 | // For each x try to put either 3 or 4 as close to mass shell as possible. |
---|
2043 | // Maximize BW_3 * BW_4 * beta_34, where latter approximate phase space. |
---|
2044 | |
---|
2045 | bool PhaseSpace2to2tauyz::constrainedM3M4() { |
---|
2046 | |
---|
2047 | // Initial values. |
---|
2048 | bool foundNonZero = false; |
---|
2049 | double wtMassMax = 0.; |
---|
2050 | double m3WtMax = 0.; |
---|
2051 | double m4WtMax = 0.; |
---|
2052 | double xMax = (mHatMax - mLower[3] - mLower[4]) / (mWidth[3] + mWidth[4]); |
---|
2053 | double xStep = THRESHOLDSTEP * min(1., xMax); |
---|
2054 | double xNow = 0.; |
---|
2055 | double wtMassXbin, wtMassMaxOld, m34, mT34Min, wtMassNow, |
---|
2056 | wtBW3Now, wtBW4Now, beta34Now; |
---|
2057 | |
---|
2058 | // Step through increasing x values. |
---|
2059 | do { |
---|
2060 | xNow += xStep; |
---|
2061 | wtMassXbin = 0.; |
---|
2062 | wtMassMaxOld = wtMassMax; |
---|
2063 | m34 = mHatMax - xNow * (mWidth[3] + mWidth[4]); |
---|
2064 | |
---|
2065 | // Study point where m3 as close as possible to on-shell. |
---|
2066 | m3 = min( mUpper[3], m34 - mLower[4]); |
---|
2067 | if (m3 > mPeak[3]) m3 = max( mLower[3], mPeak[3]); |
---|
2068 | m4 = m34 - m3; |
---|
2069 | if (m4 < mLower[4]) {m4 = mLower[4]; m3 = m34 - m4;} |
---|
2070 | |
---|
2071 | // Check that inside phase space limit set by pTmin. |
---|
2072 | mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin); |
---|
2073 | if (mT34Min < mHatMax) { |
---|
2074 | |
---|
2075 | // Breit-Wigners and beta factor give total weight. |
---|
2076 | wtMassNow = 0.; |
---|
2077 | if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4] |
---|
2078 | && m4 < mUpper[4]) { |
---|
2079 | wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) ); |
---|
2080 | wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) ); |
---|
2081 | beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4) |
---|
2082 | - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax); |
---|
2083 | wtMassNow = wtBW3Now * wtBW4Now * beta34Now; |
---|
2084 | } |
---|
2085 | |
---|
2086 | // Store new maximum, if any. |
---|
2087 | if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow; |
---|
2088 | if (wtMassNow > wtMassMax) { |
---|
2089 | foundNonZero = true; |
---|
2090 | wtMassMax = wtMassNow; |
---|
2091 | m3WtMax = m3; |
---|
2092 | m4WtMax = m4; |
---|
2093 | } |
---|
2094 | } |
---|
2095 | |
---|
2096 | // Study point where m4 as close as possible to on-shell. |
---|
2097 | m4 = min( mUpper[4], m34 - mLower[3]); |
---|
2098 | if (m4 > mPeak[4]) m4 = max( mLower[4], mPeak[4]); |
---|
2099 | m3 = m34 - m4; |
---|
2100 | if (m3 < mLower[3]) {m3 = mLower[3]; m4 = m34 - m3;} |
---|
2101 | |
---|
2102 | // Check that inside phase space limit set by pTmin. |
---|
2103 | mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin); |
---|
2104 | if (mT34Min < mHatMax) { |
---|
2105 | |
---|
2106 | // Breit-Wigners and beta factor give total weight. |
---|
2107 | wtMassNow = 0.; |
---|
2108 | if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4] |
---|
2109 | && m4 < mUpper[4]) { |
---|
2110 | wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) ); |
---|
2111 | wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) ); |
---|
2112 | beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4) |
---|
2113 | - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax); |
---|
2114 | wtMassNow = wtBW3Now * wtBW4Now * beta34Now; |
---|
2115 | } |
---|
2116 | |
---|
2117 | // Store new maximum, if any. |
---|
2118 | if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow; |
---|
2119 | if (wtMassNow > wtMassMax) { |
---|
2120 | foundNonZero = true; |
---|
2121 | wtMassMax = wtMassNow; |
---|
2122 | m3WtMax = m3; |
---|
2123 | m4WtMax = m4; |
---|
2124 | } |
---|
2125 | } |
---|
2126 | |
---|
2127 | // Continue stepping if increasing trend and more x range available. |
---|
2128 | } while ( (!foundNonZero || wtMassXbin > wtMassMaxOld) |
---|
2129 | && xNow < xMax - xStep); |
---|
2130 | |
---|
2131 | // Restore best values for subsequent maximization. Return. |
---|
2132 | m3 = m3WtMax; |
---|
2133 | m4 = m4WtMax; |
---|
2134 | return foundNonZero; |
---|
2135 | |
---|
2136 | } |
---|
2137 | |
---|
2138 | //-------------------------------------------------------------------------- |
---|
2139 | |
---|
2140 | // Special choice of m3 when mHatMax pushes it off mass shell. |
---|
2141 | // Vary x in expression m3 = mHatMax - m4 - x * Gamma3. |
---|
2142 | // Maximize BW_3 * beta_34, where latter approximate phase space. |
---|
2143 | |
---|
2144 | bool PhaseSpace2to2tauyz::constrainedM3() { |
---|
2145 | |
---|
2146 | // Initial values. |
---|
2147 | bool foundNonZero = false; |
---|
2148 | double wtMassMax = 0.; |
---|
2149 | double m3WtMax = 0.; |
---|
2150 | double mT4Min = sqrt(m4*m4 + pT2HatMin); |
---|
2151 | double xMax = (mHatMax - mLower[3] - m4) / mWidth[3]; |
---|
2152 | double xStep = THRESHOLDSTEP * min(1., xMax); |
---|
2153 | double xNow = 0.; |
---|
2154 | double wtMassNow, mT34Min, wtBW3Now, beta34Now; |
---|
2155 | |
---|
2156 | // Step through increasing x values; gives m3 unambiguously. |
---|
2157 | do { |
---|
2158 | xNow += xStep; |
---|
2159 | wtMassNow = 0.; |
---|
2160 | m3 = mHatMax - m4 - xNow * mWidth[3]; |
---|
2161 | |
---|
2162 | // Check that inside phase space limit set by pTmin. |
---|
2163 | mT34Min = sqrt(m3*m3 + pT2HatMin) + mT4Min; |
---|
2164 | if (mT34Min < mHatMax) { |
---|
2165 | |
---|
2166 | // Breit-Wigner and beta factor give total weight. |
---|
2167 | wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) ); |
---|
2168 | beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4) |
---|
2169 | - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax); |
---|
2170 | wtMassNow = wtBW3Now * beta34Now; |
---|
2171 | |
---|
2172 | // Store new maximum, if any. |
---|
2173 | if (wtMassNow > wtMassMax) { |
---|
2174 | foundNonZero = true; |
---|
2175 | wtMassMax = wtMassNow; |
---|
2176 | m3WtMax = m3; |
---|
2177 | } |
---|
2178 | } |
---|
2179 | |
---|
2180 | // Continue stepping if increasing trend and more x range available. |
---|
2181 | } while ( (!foundNonZero || wtMassNow > wtMassMax) |
---|
2182 | && xNow < xMax - xStep); |
---|
2183 | |
---|
2184 | // Restore best value for subsequent maximization. Return. |
---|
2185 | m3 = m3WtMax; |
---|
2186 | return foundNonZero; |
---|
2187 | |
---|
2188 | } |
---|
2189 | |
---|
2190 | //-------------------------------------------------------------------------- |
---|
2191 | |
---|
2192 | // Special choice of m4 when mHatMax pushes it off mass shell. |
---|
2193 | // Vary x in expression m4 = mHatMax - m3 - x * Gamma4. |
---|
2194 | // Maximize BW_4 * beta_34, where latter approximate phase space. |
---|
2195 | |
---|
2196 | bool PhaseSpace2to2tauyz::constrainedM4() { |
---|
2197 | |
---|
2198 | // Initial values. |
---|
2199 | bool foundNonZero = false; |
---|
2200 | double wtMassMax = 0.; |
---|
2201 | double m4WtMax = 0.; |
---|
2202 | double mT3Min = sqrt(m3*m3 + pT2HatMin); |
---|
2203 | double xMax = (mHatMax - mLower[4] - m3) / mWidth[4]; |
---|
2204 | double xStep = THRESHOLDSTEP * min(1., xMax); |
---|
2205 | double xNow = 0.; |
---|
2206 | double wtMassNow, mT34Min, wtBW4Now, beta34Now; |
---|
2207 | |
---|
2208 | // Step through increasing x values; gives m4 unambiguously. |
---|
2209 | do { |
---|
2210 | xNow += xStep; |
---|
2211 | wtMassNow = 0.; |
---|
2212 | m4 = mHatMax - m3 - xNow * mWidth[4]; |
---|
2213 | |
---|
2214 | // Check that inside phase space limit set by pTmin. |
---|
2215 | mT34Min = mT3Min + sqrt(m4*m4 + pT2HatMin); |
---|
2216 | if (mT34Min < mHatMax) { |
---|
2217 | |
---|
2218 | // Breit-Wigner and beta factor give total weight. |
---|
2219 | wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) ); |
---|
2220 | beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4) |
---|
2221 | - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax); |
---|
2222 | wtMassNow = wtBW4Now * beta34Now; |
---|
2223 | |
---|
2224 | // Store new maximum, if any. |
---|
2225 | if (wtMassNow > wtMassMax) { |
---|
2226 | foundNonZero = true; |
---|
2227 | wtMassMax = wtMassNow; |
---|
2228 | m4WtMax = m4; |
---|
2229 | } |
---|
2230 | } |
---|
2231 | |
---|
2232 | // Continue stepping if increasing trend and more x range available. |
---|
2233 | } while ( (!foundNonZero || wtMassNow > wtMassMax) |
---|
2234 | && xNow < xMax - xStep); |
---|
2235 | |
---|
2236 | // Restore best value for subsequent maximization. |
---|
2237 | m4 = m4WtMax; |
---|
2238 | return foundNonZero; |
---|
2239 | |
---|
2240 | } |
---|
2241 | |
---|
2242 | //========================================================================== |
---|
2243 | |
---|
2244 | // PhaseSpace2to2elastic class. |
---|
2245 | // 2 -> 2 kinematics set up for elastic scattering. |
---|
2246 | |
---|
2247 | //-------------------------------------------------------------------------- |
---|
2248 | |
---|
2249 | // Constants: could be changed here if desired, but normally should not. |
---|
2250 | // These are of technical nature, as described for each. |
---|
2251 | |
---|
2252 | // Maximum positive/negative argument for exponentiation. |
---|
2253 | const double PhaseSpace2to2elastic::EXPMAX = 50.; |
---|
2254 | |
---|
2255 | // Conversion coefficients = 1/(16pi) * (mb <-> GeV^2). |
---|
2256 | const double PhaseSpace2to2elastic::CONVERTEL = 0.0510925; |
---|
2257 | |
---|
2258 | //-------------------------------------------------------------------------- |
---|
2259 | |
---|
2260 | // Form of phase space sampling already fixed, so no optimization. |
---|
2261 | // However, need to read out relevant parameters from SigmaTotal. |
---|
2262 | |
---|
2263 | bool PhaseSpace2to2elastic::setupSampling() { |
---|
2264 | |
---|
2265 | // Find maximum = value of cross section. |
---|
2266 | sigmaNw = sigmaProcessPtr->sigmaHatWrap(); |
---|
2267 | sigmaMx = sigmaNw; |
---|
2268 | |
---|
2269 | // Squared and outgoing masses of particles. |
---|
2270 | s1 = mA * mA; |
---|
2271 | s2 = mB * mB; |
---|
2272 | m3 = mA; |
---|
2273 | m4 = mB; |
---|
2274 | |
---|
2275 | // Elastic slope. |
---|
2276 | bSlope = sigmaTotPtr->bSlopeEl(); |
---|
2277 | |
---|
2278 | // Determine maximum possible t range. |
---|
2279 | lambda12S = pow2(s - s1 - s2) - 4. * s1 * s2 ; |
---|
2280 | tLow = - lambda12S / s; |
---|
2281 | tUpp = 0; |
---|
2282 | |
---|
2283 | // Production model with Coulomb corrections need more parameters. |
---|
2284 | useCoulomb = settingsPtr->flag("SigmaTotal:setOwn") |
---|
2285 | && settingsPtr->flag("SigmaElastic:setOwn"); |
---|
2286 | if (useCoulomb) { |
---|
2287 | sigmaTot = sigmaTotPtr->sigmaTot(); |
---|
2288 | rho = settingsPtr->parm("SigmaElastic:rho"); |
---|
2289 | lambda = settingsPtr->parm("SigmaElastic:lambda"); |
---|
2290 | tAbsMin = settingsPtr->parm("SigmaElastic:tAbsMin"); |
---|
2291 | phaseCst = settingsPtr->parm("SigmaElastic:phaseConst"); |
---|
2292 | alphaEM0 = settingsPtr->parm("StandardModel:alphaEM0"); |
---|
2293 | |
---|
2294 | // Relative rate of nuclear and Coulombic parts in trials. |
---|
2295 | tUpp = -tAbsMin; |
---|
2296 | sigmaNuc = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho) / bSlope |
---|
2297 | * exp(-bSlope * tAbsMin); |
---|
2298 | sigmaCou = (useCoulomb) ? |
---|
2299 | pow2(alphaEM0) / (4. * CONVERTEL * tAbsMin) : 0.; |
---|
2300 | signCou = (idA == idB) ? 1. : -1.; |
---|
2301 | |
---|
2302 | // Dummy values. |
---|
2303 | } else { |
---|
2304 | sigmaNuc = sigmaNw; |
---|
2305 | sigmaCou = 0.; |
---|
2306 | } |
---|
2307 | |
---|
2308 | // Calculate coefficient of generation. |
---|
2309 | tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.; |
---|
2310 | |
---|
2311 | return true; |
---|
2312 | |
---|
2313 | } |
---|
2314 | |
---|
2315 | //-------------------------------------------------------------------------- |
---|
2316 | |
---|
2317 | // Select a trial kinematics phase space point. Perform full |
---|
2318 | // Monte Carlo acceptance/rejection at this stage. |
---|
2319 | |
---|
2320 | bool PhaseSpace2to2elastic::trialKin( bool, bool ) { |
---|
2321 | |
---|
2322 | // Select t according to exp(bSlope*t). |
---|
2323 | if (!useCoulomb || sigmaNuc > rndmPtr->flat() * (sigmaNuc + sigmaCou)) |
---|
2324 | tH = tUpp + log(1. + tAux * rndmPtr->flat()) / bSlope; |
---|
2325 | |
---|
2326 | // Select t according to 1/t^2. |
---|
2327 | else tH = tLow * tUpp / (tUpp + rndmPtr->flat() * (tLow - tUpp)); |
---|
2328 | |
---|
2329 | // Correction factor for ratio full/simulated. |
---|
2330 | if (useCoulomb) { |
---|
2331 | double sigmaN = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho) |
---|
2332 | * exp(bSlope * tH); |
---|
2333 | double alpEM = couplingsPtr->alphaEM(-tH); |
---|
2334 | double sigmaC = pow2(alpEM) / (4. * CONVERTEL * tH*tH); |
---|
2335 | double sigmaGen = 2. * (sigmaN + sigmaC); |
---|
2336 | double form2 = pow4(lambda/(lambda - tH)); |
---|
2337 | double phase = signCou * alpEM |
---|
2338 | * (-phaseCst - log(-0.5 * bSlope * tH)); |
---|
2339 | double sigmaCor = sigmaN + pow2(form2) * sigmaC |
---|
2340 | - signCou * alpEM * sigmaTot * (form2 / (-tH)) |
---|
2341 | * exp(0.5 * bSlope * tH) * (rho * cos(phase) + sin(phase)); |
---|
2342 | sigmaNw = sigmaMx * sigmaCor / sigmaGen; |
---|
2343 | } |
---|
2344 | |
---|
2345 | // Careful reconstruction of scattering angle. |
---|
2346 | double tRat = s * tH / lambda12S; |
---|
2347 | double cosTheta = min(1., max(-1., 1. + 2. * tRat ) ); |
---|
2348 | double sinTheta = 2. * sqrtpos( -tRat * (1. + tRat) ); |
---|
2349 | theta = asin( min(1., sinTheta)); |
---|
2350 | if (cosTheta < 0.) theta = M_PI - theta; |
---|
2351 | |
---|
2352 | return true; |
---|
2353 | |
---|
2354 | } |
---|
2355 | |
---|
2356 | //-------------------------------------------------------------------------- |
---|
2357 | |
---|
2358 | // Construct the four-vector kinematics from the trial values. |
---|
2359 | |
---|
2360 | bool PhaseSpace2to2elastic::finalKin() { |
---|
2361 | |
---|
2362 | // Particle masses. |
---|
2363 | mH[1] = mA; |
---|
2364 | mH[2] = mB; |
---|
2365 | mH[3] = m3; |
---|
2366 | mH[4] = m4; |
---|
2367 | |
---|
2368 | // Incoming particles along beam axes. |
---|
2369 | pAbs = 0.5 * sqrtpos(lambda12S) / eCM; |
---|
2370 | pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM); |
---|
2371 | pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM); |
---|
2372 | |
---|
2373 | // Outgoing particles initially along beam axes. |
---|
2374 | pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM); |
---|
2375 | pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM); |
---|
2376 | |
---|
2377 | // Then rotate them |
---|
2378 | phi = 2. * M_PI * rndmPtr->flat(); |
---|
2379 | pH[3].rot( theta, phi); |
---|
2380 | pH[4].rot( theta, phi); |
---|
2381 | |
---|
2382 | // Set some further info for completeness. |
---|
2383 | x1H = 1.; |
---|
2384 | x2H = 1.; |
---|
2385 | sH = s; |
---|
2386 | uH = 2. * (s1 + s2) - sH - tH; |
---|
2387 | mHat = eCM; |
---|
2388 | p2Abs = pAbs * pAbs; |
---|
2389 | betaZ = 0.; |
---|
2390 | pTH = pAbs * sin(theta); |
---|
2391 | |
---|
2392 | // Done. |
---|
2393 | return true; |
---|
2394 | |
---|
2395 | } |
---|
2396 | |
---|
2397 | //========================================================================== |
---|
2398 | |
---|
2399 | // PhaseSpace2to2diffractive class. |
---|
2400 | // 2 -> 2 kinematics set up for diffractive scattering. |
---|
2401 | |
---|
2402 | //-------------------------------------------------------------------------- |
---|
2403 | |
---|
2404 | // Constants: could be changed here if desired, but normally should not. |
---|
2405 | // These are of technical nature, as described for each. |
---|
2406 | |
---|
2407 | // Number of tries to find acceptable (m^2, t) set. |
---|
2408 | const int PhaseSpace2to2diffractive::NTRY = 500; |
---|
2409 | |
---|
2410 | // Maximum positive/negative argument for exponentiation. |
---|
2411 | const double PhaseSpace2to2diffractive::EXPMAX = 50.; |
---|
2412 | |
---|
2413 | // Safety margin so sum of diffractive masses not too close to eCM. |
---|
2414 | const double PhaseSpace2to2diffractive::DIFFMASSMARGIN = 0.2; |
---|
2415 | |
---|
2416 | //-------------------------------------------------------------------------- |
---|
2417 | |
---|
2418 | // Form of phase space sampling already fixed, so no optimization. |
---|
2419 | // However, need to read out relevant parameters from SigmaTotal. |
---|
2420 | |
---|
2421 | bool PhaseSpace2to2diffractive::setupSampling() { |
---|
2422 | |
---|
2423 | // Pomeron flux parametrization, and parameters of some options. |
---|
2424 | PomFlux = settingsPtr->mode("Diffraction:PomFlux"); |
---|
2425 | epsilonPF = settingsPtr->parm("Diffraction:PomFluxEpsilon"); |
---|
2426 | alphaPrimePF = settingsPtr->parm("Diffraction:PomFluxAlphaPrime"); |
---|
2427 | |
---|
2428 | // Find maximum = value of cross section. |
---|
2429 | sigmaNw = sigmaProcessPtr->sigmaHatWrap(); |
---|
2430 | sigmaMx = sigmaNw; |
---|
2431 | |
---|
2432 | // Masses of particles and minimal masses of diffractive states. |
---|
2433 | m3ElDiff = (isDiffA) ? sigmaTotPtr->mMinXB() : mA; |
---|
2434 | m4ElDiff = (isDiffB) ? sigmaTotPtr->mMinAX() : mB; |
---|
2435 | s1 = mA * mA; |
---|
2436 | s2 = mB * mB; |
---|
2437 | s3 = pow2( m3ElDiff); |
---|
2438 | s4 = pow2( m4ElDiff); |
---|
2439 | |
---|
2440 | // Determine maximum possible t range and coefficient of generation. |
---|
2441 | lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 ); |
---|
2442 | lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 ); |
---|
2443 | double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s; |
---|
2444 | double tempB = lambda12 * lambda34 / s; |
---|
2445 | double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3) |
---|
2446 | * (s1 * s4 - s2 * s3) / s; |
---|
2447 | tLow = -0.5 * (tempA + tempB); |
---|
2448 | tUpp = tempC / tLow; |
---|
2449 | |
---|
2450 | // Default for all parametrization-specific parameters. |
---|
2451 | cRes = sResXB = sResAX = sProton = bMin = bSlope = bSlope1 = bSlope2 |
---|
2452 | = probSlope1 = xIntPF = xtCorPF = mp24DL = coefDL = tAux |
---|
2453 | = tAux1 = tAux2 = 0.; |
---|
2454 | |
---|
2455 | // Schuler&Sjostrand: parameters of low-mass-resonance enhancement. |
---|
2456 | if (PomFlux == 1) { |
---|
2457 | cRes = sigmaTotPtr->cRes(); |
---|
2458 | sResXB = pow2( sigmaTotPtr->mResXB()); |
---|
2459 | sResAX = pow2( sigmaTotPtr->mResAX()); |
---|
2460 | sProton = sigmaTotPtr->sProton(); |
---|
2461 | |
---|
2462 | // Schuler&Sjostrand: lower limit diffractive slope. |
---|
2463 | if (!isDiffB) bMin = sigmaTotPtr->bMinSlopeXB(); |
---|
2464 | else if (!isDiffA) bMin = sigmaTotPtr->bMinSlopeAX(); |
---|
2465 | else bMin = sigmaTotPtr->bMinSlopeXX(); |
---|
2466 | tAux = exp( max(-EXPMAX, bMin * (tLow - tUpp)) ) - 1.; |
---|
2467 | |
---|
2468 | // Bruni&Ingelman: relative weight of two diffractive slopes. |
---|
2469 | } else if (PomFlux == 2) { |
---|
2470 | bSlope1 = 8.0; |
---|
2471 | probSlope1 = 6.38 * ( exp(max(-EXPMAX, bSlope1 * tUpp)) |
---|
2472 | - exp(max(-EXPMAX, bSlope1 * tLow)) ) / bSlope1; |
---|
2473 | bSlope2 = 3.0; |
---|
2474 | double pS2 = 0.424 * ( exp(max(-EXPMAX, bSlope2 * tUpp)) |
---|
2475 | - exp(max(-EXPMAX, bSlope2 * tLow)) ) / bSlope2; |
---|
2476 | probSlope1 /= probSlope1 + pS2; |
---|
2477 | tAux1 = exp( max(-EXPMAX, bSlope1 * (tLow - tUpp)) ) - 1.; |
---|
2478 | tAux2 = exp( max(-EXPMAX, bSlope2 * (tLow - tUpp)) ) - 1.; |
---|
2479 | |
---|
2480 | // Streng&Berger (RapGap): diffractive slope, power of mass spectrum. |
---|
2481 | } else if (PomFlux == 3) { |
---|
2482 | bSlope = 4.7; |
---|
2483 | double xPowPF = 1. - 2. * (1. + epsilonPF); |
---|
2484 | xIntPF = 2. * (1. + xPowPF); |
---|
2485 | xtCorPF = 2. * alphaPrimePF; |
---|
2486 | tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.; |
---|
2487 | |
---|
2488 | // Donnachie&Landshoff (RapGap): power of mass spectrum. |
---|
2489 | } else if (PomFlux == 4) { |
---|
2490 | mp24DL = 4. * pow2(particleDataPtr->m0(2212)); |
---|
2491 | double xPowPF = 1. - 2. * (1. + epsilonPF); |
---|
2492 | xIntPF = 2. * (1. + xPowPF); |
---|
2493 | xtCorPF = 2. * alphaPrimePF; |
---|
2494 | // Upper estimate of t dependence, for preliminary choice. |
---|
2495 | coefDL = 0.85; |
---|
2496 | tAux1 = 1. / pow3(1. - coefDL * tLow); |
---|
2497 | tAux2 = 1. / pow3(1. - coefDL * tUpp); |
---|
2498 | |
---|
2499 | // MBR model. |
---|
2500 | } else if (PomFlux == 5) { |
---|
2501 | eps = settingsPtr->parm("Diffraction:MBRepsilon"); |
---|
2502 | alph = settingsPtr->parm("Diffraction:MBRalpha"); |
---|
2503 | alph2 = alph * alph; |
---|
2504 | m2min = settingsPtr->parm("Diffraction:MBRm2Min"); |
---|
2505 | dyminSD = settingsPtr->parm("Diffraction:MBRdyminSD"); |
---|
2506 | dyminDD = settingsPtr->parm("Diffraction:MBRdyminDD"); |
---|
2507 | dyminSigSD = settingsPtr->parm("Diffraction:MBRdyminSigSD"); |
---|
2508 | dyminSigDD = settingsPtr->parm("Diffraction:MBRdyminSigDD"); |
---|
2509 | |
---|
2510 | // Max f(dy) for Von Neumann method, from SigmaTot. |
---|
2511 | sdpmax= sigmaTotPtr->sdpMax(); |
---|
2512 | ddpmax= sigmaTotPtr->ddpMax(); |
---|
2513 | } |
---|
2514 | |
---|
2515 | // Done. |
---|
2516 | return true; |
---|
2517 | |
---|
2518 | } |
---|
2519 | |
---|
2520 | //-------------------------------------------------------------------------- |
---|
2521 | |
---|
2522 | // Select a trial kinematics phase space point. Perform full |
---|
2523 | // Monte Carlo acceptance/rejection at this stage. |
---|
2524 | |
---|
2525 | bool PhaseSpace2to2diffractive::trialKin( bool, bool ) { |
---|
2526 | |
---|
2527 | // Loop over attempts to set up masses and t consistently. |
---|
2528 | for (int loop = 0; ; ++loop) { |
---|
2529 | if (loop == NTRY) { |
---|
2530 | infoPtr->errorMsg("Error in PhaseSpace2to2diffractive::trialKin: " |
---|
2531 | " quit after repeated tries"); |
---|
2532 | return false; |
---|
2533 | } |
---|
2534 | |
---|
2535 | // Schuler and Sjostrand: |
---|
2536 | if (PomFlux == 1) { |
---|
2537 | |
---|
2538 | // Select diffractive mass(es) according to dm^2/m^2. |
---|
2539 | m3 = (isDiffA) ? m3ElDiff * pow( max(mA, eCM - m4ElDiff) / m3ElDiff, |
---|
2540 | rndmPtr->flat()) : m3ElDiff; |
---|
2541 | m4 = (isDiffB) ? m4ElDiff * pow( max(mB, eCM - m3ElDiff) / m4ElDiff, |
---|
2542 | rndmPtr->flat()) : m4ElDiff; |
---|
2543 | if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue; |
---|
2544 | s3 = m3 * m3; |
---|
2545 | s4 = m4 * m4; |
---|
2546 | |
---|
2547 | // Additional mass factors, including resonance enhancement. |
---|
2548 | if (isDiffA && !isDiffB) { |
---|
2549 | double facXB = (1. - s3 / s) |
---|
2550 | * (1. + cRes * sResXB / (sResXB + s3)); |
---|
2551 | if (facXB < rndmPtr->flat() * (1. + cRes)) continue; |
---|
2552 | } else if (isDiffB && !isDiffA) { |
---|
2553 | double facAX = (1. - s4 / s) |
---|
2554 | * (1. + cRes * sResAX / (sResAX + s4)); |
---|
2555 | if (facAX < rndmPtr->flat() * (1. + cRes)) continue; |
---|
2556 | } else { |
---|
2557 | double facXX = (1. - pow2(m3 + m4) / s) |
---|
2558 | * (s * sProton / (s * sProton + s3 * s4)) |
---|
2559 | * (1. + cRes * sResXB / (sResXB + s3)) |
---|
2560 | * (1. + cRes * sResAX / (sResAX + s4)); |
---|
2561 | if (facXX < rndmPtr->flat() * pow2(1. + cRes)) continue; |
---|
2562 | } |
---|
2563 | |
---|
2564 | // Select t according to exp(bMin*t) and correct to right slope. |
---|
2565 | tH = tUpp + log(1. + tAux * rndmPtr->flat()) / bMin; |
---|
2566 | double bDiff = 0.; |
---|
2567 | if (isDiffA && !isDiffB) bDiff = sigmaTotPtr->bSlopeXB(s3) - bMin; |
---|
2568 | else if (!isDiffA) bDiff = sigmaTotPtr->bSlopeAX(s4) - bMin; |
---|
2569 | else bDiff = sigmaTotPtr->bSlopeXX(s3, s4) - bMin; |
---|
2570 | bDiff = max(0., bDiff); |
---|
2571 | if (exp( max(-EXPMAX, bDiff * (tH - tUpp)) ) < rndmPtr->flat()) continue; |
---|
2572 | |
---|
2573 | // Bruni and Ingelman: |
---|
2574 | } else if (PomFlux == 2) { |
---|
2575 | |
---|
2576 | // Select diffractive mass(es) according to dm^2/m^2. |
---|
2577 | m3 = (isDiffA) ? m3ElDiff * pow( max(mA, eCM - m4ElDiff) / m3ElDiff, |
---|
2578 | rndmPtr->flat()) : m3ElDiff; |
---|
2579 | m4 = (isDiffB) ? m4ElDiff * pow( max(mB, eCM - m3ElDiff) / m4ElDiff, |
---|
2580 | rndmPtr->flat()) : m4ElDiff; |
---|
2581 | if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue; |
---|
2582 | s3 = m3 * m3; |
---|
2583 | s4 = m4 * m4; |
---|
2584 | |
---|
2585 | // Select t according to exp(bSlope*t) with two possible slopes. |
---|
2586 | tH = (rndmPtr->flat() < probSlope1) |
---|
2587 | ? tUpp + log(1. + tAux1 * rndmPtr->flat()) / bSlope1 |
---|
2588 | : tUpp + log(1. + tAux2 * rndmPtr->flat()) / bSlope2; |
---|
2589 | |
---|
2590 | // Streng and Berger et al. (RapGap): |
---|
2591 | } else if (PomFlux == 3) { |
---|
2592 | |
---|
2593 | // Select diffractive mass(es) according to dm^2/(m^2)^(1 + 2 epsilon). |
---|
2594 | m3 = m3ElDiff; |
---|
2595 | m4 = m4ElDiff; |
---|
2596 | if (isDiffA) { |
---|
2597 | double s3MinPow = pow( m3ElDiff, xIntPF ); |
---|
2598 | double s3MaxPow = pow( max(mA, eCM - m4ElDiff), xIntPF ); |
---|
2599 | m3 = pow( s3MinPow + rndmPtr->flat() * (s3MaxPow - s3MinPow), |
---|
2600 | 1. / xIntPF ); |
---|
2601 | } |
---|
2602 | if (isDiffB) { |
---|
2603 | double s4MinPow = pow( m4ElDiff, xIntPF ); |
---|
2604 | double s4MaxPow = pow( max(mB, eCM - m3ElDiff), xIntPF ); |
---|
2605 | m4 = pow( s4MinPow + rndmPtr->flat() * (s4MaxPow - s4MinPow), |
---|
2606 | 1. / xIntPF ); |
---|
2607 | } |
---|
2608 | if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue; |
---|
2609 | s3 = m3 * m3; |
---|
2610 | s4 = m4 * m4; |
---|
2611 | |
---|
2612 | // Select t according to exponential and weigh by x_P^(2 alpha' |t|). |
---|
2613 | tH = tUpp + log(1. + tAux * rndmPtr->flat()) / bSlope; |
---|
2614 | if ( isDiffA && pow( s3 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() ) |
---|
2615 | continue; |
---|
2616 | if ( isDiffB && pow( s4 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() ) |
---|
2617 | continue; |
---|
2618 | |
---|
2619 | // Donnachie and Landshoff (RapGap): |
---|
2620 | } else if (PomFlux == 4) { |
---|
2621 | |
---|
2622 | // Select diffractive mass(es) according to dm^2/(m^2)^(1 + 2 epsilon). |
---|
2623 | m3 = m3ElDiff; |
---|
2624 | m4 = m4ElDiff; |
---|
2625 | if (isDiffA) { |
---|
2626 | double s3MinPow = pow( m3ElDiff, xIntPF ); |
---|
2627 | double s3MaxPow = pow( max(mA, eCM - m4ElDiff), xIntPF ); |
---|
2628 | m3 = pow( s3MinPow + rndmPtr->flat() * (s3MaxPow - s3MinPow), |
---|
2629 | 1. / xIntPF ); |
---|
2630 | } |
---|
2631 | if (isDiffB) { |
---|
2632 | double s4MinPow = pow( m4ElDiff, xIntPF ); |
---|
2633 | double s4MaxPow = pow( max(mB, eCM - m3ElDiff), xIntPF ); |
---|
2634 | m4 = pow( s4MinPow + rndmPtr->flat() * (s4MaxPow - s4MinPow), |
---|
2635 | 1. / xIntPF ); |
---|
2636 | } |
---|
2637 | if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue; |
---|
2638 | s3 = m3 * m3; |
---|
2639 | s4 = m4 * m4; |
---|
2640 | |
---|
2641 | // Select t according to power and weigh by x_P^(2 alpha' |t|). |
---|
2642 | tH = - (1. / pow( tAux1 + rndmPtr->flat() * (tAux2 - tAux1), 1./3.) |
---|
2643 | - 1.) / coefDL; |
---|
2644 | double wDL = pow2( (mp24DL - 2.8 * tH) / (mp24DL - tH) ) |
---|
2645 | / pow4( 1. - tH / 0.7); |
---|
2646 | double wMX = 1. / pow4( 1. - coefDL * tH); |
---|
2647 | if (wDL < rndmPtr->flat() * wMX) continue; |
---|
2648 | if ( isDiffA && pow( s3 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() ) |
---|
2649 | continue; |
---|
2650 | if ( isDiffB && pow( s4 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() ) |
---|
2651 | continue; |
---|
2652 | |
---|
2653 | // MBR model: |
---|
2654 | } else if (PomFlux == 5) { |
---|
2655 | m3 = mA; |
---|
2656 | m4 = mB; |
---|
2657 | double xi, P, yRnd, dy; |
---|
2658 | |
---|
2659 | // MBR double diffractive. |
---|
2660 | if (isDiffA && isDiffB) { |
---|
2661 | dymin0 = 0.; |
---|
2662 | dymax = log(s/pow2(m2min)); |
---|
2663 | |
---|
2664 | // Von Neumann method to generate dy, uses ddpmax from SigmaTot. |
---|
2665 | do { |
---|
2666 | dy = dymin0 + (dymax - dymin0) * rndmPtr->flat(); |
---|
2667 | P = (dymax - dy) * exp(eps*dy) * ( exp(-2. * alph * dy * exp(-dy)) |
---|
2668 | - exp(-2. * alph * dy * exp(dy)) ) / dy; |
---|
2669 | // Suppress smaller gap, smooth transition to non-diffractive. |
---|
2670 | P *= 0.5 * (1 + erf( ( dy - dyminDD) / dyminSigDD ) ); |
---|
2671 | if (P > ddpmax) { |
---|
2672 | ostringstream osWarn; |
---|
2673 | osWarn << "ddpmax = " << scientific << setprecision(3) |
---|
2674 | << ddpmax << " " << P << " " << dy << endl; |
---|
2675 | infoPtr->errorMsg("Warning in PhaseSpace2to2diffractive::" |
---|
2676 | "trialKin for double diffraction:", osWarn.str()); |
---|
2677 | } |
---|
2678 | yRnd = ddpmax * rndmPtr->flat(); |
---|
2679 | } while (yRnd > P); |
---|
2680 | |
---|
2681 | double y0max = (dymax - dy)/2.; |
---|
2682 | double y0min = -y0max; |
---|
2683 | double y0 = y0min + (y0max - y0min) * rndmPtr->flat(); |
---|
2684 | am1 = sqrt( eCM * exp( -y0 - dy/2. ) ); |
---|
2685 | am2 = sqrt( eCM * exp( y0 - dy/2. ) ); |
---|
2686 | |
---|
2687 | // Generate 4-momentum transfer, t from exp. |
---|
2688 | double b = 2. * alph * dy; |
---|
2689 | tUpp = -exp( -dy ); |
---|
2690 | tLow = -exp( dy ); |
---|
2691 | tAux = exp( b * (tLow - tUpp) ) - 1.; |
---|
2692 | t = tUpp + log(1. + tAux * rndmPtr->flat()) / b; |
---|
2693 | m3 = am1; |
---|
2694 | m4 = am2; |
---|
2695 | tH = t; |
---|
2696 | |
---|
2697 | // MBR single diffractive. |
---|
2698 | } else if (isDiffA || isDiffB) { |
---|
2699 | dymin0 = 0.; |
---|
2700 | dymax = log(s/m2min); |
---|
2701 | |
---|
2702 | // Von Neumann method to generate dy, uses sdpmax from SigmaTot. |
---|
2703 | do { |
---|
2704 | dy = dymin0 + (dymax - dymin0) * rndmPtr->flat(); |
---|
2705 | P = exp(eps * dy) * ( (FFA1 / (FFB1 + 2. * alph * dy) ) |
---|
2706 | + (FFA2 / (FFB2 + 2. * alph * dy) ) ); |
---|
2707 | // Suppress smaller gap. |
---|
2708 | P *= 0.5 * (1. + erf( (dy - dyminSD) / dyminSigSD) ); |
---|
2709 | if (P > sdpmax) { |
---|
2710 | ostringstream osWarn; |
---|
2711 | osWarn << "sdpmax = " << scientific << setprecision(3) |
---|
2712 | << sdpmax << " " << P << " " << dy << endl; |
---|
2713 | infoPtr->errorMsg("Warning in PhaseSpace2to2diffractive::" |
---|
2714 | "trialKin for single diffraction:", osWarn.str()); |
---|
2715 | } |
---|
2716 | yRnd = sdpmax * rndmPtr->flat(); |
---|
2717 | } while (yRnd > P); |
---|
2718 | xi = exp( -dy ); |
---|
2719 | amx = sqrt( xi * s ); |
---|
2720 | |
---|
2721 | // Generate 4-momentum transfer, t. First exponent, then FF*exp. |
---|
2722 | double tmin = -s1 * xi * xi / (1 - xi); |
---|
2723 | do { |
---|
2724 | t = tmin + log(1. - rndmPtr->flat()); |
---|
2725 | double pFF = (4. * s1 - 2.8 * t) / ( (4. * s1 - t) |
---|
2726 | * pow2(1. - t / 0.71) ); |
---|
2727 | P = pow2(pFF) * exp(2. * alph * dy * t); |
---|
2728 | yRnd = exp(t) * rndmPtr->flat(); |
---|
2729 | } while (yRnd > P); |
---|
2730 | if(isDiffA) m3 = amx; |
---|
2731 | if(isDiffB) m4 = amx; |
---|
2732 | tH = t; |
---|
2733 | } |
---|
2734 | |
---|
2735 | // End of MBR model code. |
---|
2736 | s3 = m3 * m3; |
---|
2737 | s4 = m4 * m4; |
---|
2738 | } |
---|
2739 | |
---|
2740 | // Check whether m^2 and t choices are consistent. |
---|
2741 | lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 ); |
---|
2742 | double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s; |
---|
2743 | double tempB = lambda12 * lambda34 / s; |
---|
2744 | double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3) |
---|
2745 | * (s1 * s4 - s2 * s3) / s; |
---|
2746 | double tLowNow = -0.5 * (tempA + tempB); |
---|
2747 | double tUppNow = tempC / tLowNow; |
---|
2748 | if (tH < tLowNow || tH > tUppNow) continue; |
---|
2749 | |
---|
2750 | // Careful reconstruction of scattering angle. |
---|
2751 | double cosTheta = min(1., max(-1., (tempA + 2. * tH) / tempB)); |
---|
2752 | double sinTheta = 2. * sqrtpos( -(tempC + tempA * tH + tH * tH) ) |
---|
2753 | / tempB; |
---|
2754 | theta = asin( min(1., sinTheta)); |
---|
2755 | if (cosTheta < 0.) theta = M_PI - theta; |
---|
2756 | |
---|
2757 | // Found acceptable kinematics, so no more looping. Done |
---|
2758 | break; |
---|
2759 | } |
---|
2760 | return true; |
---|
2761 | |
---|
2762 | } |
---|
2763 | |
---|
2764 | //-------------------------------------------------------------------------- |
---|
2765 | |
---|
2766 | // Construct the four-vector kinematics from the trial values. |
---|
2767 | |
---|
2768 | bool PhaseSpace2to2diffractive::finalKin() { |
---|
2769 | |
---|
2770 | // Particle masses; incoming always on mass shell. |
---|
2771 | mH[1] = mA; |
---|
2772 | mH[2] = mB; |
---|
2773 | mH[3] = m3; |
---|
2774 | mH[4] = m4; |
---|
2775 | |
---|
2776 | // Incoming particles along beam axes. |
---|
2777 | pAbs = 0.5 * lambda12 / eCM; |
---|
2778 | pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM); |
---|
2779 | pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM); |
---|
2780 | |
---|
2781 | // Outgoing particles initially along beam axes. |
---|
2782 | pAbs = 0.5 * lambda34 / eCM; |
---|
2783 | pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s3 - s4) / eCM); |
---|
2784 | pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s4 - s3) / eCM); |
---|
2785 | |
---|
2786 | // Then rotate them |
---|
2787 | phi = 2. * M_PI * rndmPtr->flat(); |
---|
2788 | pH[3].rot( theta, phi); |
---|
2789 | pH[4].rot( theta, phi); |
---|
2790 | |
---|
2791 | // Set some further info for completeness (but Info can be for subprocess). |
---|
2792 | x1H = 1.; |
---|
2793 | x2H = 1.; |
---|
2794 | sH = s; |
---|
2795 | uH = s1 + s2 + s3 + s4 - sH - tH; |
---|
2796 | mHat = eCM; |
---|
2797 | p2Abs = pAbs * pAbs; |
---|
2798 | betaZ = 0.; |
---|
2799 | pTH = pAbs * sin(theta); |
---|
2800 | |
---|
2801 | // Done. |
---|
2802 | return true; |
---|
2803 | |
---|
2804 | } |
---|
2805 | |
---|
2806 | //========================================================================== |
---|
2807 | |
---|
2808 | // PhaseSpace2to3diffractive class. |
---|
2809 | // 2 -> 3 kinematics set up for central diffractive scattering. |
---|
2810 | |
---|
2811 | //-------------------------------------------------------------------------- |
---|
2812 | |
---|
2813 | // Constants: could be changed here if desired, but normally should not. |
---|
2814 | // These are of technical nature, as described for each. |
---|
2815 | |
---|
2816 | // Number of tries to find acceptable (m^2, t1, t2) set. |
---|
2817 | const int PhaseSpace2to3diffractive::NTRY = 500; |
---|
2818 | const int PhaseSpace2to3diffractive::NINTEG2 = 40; |
---|
2819 | |
---|
2820 | // Maximum positive/negative argument for exponentiation. |
---|
2821 | const double PhaseSpace2to3diffractive::EXPMAX = 50.; |
---|
2822 | |
---|
2823 | // Minimal mass of central diffractive system. |
---|
2824 | const double PhaseSpace2to3diffractive::DIFFMASSMIN = 0.8; |
---|
2825 | |
---|
2826 | // Safety margin so sum of diffractive masses not too close to eCM. |
---|
2827 | const double PhaseSpace2to3diffractive::DIFFMASSMARGIN = 0.2; |
---|
2828 | |
---|
2829 | //-------------------------------------------------------------------------- |
---|
2830 | |
---|
2831 | // Set up for phase space sampling. |
---|
2832 | |
---|
2833 | bool PhaseSpace2to3diffractive::setupSampling() { |
---|
2834 | |
---|
2835 | // Pomeron flux parametrization, and parameters of some options. |
---|
2836 | PomFlux = settingsPtr->mode("Diffraction:PomFlux"); |
---|
2837 | epsilonPF = settingsPtr->parm("Diffraction:PomFluxEpsilon"); |
---|
2838 | alphaPrimePF = settingsPtr->parm("Diffraction:PomFluxAlphaPrime"); |
---|
2839 | |
---|
2840 | // Find maximum = value of cross section. |
---|
2841 | sigmaNw = sigmaProcessPtr->sigmaHatWrap(); |
---|
2842 | sigmaMx = sigmaNw; |
---|
2843 | |
---|
2844 | // Squared masses of particles and minimal mass of diffractive states. |
---|
2845 | s1 = mA * mA; |
---|
2846 | s2 = mB * mB; |
---|
2847 | m5min = sigmaTotPtr->mMinAXB(); |
---|
2848 | s5min = m5min * m5min; |
---|
2849 | |
---|
2850 | // Loop over two cases: s4 = (X + B)^2 and s3 = (A + X)^2. |
---|
2851 | for (int i = 0; i < 2; ++i) { |
---|
2852 | s3 = (i == 0) ? s1 : pow2(mA + m5min); |
---|
2853 | s4 = (i == 0) ? pow2(mB + m5min) : s2; |
---|
2854 | |
---|
2855 | // Determine maximum possible t range and coefficient of generation. |
---|
2856 | double lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 ); |
---|
2857 | double lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 ); |
---|
2858 | double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s; |
---|
2859 | double tempB = lambda12 * lambda34 / s; |
---|
2860 | double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3) |
---|
2861 | * (s1 * s4 - s2 * s3) / s; |
---|
2862 | tLow[i] = -0.5 * (tempA + tempB); |
---|
2863 | tUpp[i] = tempC / tLow[i]; |
---|
2864 | } |
---|
2865 | s3 = s1; |
---|
2866 | s4 = s2; |
---|
2867 | |
---|
2868 | // Default for all parametrization-specific parameters. |
---|
2869 | bSlope1 = bSlope2 = bSlope = xIntPF = xIntInvPF = xtCorPF = mp24DL |
---|
2870 | = coefDL = 0.; |
---|
2871 | for (int i = 0; i < 2; ++i) |
---|
2872 | bMin[i] = tAux[i] = probSlope1[i] = tAux1[i] = tAux2[i] = 0.; |
---|
2873 | |
---|
2874 | // Schuler&Sjostrand: lower limit diffractive slope. |
---|
2875 | if (PomFlux == 1) { |
---|
2876 | bMin[0] = sigmaTotPtr->bMinSlopeAX(); |
---|
2877 | tAux[0] = exp( max(-EXPMAX, bMin[0] * (tLow[0] - tUpp[0])) ) - 1.; |
---|
2878 | bMin[1] = sigmaTotPtr->bMinSlopeXB(); |
---|
2879 | tAux[1] = exp( max(-EXPMAX, bMin[1] * (tLow[1] - tUpp[1])) ) - 1.; |
---|
2880 | |
---|
2881 | // Bruni&Ingelman: relative weight of two diffractive slopes. |
---|
2882 | } else if (PomFlux == 2) { |
---|
2883 | bSlope1 = 8.0; |
---|
2884 | bSlope2 = 3.0; |
---|
2885 | for (int i = 0; i < 2; ++i) { |
---|
2886 | probSlope1[i] = 6.38 * ( exp(max(-EXPMAX, bSlope1 * tUpp[i])) |
---|
2887 | - exp(max(-EXPMAX, bSlope1 * tLow[i])) ) / bSlope1; |
---|
2888 | double pS2 = 0.424 * ( exp(max(-EXPMAX, bSlope2 * tUpp[i])) |
---|
2889 | - exp(max(-EXPMAX, bSlope2 * tLow[i])) ) / bSlope2; |
---|
2890 | probSlope1[i] /= probSlope1[i] + pS2; |
---|
2891 | tAux1[i] = exp( max(-EXPMAX, bSlope1 * (tLow[i] - tUpp[i])) ) - 1.; |
---|
2892 | tAux2[i] = exp( max(-EXPMAX, bSlope2 * (tLow[i] - tUpp[i])) ) - 1.; |
---|
2893 | } |
---|
2894 | |
---|
2895 | // Streng&Berger (RapGap): diffractive slope, power of mass spectrum. |
---|
2896 | } else if (PomFlux == 3) { |
---|
2897 | bSlope = 4.7; |
---|
2898 | double xPowPF = 1. - 2. * (1. + epsilonPF); |
---|
2899 | xIntPF = 1. + xPowPF; |
---|
2900 | xIntInvPF = 1. / xIntPF; |
---|
2901 | xtCorPF = 2. * alphaPrimePF; |
---|
2902 | tAux[0] = exp( max(-EXPMAX, bSlope * (tLow[0] - tUpp[0])) ) - 1.; |
---|
2903 | tAux[1] = exp( max(-EXPMAX, bSlope * (tLow[1] - tUpp[1])) ) - 1.; |
---|
2904 | |
---|
2905 | // Donnachie&Landshoff (RapGap): power of mass spectrum. |
---|
2906 | } else if (PomFlux == 4) { |
---|
2907 | mp24DL = 4. * pow2(particleDataPtr->m0(2212)); |
---|
2908 | double xPowPF = 1. - 2. * (1. + epsilonPF); |
---|
2909 | xIntPF = 1. + xPowPF; |
---|
2910 | xIntInvPF = 1. / xIntPF; |
---|
2911 | xtCorPF = 2. * alphaPrimePF; |
---|
2912 | // Upper estimate of t dependence, for preliminary choice. |
---|
2913 | coefDL = 0.85; |
---|
2914 | tAux1[0] = 1. / pow3(1. - coefDL * tLow[0]); |
---|
2915 | tAux2[0] = 1. / pow3(1. - coefDL * tUpp[0]); |
---|
2916 | tAux1[1] = 1. / pow3(1. - coefDL * tLow[1]); |
---|
2917 | tAux2[1] = 1. / pow3(1. - coefDL * tUpp[1]); |
---|
2918 | |
---|
2919 | // Setup for the MBR model. |
---|
2920 | } else if (PomFlux == 5) { |
---|
2921 | epsMBR = settingsPtr->parm("Diffraction:MBRepsilon"); |
---|
2922 | alphMBR = settingsPtr->parm("Diffraction:MBRalpha"); |
---|
2923 | m2minMBR = settingsPtr->parm("Diffraction:MBRm2Min"); |
---|
2924 | dyminMBR = settingsPtr->parm("Diffraction:MBRdyminCD"); |
---|
2925 | dyminSigMBR = settingsPtr->parm("Diffraction:MBRdyminSigCD"); |
---|
2926 | dyminInvMBR = sqrt(2.) / dyminSigMBR; |
---|
2927 | // Max f(dy) for Von Neumann method, dpepmax from SigmaTot. |
---|
2928 | dpepmax = sigmaTotPtr->dpepMax(); |
---|
2929 | } |
---|
2930 | |
---|
2931 | // Done. |
---|
2932 | return true; |
---|
2933 | |
---|
2934 | } |
---|
2935 | |
---|
2936 | //-------------------------------------------------------------------------- |
---|
2937 | |
---|
2938 | // Select a trial kinematics phase space point. Perform full |
---|
2939 | // Monte Carlo acceptance/rejection at this stage. |
---|
2940 | |
---|
2941 | bool PhaseSpace2to3diffractive::trialKin( bool, bool ) { |
---|
2942 | |
---|
2943 | // Trivial kinematics of incoming hadrons. |
---|
2944 | double lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 ); |
---|
2945 | pAbs = 0.5 * lambda12 / eCM; |
---|
2946 | p1.p( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM); |
---|
2947 | p2.p( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM); |
---|
2948 | |
---|
2949 | // Loop over attempts to set up mass, t1, t2 consistently. |
---|
2950 | for (int loop = 0; ; ++loop) { |
---|
2951 | if (loop == NTRY) { |
---|
2952 | infoPtr->errorMsg("Error in PhaseSpace2to3diffractive::trialKin: " |
---|
2953 | " quit after repeated tries"); |
---|
2954 | return false; |
---|
2955 | } |
---|
2956 | double xi1 = 0.; |
---|
2957 | double xi2 = 0.; |
---|
2958 | double tVal[2] = { 0., 0.}; |
---|
2959 | |
---|
2960 | // Schuler and Sjostrand: |
---|
2961 | if (PomFlux == 1) { |
---|
2962 | |
---|
2963 | // Select mass according to dxi_1/xi_1 * dxi_2/xi_2 * (1 - m^2/s). |
---|
2964 | do { |
---|
2965 | xi1 = pow( s5min / s, rndmPtr->flat()); |
---|
2966 | xi2 = pow( s5min / s, rndmPtr->flat()); |
---|
2967 | s5 = xi1 * xi2 * s; |
---|
2968 | } while (s5 < s5min || xi1 * xi2 > rndmPtr->flat()); |
---|
2969 | if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue; |
---|
2970 | |
---|
2971 | // Select t according to exp(bMin*t) and correct to right slope. |
---|
2972 | bool tryAgain = false; |
---|
2973 | for (int i = 0; i < 2; ++i) { |
---|
2974 | tVal[i] = tUpp[i] + log(1. + tAux[i] * rndmPtr->flat()) / bMin[i]; |
---|
2975 | double bDiff = (i == 0) ? sigmaTotPtr->bSlopeAX(s2 + xi1 * s) |
---|
2976 | : sigmaTotPtr->bSlopeXB(s1 + xi2 * s); |
---|
2977 | bDiff = max(0., bDiff - bMin[i]); |
---|
2978 | if (exp( max(-EXPMAX, bDiff * (tVal[i] - tUpp[i])) ) |
---|
2979 | < rndmPtr->flat()) tryAgain = true; |
---|
2980 | } |
---|
2981 | if (tryAgain) continue; |
---|
2982 | |
---|
2983 | // Bruni and Ingelman: |
---|
2984 | } else if (PomFlux == 2) { |
---|
2985 | |
---|
2986 | // Select mass according to dxi_1/xi_1 * dxi_2/xi_2. |
---|
2987 | do { |
---|
2988 | xi1 = pow( s5min / s, rndmPtr->flat()); |
---|
2989 | xi2 = pow( s5min / s, rndmPtr->flat()); |
---|
2990 | s5 = xi1 * xi2 * s; |
---|
2991 | } while (s5 < s5min); |
---|
2992 | if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue; |
---|
2993 | |
---|
2994 | // Select t according to exp(bSlope*t) with two possible slopes. |
---|
2995 | for (int i = 0; i < 2; ++i) |
---|
2996 | tVal[i] = (rndmPtr->flat() < probSlope1[i]) |
---|
2997 | ? tUpp[i] + log(1. + tAux1[i] * rndmPtr->flat()) / bSlope1 |
---|
2998 | : tUpp[i] + log(1. + tAux2[i] * rndmPtr->flat()) / bSlope2; |
---|
2999 | |
---|
3000 | // Streng and Berger et al. (RapGap): |
---|
3001 | } else if (PomFlux == 3) { |
---|
3002 | |
---|
3003 | // Select mass by dxi_1 * dxi_2 / (xi_1 * xi_2)^(1 + 2 epsilon). |
---|
3004 | double sMinPow = pow( s5min / s, xIntPF); |
---|
3005 | do { |
---|
3006 | xi1 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF ); |
---|
3007 | xi2 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF ); |
---|
3008 | s5 = xi1 * xi2 * s; |
---|
3009 | } while (s5 < s5min); |
---|
3010 | if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue; |
---|
3011 | |
---|
3012 | // Select t according to exponential and weigh by x_P^(2 alpha' |t|). |
---|
3013 | bool tryAgain = false; |
---|
3014 | for (int i = 0; i < 2; ++i) { |
---|
3015 | tVal[i] = tUpp[i] + log(1. + tAux[i] * rndmPtr->flat()) / bSlope; |
---|
3016 | double xi = (i == 0) ? xi1 : xi2; |
---|
3017 | if ( pow( xi, xtCorPF * abs(tVal[i]) ) < rndmPtr->flat() ) |
---|
3018 | tryAgain = true; |
---|
3019 | } |
---|
3020 | if (tryAgain) continue; |
---|
3021 | |
---|
3022 | // Donnachie and Landshoff (RapGap): |
---|
3023 | } else if (PomFlux == 4) { |
---|
3024 | |
---|
3025 | // Select mass by dxi_1 * dxi_2 / (xi_1 * xi_2)^(1 + 2 epsilon). |
---|
3026 | double sMinPow = pow( s5min / s, xIntPF); |
---|
3027 | do { |
---|
3028 | xi1 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF ); |
---|
3029 | xi2 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF ); |
---|
3030 | s5 = xi1 * xi2 * s; |
---|
3031 | } while (s5 < s5min); |
---|
3032 | if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue; |
---|
3033 | |
---|
3034 | // Select t according to power and weigh by x_P^(2 alpha' |t|). |
---|
3035 | bool tryAgain = false; |
---|
3036 | for (int i = 0; i < 2; ++i) { |
---|
3037 | tVal[i] = - (1. / pow( tAux1[i] + rndmPtr->flat() |
---|
3038 | * (tAux2[i] - tAux1[i]), 1./3.) - 1.) / coefDL; |
---|
3039 | double wDL = pow2( (mp24DL - 2.8 * tVal[i]) / (mp24DL - tVal[i]) ) |
---|
3040 | / pow4( 1. - tVal[i] / 0.7); |
---|
3041 | double wMX = 1. / pow4( 1. - coefDL * tVal[i]); |
---|
3042 | if (wDL < rndmPtr->flat() * wMX) tryAgain = true; |
---|
3043 | double xi = (i == 0) ? xi1 : xi2; |
---|
3044 | if ( pow( xi, xtCorPF * abs(tVal[i]) ) < rndmPtr->flat() ) |
---|
3045 | tryAgain = true; |
---|
3046 | } |
---|
3047 | if (tryAgain) continue; |
---|
3048 | |
---|
3049 | // The MBR model (PomFlux == 5). |
---|
3050 | } else { |
---|
3051 | double dymin0 = 0.; |
---|
3052 | double dymax = log(s/m2minMBR); |
---|
3053 | double f1, f2, step2, dy, yc, ycmin, ycmax, dy1, dy2, |
---|
3054 | P, P1, P2, yRnd, yRnd1, yRnd2; |
---|
3055 | |
---|
3056 | // Von Neumann method to generate dy, uses dpepmax from SigmaTot. |
---|
3057 | do { |
---|
3058 | dy = dymin0 + (dymax - dymin0) * rndmPtr->flat(); |
---|
3059 | P = 0.; |
---|
3060 | step2 = (dy - dymin0) / NINTEG2; |
---|
3061 | for (int j = 0; j < NINTEG2 ; ++j) { |
---|
3062 | yc = -(dy - dymin0) / 2. + (j + 0.5) * step2; |
---|
3063 | dy1 = 0.5 * dy - yc; |
---|
3064 | dy2 = 0.5 * dy + yc; |
---|
3065 | f1 = exp(epsMBR * dy1) * ( (FFA1 / (FFB1 + 2. * alphMBR * dy1) ) |
---|
3066 | + (FFA2 / (FFB2 + 2. * alphMBR * dy1) ) ); |
---|
3067 | f2 = exp(epsMBR * dy2) * ( (FFA1 / (FFB1 + 2. * alphMBR * dy2) ) |
---|
3068 | + (FFA2 / (FFB2 + 2. * alphMBR * dy2) ) ); |
---|
3069 | f1 *= 0.5 * (1. + erf( (dy1 - 0.5 * dyminMBR) * dyminInvMBR )); |
---|
3070 | f2 *= 0.5 * (1. + erf( (dy2 - 0.5 * dyminMBR) * dyminInvMBR )); |
---|
3071 | P += f1 * f2 * step2; |
---|
3072 | } |
---|
3073 | if (P > dpepmax) { |
---|
3074 | ostringstream osWarn; |
---|
3075 | osWarn << "dpepmax = " << scientific << setprecision(3) |
---|
3076 | << dpepmax << " " << P << " " << dy << endl; |
---|
3077 | infoPtr->errorMsg("Warning in PhaseSpace2to2diffractive::" |
---|
3078 | "trialKin for central diffraction:", osWarn.str()); |
---|
3079 | } |
---|
3080 | yRnd = dpepmax * rndmPtr->flat(); |
---|
3081 | |
---|
3082 | // Generate dyc. |
---|
3083 | ycmax = (dy - dymin0) / 2.; |
---|
3084 | ycmin = -ycmax; |
---|
3085 | yc = ycmin + (ycmax - ycmin) * rndmPtr->flat(); |
---|
3086 | |
---|
3087 | // xi1, xi2 from dy and dy0. |
---|
3088 | dy1 = 0.5 * dy + yc; |
---|
3089 | dy2 = 0.5 * dy - yc; |
---|
3090 | P1 = 0.5 * (1. + erf( (dy1 - 0.5 * dyminMBR) * dyminInvMBR )); |
---|
3091 | P2 = 0.5 * (1. + erf( (dy2 - 0.5 * dyminMBR) * dyminInvMBR )); |
---|
3092 | yRnd1 = rndmPtr->flat(); |
---|
3093 | yRnd2 = rndmPtr->flat(); |
---|
3094 | } while( !(yRnd < P && yRnd1 < P1 && yRnd2 < P2) ); |
---|
3095 | xi1 = exp( -dy1 ); |
---|
3096 | xi2 = exp( -dy2 ); |
---|
3097 | |
---|
3098 | // Generate t1 at vertex1. First exponent, then FF*exp. |
---|
3099 | double tmin = -s1 * xi1 * xi1 / (1. - xi1); |
---|
3100 | do { |
---|
3101 | t1 = tmin + log(1. - rndmPtr->flat()); |
---|
3102 | double pFF = (4. * s1 - 2.8 * t1) / ( (4. * s1 - t1) |
---|
3103 | * pow2(1. - t1 / 0.71)); |
---|
3104 | P = pow2(pFF) * exp(2. * alphMBR * dy1 * t1); |
---|
3105 | yRnd = exp(t1) * rndmPtr->flat(); |
---|
3106 | } while (yRnd > P); |
---|
3107 | |
---|
3108 | // Generate t2 at vertex2. First exponent, then FF*exp. |
---|
3109 | tmin = -s2 * xi2 * xi2 / (1. - xi2); |
---|
3110 | do { |
---|
3111 | t2 = tmin + log(1. - rndmPtr->flat()); |
---|
3112 | double pFF = (4. * s2 - 2.8 * t2) / ((4. * s2 - t2) |
---|
3113 | * pow2(1. - t2 / 0.71)); |
---|
3114 | P = pow2(pFF) * exp(2. * alphMBR * dy2 * t2); |
---|
3115 | yRnd = exp(t2) * rndmPtr->flat(); |
---|
3116 | } while (yRnd > P); |
---|
3117 | } |
---|
3118 | |
---|
3119 | // Checks and kinematics construction four first options. |
---|
3120 | double pz3 = 0.; |
---|
3121 | double pz4 = 0.; |
---|
3122 | double pT3 = 0.; |
---|
3123 | double pT4 = 0.; |
---|
3124 | if (PomFlux < 5) { |
---|
3125 | |
---|
3126 | // Check whether m^2 (i.e. xi) and t choices are consistent. |
---|
3127 | bool tryAgain = false; |
---|
3128 | for (int i = 0; i < 2; ++i) { |
---|
3129 | double sx1 = (i == 0) ? s1 : s2; |
---|
3130 | double sx2 = (i == 0) ? s2 : s1; |
---|
3131 | double sx3 = sx1; |
---|
3132 | double sx4 = (i == 0) ? s2 + xi1 * s : s1 + xi2 * s; |
---|
3133 | if (sqrt(sx3) + sqrt(sx4) + DIFFMASSMARGIN > eCM) tryAgain = true; |
---|
3134 | double lambda34 = sqrtpos( pow2( s - sx3 - sx4) - 4. * sx3 * sx4 ); |
---|
3135 | double tempA = s - (sx1 + sx2 + sx3 + sx4) |
---|
3136 | + (sx1 - sx2) * (sx3 - sx4) / s; |
---|
3137 | double tempB = lambda12 * lambda34 / s; |
---|
3138 | double tempC = (sx3 - sx1) * (sx4 - sx2) + (sx1 + sx4 - sx2 - sx3) |
---|
3139 | * (sx1 * sx4 - sx2 * sx3) / s; |
---|
3140 | double tLowNow = -0.5 * (tempA + tempB); |
---|
3141 | double tUppNow = tempC / tLowNow; |
---|
3142 | if (tVal[i] < tLowNow || tVal[i] > tUppNow) tryAgain = true; |
---|
3143 | |
---|
3144 | // Careful reconstruction of scattering angle. |
---|
3145 | double cosTheta = min(1., max(-1., (tempA + 2. * tVal[i]) / tempB)); |
---|
3146 | double sinTheta = 2. * sqrtpos( -(tempC + tempA * tVal[i] |
---|
3147 | + tVal[i] * tVal[i]) ) / tempB; |
---|
3148 | theta = asin( min(1., sinTheta)); |
---|
3149 | if (cosTheta < 0.) theta = M_PI - theta; |
---|
3150 | double pAbs34 = 0.5 * lambda34 / eCM; |
---|
3151 | if (i == 0) { |
---|
3152 | pz3 = pAbs34 * cos(theta); |
---|
3153 | pT3 = pAbs34 * sin(theta); |
---|
3154 | } else { |
---|
3155 | pz4 = -pAbs34 * cos(theta); |
---|
3156 | pT4 = pAbs34 * sin(theta); |
---|
3157 | } |
---|
3158 | } |
---|
3159 | if (tryAgain) continue; |
---|
3160 | t1 = tVal[0]; |
---|
3161 | t2 = tVal[1]; |
---|
3162 | |
---|
3163 | // Kinematics construction in the MBR model. |
---|
3164 | } else { |
---|
3165 | pz3 = pAbs * (1. - xi1); |
---|
3166 | pz4 = -pAbs * (1. - xi2); |
---|
3167 | pT3 = sqrt( (1. - xi1) * abs(t1) - s1 * pow2(xi1) ); |
---|
3168 | pT4 = sqrt( (1. - xi2) * abs(t2) - s2 * pow2(xi2) ); |
---|
3169 | } |
---|
3170 | |
---|
3171 | // Common final steps of kinematics. |
---|
3172 | double phi3 = 2. * M_PI * rndmPtr->flat(); |
---|
3173 | double phi4 = 2. * M_PI * rndmPtr->flat(); |
---|
3174 | p3.p( pT3 * cos(phi3), pT3 * sin(phi3), pz3, |
---|
3175 | sqrt(pz3 * pz3 + pT3 * pT3 + s1) ); |
---|
3176 | p4.p( pT4 * cos(phi4), pT4 * sin(phi4), pz4, |
---|
3177 | sqrt(pz4 * pz4 + pT4 * pT4 + s2) ); |
---|
3178 | |
---|
3179 | // Central dissociated system, from Pomeron-Pomeron 4 vectors. |
---|
3180 | p5 = (p1 - p3) + (p2 - p4); |
---|
3181 | mHat = p5.mCalc(); |
---|
3182 | |
---|
3183 | // If acceptable diffractive mass then no more looping. |
---|
3184 | if (mHat > DIFFMASSMIN) break; |
---|
3185 | } |
---|
3186 | return true; |
---|
3187 | |
---|
3188 | } |
---|
3189 | |
---|
3190 | //-------------------------------------------------------------------------- |
---|
3191 | |
---|
3192 | // Construct the four-vector kinematics from the trial values. |
---|
3193 | |
---|
3194 | bool PhaseSpace2to3diffractive::finalKin() { |
---|
3195 | |
---|
3196 | // Particle four-momenta and masses. |
---|
3197 | pH[1] = p1; |
---|
3198 | pH[2] = p2; |
---|
3199 | pH[3] = p3; |
---|
3200 | pH[4] = p4; |
---|
3201 | pH[5] = p5; |
---|
3202 | mH[1] = mA; |
---|
3203 | mH[2] = mB; |
---|
3204 | mH[3] = mA; |
---|
3205 | mH[4] = mB; |
---|
3206 | mH[5] = mHat; |
---|
3207 | |
---|
3208 | // Set some further info for completeness (but Info can be for subprocess). |
---|
3209 | x1H = 1.; |
---|
3210 | x2H = 1.; |
---|
3211 | sH = s; |
---|
3212 | tH = (p1 - p3).m2Calc(); |
---|
3213 | uH = (p2 - p4).m2Calc(); |
---|
3214 | mHat = eCM; |
---|
3215 | p2Abs = pAbs * pAbs; |
---|
3216 | betaZ = 0.; |
---|
3217 | // Store average pT of three final particles for documentation. |
---|
3218 | pTH = (p3.pT() + p4.pT() + p5.pT()) / 3.; |
---|
3219 | |
---|
3220 | // Done. |
---|
3221 | return true; |
---|
3222 | |
---|
3223 | } |
---|
3224 | |
---|
3225 | //========================================================================== |
---|
3226 | |
---|
3227 | // PhaseSpace2to3tauycyl class. |
---|
3228 | // 2 -> 3 kinematics for normal subprocesses. |
---|
3229 | |
---|
3230 | //-------------------------------------------------------------------------- |
---|
3231 | |
---|
3232 | // Constants: could be changed here if desired, but normally should not. |
---|
3233 | // These are of technical nature, as described for each. |
---|
3234 | |
---|
3235 | // Number of Newton-Raphson iterations of kinematics when masses introduced. |
---|
3236 | const int PhaseSpace2to3tauycyl::NITERNR = 5; |
---|
3237 | |
---|
3238 | //-------------------------------------------------------------------------- |
---|
3239 | |
---|
3240 | // Set up for fixed or Breit-Wigner mass selection. |
---|
3241 | |
---|
3242 | bool PhaseSpace2to3tauycyl::setupMasses() { |
---|
3243 | |
---|
3244 | // Treat Z0 as such or as gamma*/Z0 |
---|
3245 | gmZmode = gmZmodeGlobal; |
---|
3246 | int gmZmodeProc = sigmaProcessPtr->gmZmode(); |
---|
3247 | if (gmZmodeProc >= 0) gmZmode = gmZmodeProc; |
---|
3248 | |
---|
3249 | // Set sHat limits - based on global limits only. |
---|
3250 | mHatMin = mHatGlobalMin; |
---|
3251 | sHatMin = mHatMin*mHatMin; |
---|
3252 | mHatMax = eCM; |
---|
3253 | if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax); |
---|
3254 | sHatMax = mHatMax*mHatMax; |
---|
3255 | |
---|
3256 | // Masses and widths of resonances. |
---|
3257 | setupMass1(3); |
---|
3258 | setupMass1(4); |
---|
3259 | setupMass1(5); |
---|
3260 | |
---|
3261 | // Reduced mass range - do not make it as fancy as in two-body case. |
---|
3262 | if (useBW[3]) mUpper[3] -= (mPeak[4] + mPeak[5]); |
---|
3263 | if (useBW[4]) mUpper[4] -= (mPeak[3] + mPeak[5]); |
---|
3264 | if (useBW[5]) mUpper[5] -= (mPeak[3] + mPeak[4]); |
---|
3265 | |
---|
3266 | // If closed phase space then unallowed process. |
---|
3267 | bool physical = true; |
---|
3268 | if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false; |
---|
3269 | if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false; |
---|
3270 | if (useBW[5] && mUpper[5] < mLower[5] + MASSMARGIN) physical = false; |
---|
3271 | if (!useBW[3] && !useBW[4] && !useBW[5] && mHatMax < mPeak[3] |
---|
3272 | + mPeak[4] + mPeak[5] + MASSMARGIN) physical = false; |
---|
3273 | if (!physical) return false; |
---|
3274 | |
---|
3275 | // No extra pT precautions in massless limit - assumed fixed by ME's. |
---|
3276 | pTHatMin = pTHatGlobalMin; |
---|
3277 | pT2HatMin = pTHatMin * pTHatMin; |
---|
3278 | pTHatMax = pTHatGlobalMax; |
---|
3279 | pT2HatMax = pTHatMax * pTHatMax; |
---|
3280 | |
---|
3281 | // Prepare to select m3 by BW + flat + 1/s_3. |
---|
3282 | if (useBW[3]) { |
---|
3283 | double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5]) |
---|
3284 | * mWidth[3] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5])); |
---|
3285 | double distToThreshB = (mHatMax - mPeak[3] - mMin[4] - mMin[5]) |
---|
3286 | / mWidth[3]; |
---|
3287 | double distToThresh = min( distToThreshA, distToThreshB); |
---|
3288 | setupMass2(3, distToThresh); |
---|
3289 | } |
---|
3290 | |
---|
3291 | // Prepare to select m4 by BW + flat + 1/s_3. |
---|
3292 | if (useBW[4]) { |
---|
3293 | double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5]) |
---|
3294 | * mWidth[4] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5])); |
---|
3295 | double distToThreshB = (mHatMax - mPeak[4] - mMin[3] - mMin[5]) |
---|
3296 | / mWidth[4]; |
---|
3297 | double distToThresh = min( distToThreshA, distToThreshB); |
---|
3298 | setupMass2(4, distToThresh); |
---|
3299 | } |
---|
3300 | |
---|
3301 | // Prepare to select m5 by BW + flat + 1/s_3. |
---|
3302 | if (useBW[5]) { |
---|
3303 | double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5]) |
---|
3304 | * mWidth[5] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5])); |
---|
3305 | double distToThreshB = (mHatMax - mPeak[5] - mMin[3] - mMin[4]) |
---|
3306 | / mWidth[5]; |
---|
3307 | double distToThresh = min( distToThreshA, distToThreshB); |
---|
3308 | setupMass2(5, distToThresh); |
---|
3309 | } |
---|
3310 | |
---|
3311 | // Initialization masses. For now give up when constrained phase space. |
---|
3312 | m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3]; |
---|
3313 | m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4]; |
---|
3314 | m5 = (useBW[5]) ? min(mPeak[5], mUpper[5]) : mPeak[5]; |
---|
3315 | if (m3 + m4 + m5 + MASSMARGIN > mHatMax) physical = false; |
---|
3316 | s3 = m3*m3; |
---|
3317 | s4 = m4*m4; |
---|
3318 | s5 = m5*m5; |
---|
3319 | |
---|
3320 | // Correct selected mass-spectrum to running-width Breit-Wigner. |
---|
3321 | // Extra safety margin for maximum search. |
---|
3322 | wtBW = 1.; |
---|
3323 | if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX; |
---|
3324 | if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX; |
---|
3325 | if (useBW[5]) wtBW *= weightMass(5) * EXTRABWWTMAX; |
---|
3326 | |
---|
3327 | // Done. |
---|
3328 | return physical; |
---|
3329 | |
---|
3330 | } |
---|
3331 | |
---|
3332 | //-------------------------------------------------------------------------- |
---|
3333 | |
---|
3334 | // Select Breit-Wigner-distributed or fixed masses. |
---|
3335 | |
---|
3336 | bool PhaseSpace2to3tauycyl::trialMasses() { |
---|
3337 | |
---|
3338 | // By default vanishing cross section. |
---|
3339 | sigmaNw = 0.; |
---|
3340 | wtBW = 1.; |
---|
3341 | |
---|
3342 | // Pick m3, m4 and m5 independently. |
---|
3343 | trialMass(3); |
---|
3344 | trialMass(4); |
---|
3345 | trialMass(5); |
---|
3346 | |
---|
3347 | // If outside phase space then reject event. |
---|
3348 | if (m3 + m4 + m5 + MASSMARGIN > mHatMax) return false; |
---|
3349 | |
---|
3350 | // Correct selected mass-spectrum to running-width Breit-Wigner. |
---|
3351 | if (useBW[3]) wtBW *= weightMass(3); |
---|
3352 | if (useBW[4]) wtBW *= weightMass(4); |
---|
3353 | if (useBW[5]) wtBW *= weightMass(5); |
---|
3354 | |
---|
3355 | // Done. |
---|
3356 | return true; |
---|
3357 | |
---|
3358 | } |
---|
3359 | |
---|
3360 | //-------------------------------------------------------------------------- |
---|
3361 | |
---|
3362 | // Construct the four-vector kinematics from the trial values. |
---|
3363 | |
---|
3364 | bool PhaseSpace2to3tauycyl::finalKin() { |
---|
3365 | |
---|
3366 | // Assign masses to particles assumed massless in matrix elements. |
---|
3367 | int id3 = sigmaProcessPtr->id(3); |
---|
3368 | int id4 = sigmaProcessPtr->id(4); |
---|
3369 | int id5 = sigmaProcessPtr->id(5); |
---|
3370 | if (idMass[3] == 0) { m3 = particleDataPtr->m0(id3); s3 = m3*m3; } |
---|
3371 | if (idMass[4] == 0) { m4 = particleDataPtr->m0(id4); s4 = m4*m4; } |
---|
3372 | if (idMass[5] == 0) { m5 = particleDataPtr->m0(id5); s5 = m5*m5; } |
---|
3373 | |
---|
3374 | // Check that phase space still open after new mass assignment. |
---|
3375 | if (m3 + m4 + m5 + MASSMARGIN > mHat) { |
---|
3376 | infoPtr->errorMsg("Warning in PhaseSpace2to3tauycyl::finalKin: " |
---|
3377 | "failed after mass assignment"); |
---|
3378 | return false; |
---|
3379 | } |
---|
3380 | |
---|
3381 | // Particle masses; incoming always on mass shell. |
---|
3382 | mH[1] = 0.; |
---|
3383 | mH[2] = 0.; |
---|
3384 | mH[3] = m3; |
---|
3385 | mH[4] = m4; |
---|
3386 | mH[5] = m5; |
---|
3387 | |
---|
3388 | // Incoming partons along beam axes. |
---|
3389 | pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H); |
---|
3390 | pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H); |
---|
3391 | |
---|
3392 | // Begin three-momentum rescaling to compensate for masses. |
---|
3393 | if (idMass[3] == 0 || idMass[4] == 0 || idMass[5] == 0) { |
---|
3394 | double p3S = p3cm.pAbs2(); |
---|
3395 | double p4S = p4cm.pAbs2(); |
---|
3396 | double p5S = p5cm.pAbs2(); |
---|
3397 | double fac = 1.; |
---|
3398 | double e3, e4, e5, value, deriv; |
---|
3399 | |
---|
3400 | // Iterate rescaling solution five times, using Newton-Raphson. |
---|
3401 | for (int i = 0; i < NITERNR; ++i) { |
---|
3402 | e3 = sqrt(s3 + fac * p3S); |
---|
3403 | e4 = sqrt(s4 + fac * p4S); |
---|
3404 | e5 = sqrt(s5 + fac * p5S); |
---|
3405 | value = e3 + e4 + e5 - mHat; |
---|
3406 | deriv = 0.5 * (p3S / e3 + p4S / e4 + p5S / e5); |
---|
3407 | fac -= value / deriv; |
---|
3408 | } |
---|
3409 | |
---|
3410 | // Rescale momenta appropriately. |
---|
3411 | double facRoot = sqrt(fac); |
---|
3412 | p3cm.rescale3( facRoot ); |
---|
3413 | p4cm.rescale3( facRoot ); |
---|
3414 | p5cm.rescale3( facRoot ); |
---|
3415 | p3cm.e( sqrt(s3 + fac * p3S) ); |
---|
3416 | p4cm.e( sqrt(s4 + fac * p4S) ); |
---|
3417 | p5cm.e( sqrt(s5 + fac * p5S) ); |
---|
3418 | } |
---|
3419 | |
---|
3420 | // Outgoing partons initially in collision CM frame along beam axes. |
---|
3421 | pH[3] = p3cm; |
---|
3422 | pH[4] = p4cm; |
---|
3423 | pH[5] = p5cm; |
---|
3424 | |
---|
3425 | // Then boost them to overall CM frame |
---|
3426 | betaZ = (x1H - x2H)/(x1H + x2H); |
---|
3427 | pH[3].rot( theta, phi); |
---|
3428 | pH[4].rot( theta, phi); |
---|
3429 | pH[3].bst( 0., 0., betaZ); |
---|
3430 | pH[4].bst( 0., 0., betaZ); |
---|
3431 | pH[5].bst( 0., 0., betaZ); |
---|
3432 | |
---|
3433 | // Store average pT of three final particles for documentation. |
---|
3434 | pTH = (p3cm.pT() + p4cm.pT() + p5cm.pT()) / 3.; |
---|
3435 | |
---|
3436 | // Done. |
---|
3437 | return true; |
---|
3438 | } |
---|
3439 | |
---|
3440 | //========================================================================== |
---|
3441 | |
---|
3442 | // The PhaseSpace2to3yyycyl class. |
---|
3443 | // Phase space for 2 -> 3 QCD processes, 1 + 2 -> 3 + 4 + 5 set up in |
---|
3444 | // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut. |
---|
3445 | // Note: here cout is used for output, not os. Change?? |
---|
3446 | |
---|
3447 | //-------------------------------------------------------------------------- |
---|
3448 | |
---|
3449 | // Sample the phase space of the process. |
---|
3450 | |
---|
3451 | bool PhaseSpace2to3yyycyl::setupSampling() { |
---|
3452 | |
---|
3453 | // Phase space cuts specifically for 2 -> 3 QCD processes. |
---|
3454 | pTHat3Min = settingsPtr->parm("PhaseSpace:pTHat3Min"); |
---|
3455 | pTHat3Max = settingsPtr->parm("PhaseSpace:pTHat3Max"); |
---|
3456 | pTHat5Min = settingsPtr->parm("PhaseSpace:pTHat5Min"); |
---|
3457 | pTHat5Max = settingsPtr->parm("PhaseSpace:pTHat5Max"); |
---|
3458 | RsepMin = settingsPtr->parm("PhaseSpace:RsepMin"); |
---|
3459 | R2sepMin = pow2(RsepMin); |
---|
3460 | |
---|
3461 | // If both beams are baryons then softer PDF's than for mesons/Pomerons. |
---|
3462 | hasBaryonBeams = ( beamAPtr->isBaryon() && beamBPtr->isBaryon() ); |
---|
3463 | |
---|
3464 | // Work with massless partons. |
---|
3465 | for (int i = 0; i < 6; ++i) mH[i] = 0.; |
---|
3466 | |
---|
3467 | // Constrain to possible cuts at current CM energy and check consistency. |
---|
3468 | pT3Min = pTHat3Min; |
---|
3469 | pT3Max = pTHat3Max; |
---|
3470 | if (pT3Max < pT3Min) pT3Max = 0.5 * eCM; |
---|
3471 | pT5Min = pTHat5Min; |
---|
3472 | pT5Max = pTHat5Max; |
---|
3473 | if (pT5Max < pT5Min) pT5Max = 0.5 * eCM; |
---|
3474 | if (pT5Max > pT3Max || pT5Min > pT3Min || pT3Min + 2. * pT5Min > eCM) { |
---|
3475 | infoPtr->errorMsg("Error in PhaseSpace2to3yyycyl::setupSampling: " |
---|
3476 | "inconsistent pT limits in 3-body phase space"); |
---|
3477 | return false; |
---|
3478 | } |
---|
3479 | |
---|
3480 | // Loop over some configurations where cross section could be maximal. |
---|
3481 | // In all cases all sum p_z = 0, for maximal PDF weights. |
---|
3482 | // Also pT3 and R45 are minimal, while pT5 may vary. |
---|
3483 | sigmaMx = 0.; |
---|
3484 | double pT5EffMax = min( pT5Max, 0.5 * pT3Min / cos(0.5 * RsepMin) ); |
---|
3485 | double pT3EffMin = max( pT3Min, 2.0 * pT5Min * cos(0.5 * RsepMin) ) ; |
---|
3486 | double sinhR = sinh(0.5 * RsepMin); |
---|
3487 | double coshR = cosh(0.5 * RsepMin); |
---|
3488 | for (int iStep = 0; iStep < 120; ++iStep) { |
---|
3489 | |
---|
3490 | // First kind: |phi4 - phi5| = R, all p_z = 0, i.e. separation in phi. |
---|
3491 | if (iStep < 10) { |
---|
3492 | pT3 = pT3EffMin; |
---|
3493 | pT5 = pT5Min * pow( pT5EffMax / pT5Min, iStep / 9.); |
---|
3494 | double pTRat = pT5 / pT3; |
---|
3495 | double sin2Rsep = pow2( sin(RsepMin) ); |
---|
3496 | double cosPhi35 = - cos(RsepMin) * sqrtpos(1. - sin2Rsep |
---|
3497 | * pow2(pTRat)) - sin2Rsep * pTRat; |
---|
3498 | cosPhi35 = max( cosPhi35, cos(M_PI - 0.5 * RsepMin) ); |
---|
3499 | double sinPhi35 = sqrt(1. - pow2(cosPhi35)); |
---|
3500 | pT4 = sqrt( pow2(pT3) + pow2(pT5) + 2. * pT3 * pT5 * cosPhi35); |
---|
3501 | p3cm = pT3 * Vec4( 1., 0., 0., 1.); |
---|
3502 | p4cm = Vec4(-pT3 - pT5 * cosPhi35, -pT5 * sinPhi35, 0., pT4); |
---|
3503 | p5cm = pT5 * Vec4( cosPhi35, sinPhi35, 0., 1.); |
---|
3504 | |
---|
3505 | // Second kind: |y4 - y5| = R, phi4 = phi5, i.e. separation in y. |
---|
3506 | } else { |
---|
3507 | pT5 = pT5Min * pow( pT5Max / pT5Min, iStep%10 / 9. ); |
---|
3508 | pT3 = max( pT3Min, 2. * pT5); |
---|
3509 | pT4 = pT3 - pT5; |
---|
3510 | p4cm = pT4 * Vec4( -1., 0., sinhR, coshR ); |
---|
3511 | p5cm = pT5 * Vec4( -1., 0., -sinhR, coshR ); |
---|
3512 | y3 = -1.2 + 0.2 * (iStep/10); |
---|
3513 | p3cm = pT3 * Vec4( 1., 0., sinh(y3), cosh(y3)); |
---|
3514 | betaZ = (p3cm.pz() + p4cm.pz() + p5cm.pz()) |
---|
3515 | / (p3cm.e() + p4cm.e() + p5cm.e()); |
---|
3516 | p3cm.bst( 0., 0., -betaZ); |
---|
3517 | p4cm.bst( 0., 0., -betaZ); |
---|
3518 | p5cm.bst( 0., 0., -betaZ); |
---|
3519 | } |
---|
3520 | |
---|
3521 | // Find cross section in chosen phase space point. |
---|
3522 | pInSum = p3cm + p4cm + p5cm; |
---|
3523 | x1H = pInSum.e() / eCM; |
---|
3524 | x2H = x1H; |
---|
3525 | sH = pInSum.m2Calc(); |
---|
3526 | sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm, |
---|
3527 | 0., 0., 0., 1., 1., 1.); |
---|
3528 | sigmaNw = sigmaProcessPtr->sigmaPDF(); |
---|
3529 | |
---|
3530 | // Multiply by Jacobian. |
---|
3531 | double flux = 1. /(8. * pow2(sH) * pow5(2. * M_PI)); |
---|
3532 | double pTRng = pow2(M_PI) |
---|
3533 | * pow4(pT3) * (1./pow2(pT3Min) - 1./pow2(pT3Max)) |
---|
3534 | * pow2(pT5) * 2.* log(pT5Max/pT5Min); |
---|
3535 | double yRng = 8. * log(eCM / pT3) * log(eCM / pT4) * log(eCM / pT5); |
---|
3536 | sigmaNw *= SAFETYMARGIN * flux * pTRng * yRng; |
---|
3537 | |
---|
3538 | // Update to largest maximum. |
---|
3539 | if (showSearch && sigmaNw > sigmaMx) cout << "\n New sigmamax is " |
---|
3540 | << scientific << setprecision(3) << sigmaNw << " for x1 = " << x1H |
---|
3541 | << " x2 = " << x2H << " sH = " << sH << endl << " p3 = " << p3cm |
---|
3542 | << " p4 = " << p4cm << " p5 = " << p5cm; |
---|
3543 | if (sigmaNw > sigmaMx) sigmaMx = sigmaNw; |
---|
3544 | } |
---|
3545 | sigmaPos = sigmaMx; |
---|
3546 | |
---|
3547 | // Done. |
---|
3548 | return true; |
---|
3549 | } |
---|
3550 | |
---|
3551 | //-------------------------------------------------------------------------- |
---|
3552 | |
---|
3553 | // Sample the phase space of the process. |
---|
3554 | |
---|
3555 | bool PhaseSpace2to3yyycyl::trialKin(bool inEvent, bool) { |
---|
3556 | |
---|
3557 | // Allow for possibility that energy varies from event to event. |
---|
3558 | if (doEnergySpread) { |
---|
3559 | eCM = infoPtr->eCM(); |
---|
3560 | s = eCM * eCM; |
---|
3561 | } |
---|
3562 | sigmaNw = 0.; |
---|
3563 | |
---|
3564 | // Constrain to possible cuts at current CM energy and check consistency. |
---|
3565 | pT3Min = pTHat3Min; |
---|
3566 | pT3Max = pTHat3Max; |
---|
3567 | if (pT3Max < pT3Min) pT3Max = 0.5 * eCM; |
---|
3568 | pT5Min = pTHat5Min; |
---|
3569 | pT5Max = pTHat5Max; |
---|
3570 | if (pT5Max < pT5Min) pT5Max = 0.5 * eCM; |
---|
3571 | if (pT5Max > pT3Max || pT5Min > pT3Min || pT3Min + 2. * pT5Min > eCM) { |
---|
3572 | infoPtr->errorMsg("Error in PhaseSpace2to3yyycyl::trialKin: " |
---|
3573 | "inconsistent pT limits in 3-body phase space"); |
---|
3574 | return false; |
---|
3575 | } |
---|
3576 | |
---|
3577 | // Pick pT3 according to d^2(pT3)/pT3^4 and pT5 to d^2(pT5)/pT5^2. |
---|
3578 | pT3 = pT3Min * pT3Max / sqrt( pow2(pT3Min) + |
---|
3579 | rndmPtr->flat() * (pow2(pT3Max) - pow2(pT3Min)) ); |
---|
3580 | pT5Max = min(pT5Max, pT3); |
---|
3581 | if (pT5Max < pT5Min) return false; |
---|
3582 | pT5 = pT5Min * pow( pT5Max / pT5Min, rndmPtr->flat() ); |
---|
3583 | |
---|
3584 | // Pick azimuthal angles flat and reconstruct pT4, between pT3 and pT5. |
---|
3585 | phi3 = 2. * M_PI * rndmPtr->flat(); |
---|
3586 | phi5 = 2. * M_PI * rndmPtr->flat(); |
---|
3587 | pT4 = sqrt( pow2(pT3) + pow2(pT5) + 2. * pT3 * pT5 * cos(phi3 - phi5) ); |
---|
3588 | if (pT4 > pT3 || pT4 < pT5) return false; |
---|
3589 | phi4 = atan2( -(pT3 * sin(phi3) + pT5 * sin(phi5)), |
---|
3590 | -(pT3 * cos(phi3) + pT5 * cos(phi5)) ); |
---|
3591 | |
---|
3592 | // Pick rapidities flat in allowed ranges. |
---|
3593 | y3Max = log(eCM / pT3); |
---|
3594 | y4Max = log(eCM / pT4); |
---|
3595 | y5Max = log(eCM / pT5); |
---|
3596 | y3 = y3Max * (2. * rndmPtr->flat() - 1.); |
---|
3597 | y4 = y4Max * (2. * rndmPtr->flat() - 1.); |
---|
3598 | y5 = y5Max * (2. * rndmPtr->flat() - 1.); |
---|
3599 | |
---|
3600 | // Reject some events at large rapidities to improve efficiency. |
---|
3601 | // (Works for baryons, not pions or Pomerons if they have hard PDF's.) |
---|
3602 | double WTy = (hasBaryonBeams) ? (1. - pow2(y3/y3Max)) |
---|
3603 | * (1. - pow2(y4/y4Max)) * (1. - pow2(y5/y5Max)) : 1.; |
---|
3604 | if (WTy < rndmPtr->flat()) return false; |
---|
3605 | |
---|
3606 | // Check that any pair separated more then RsepMin in (y, phi) space. |
---|
3607 | dphi = abs(phi3 - phi4); |
---|
3608 | if (dphi > M_PI) dphi = 2. * M_PI - dphi; |
---|
3609 | if (pow2(y3 - y4) + pow2(dphi) < R2sepMin) return false; |
---|
3610 | dphi = abs(phi3 - phi5); |
---|
3611 | if (dphi > M_PI) dphi = 2. * M_PI - dphi; |
---|
3612 | if (pow2(y3 - y5) + pow2(dphi) < R2sepMin) return false; |
---|
3613 | dphi = abs(phi4 - phi5); |
---|
3614 | if (dphi > M_PI) dphi = 2. * M_PI - dphi; |
---|
3615 | if (pow2(y4 - y5) + pow2(dphi) < R2sepMin) return false; |
---|
3616 | |
---|
3617 | // Reconstruct all four-vectors. |
---|
3618 | pH[3] = pT3 * Vec4( cos(phi3), sin(phi3), sinh(y3), cosh(y3) ); |
---|
3619 | pH[4] = pT4 * Vec4( cos(phi4), sin(phi4), sinh(y4), cosh(y4) ); |
---|
3620 | pH[5] = pT5 * Vec4( cos(phi5), sin(phi5), sinh(y5), cosh(y5) ); |
---|
3621 | pInSum = pH[3] + pH[4] + pH[5]; |
---|
3622 | |
---|
3623 | // Check that x values physical and sHat in allowed range. |
---|
3624 | x1H = (pInSum.e() + pInSum.pz()) / eCM; |
---|
3625 | x2H = (pInSum.e() - pInSum.pz()) / eCM; |
---|
3626 | if (x1H >= 1. || x2H >= 1.) return false; |
---|
3627 | sH = pInSum.m2Calc(); |
---|
3628 | if ( sH < pow2(mHatGlobalMin) || |
---|
3629 | (mHatGlobalMax > mHatGlobalMin && sH > pow2(mHatGlobalMax)) ) |
---|
3630 | return false; |
---|
3631 | |
---|
3632 | // Boost four-vectors to rest frame of collision. |
---|
3633 | betaZ = (x1H - x2H)/(x1H + x2H); |
---|
3634 | p3cm = pH[3]; p3cm.bst( 0., 0., -betaZ); |
---|
3635 | p4cm = pH[4]; p4cm.bst( 0., 0., -betaZ); |
---|
3636 | p5cm = pH[5]; p5cm.bst( 0., 0., -betaZ); |
---|
3637 | |
---|
3638 | // Find cross section in chosen phase space point. |
---|
3639 | sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm, |
---|
3640 | 0., 0., 0., 1., 1., 1.); |
---|
3641 | sigmaNw = sigmaProcessPtr->sigmaPDF(); |
---|
3642 | |
---|
3643 | // Multiply by Jacobian. Correct for rejection of large rapidities. |
---|
3644 | double flux = 1. /(8. * pow2(sH) * pow5(2. * M_PI)); |
---|
3645 | double yRng = 8. * y3Max * y4Max * y5Max; |
---|
3646 | double pTRng = pow2(M_PI) |
---|
3647 | * pow4(pT3) * (1./pow2(pT3Min) - 1./pow2(pT3Max)) |
---|
3648 | * pow2(pT5) * 2.* log(pT5Max/pT5Min); |
---|
3649 | sigmaNw *= flux * yRng * pTRng / WTy; |
---|
3650 | |
---|
3651 | // Allow possibility for user to modify cross section. |
---|
3652 | if (canModifySigma) sigmaNw |
---|
3653 | *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, inEvent); |
---|
3654 | if (canBiasSelection) sigmaNw |
---|
3655 | *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, inEvent); |
---|
3656 | if (canBias2Sel) sigmaNw *= pow( pTH / bias2SelRef, bias2SelPow); |
---|
3657 | |
---|
3658 | // Check if maximum violated. |
---|
3659 | newSigmaMx = false; |
---|
3660 | if (sigmaNw > sigmaMx) { |
---|
3661 | infoPtr->errorMsg("Warning in PhaseSpace2to3yyycyl::trialKin: " |
---|
3662 | "maximum for cross section violated"); |
---|
3663 | |
---|
3664 | // Violation strategy 1: increase maximum (always during initialization). |
---|
3665 | if (increaseMaximum || !inEvent) { |
---|
3666 | double violFact = SAFETYMARGIN * sigmaNw / sigmaMx; |
---|
3667 | sigmaMx = SAFETYMARGIN * sigmaNw; |
---|
3668 | newSigmaMx = true; |
---|
3669 | if (showViolation) { |
---|
3670 | if (violFact < 9.99) cout << fixed; |
---|
3671 | else cout << scientific; |
---|
3672 | cout << " PYTHIA Maximum for " << sigmaProcessPtr->name() |
---|
3673 | << " increased by factor " << setprecision(3) << violFact |
---|
3674 | << " to " << scientific << sigmaMx << endl; |
---|
3675 | } |
---|
3676 | |
---|
3677 | // Violation strategy 2: weight event (done in ProcessContainer). |
---|
3678 | } else if (showViolation && sigmaNw > sigmaPos) { |
---|
3679 | double violFact = sigmaNw / sigmaMx; |
---|
3680 | if (violFact < 9.99) cout << fixed; |
---|
3681 | else cout << scientific; |
---|
3682 | cout << " PYTHIA Maximum for " << sigmaProcessPtr->name() |
---|
3683 | << " exceeded by factor " << setprecision(3) << violFact << endl; |
---|
3684 | sigmaPos = sigmaNw; |
---|
3685 | } |
---|
3686 | } |
---|
3687 | |
---|
3688 | // Check if negative cross section. |
---|
3689 | if (sigmaNw < sigmaNeg) { |
---|
3690 | infoPtr->errorMsg("Warning in PhaseSpace2to3yyycyl::trialKin:" |
---|
3691 | " negative cross section set 0", "for " + sigmaProcessPtr->name() ); |
---|
3692 | sigmaNeg = sigmaNw; |
---|
3693 | |
---|
3694 | // Optional printout of (all) violations. |
---|
3695 | if (showViolation) cout << " PYTHIA Negative minimum for " |
---|
3696 | << sigmaProcessPtr->name() << " changed to " << scientific |
---|
3697 | << setprecision(3) << sigmaNeg << endl; |
---|
3698 | } |
---|
3699 | if (sigmaNw < 0.) sigmaNw = 0.; |
---|
3700 | |
---|
3701 | |
---|
3702 | // Done. |
---|
3703 | return true; |
---|
3704 | } |
---|
3705 | |
---|
3706 | //-------------------------------------------------------------------------- |
---|
3707 | |
---|
3708 | // Construct the final kinematics of the process: not much left |
---|
3709 | |
---|
3710 | bool PhaseSpace2to3yyycyl::finalKin() { |
---|
3711 | |
---|
3712 | // Work with massless partons. |
---|
3713 | for (int i = 0; i < 6; ++i) mH[i] = 0.; |
---|
3714 | |
---|
3715 | // Ibncoming partons to collision. |
---|
3716 | pH[1] = 0.5 * (pInSum.e() + pInSum.pz()) * Vec4( 0., 0., 1., 1.); |
---|
3717 | pH[2] = 0.5 * (pInSum.e() - pInSum.pz()) * Vec4( 0., 0., -1., 1.); |
---|
3718 | |
---|
3719 | // Some quantities meaningless for 2 -> 3. pT devined as average value. |
---|
3720 | tH = 0.; |
---|
3721 | uH = 0.; |
---|
3722 | pTH = (pH[3].pT() + pH[4].pT() + pH[5].pT()) / 3.; |
---|
3723 | theta = 0.; |
---|
3724 | phi = 0.; |
---|
3725 | |
---|
3726 | return true; |
---|
3727 | } |
---|
3728 | |
---|
3729 | |
---|
3730 | //========================================================================== |
---|
3731 | |
---|
3732 | // The PhaseSpaceLHA class. |
---|
3733 | // A derived class for Les Houches events. |
---|
3734 | // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks. |
---|
3735 | |
---|
3736 | //-------------------------------------------------------------------------- |
---|
3737 | |
---|
3738 | // Constants: could be changed here if desired, but normally should not. |
---|
3739 | // These are of technical nature, as described for each. |
---|
3740 | |
---|
3741 | // LHA convention with cross section in pb forces conversion to mb. |
---|
3742 | const double PhaseSpaceLHA::CONVERTPB2MB = 1e-9; |
---|
3743 | |
---|
3744 | //-------------------------------------------------------------------------- |
---|
3745 | |
---|
3746 | // Find maximal cross section for comparison with internal processes. |
---|
3747 | |
---|
3748 | bool PhaseSpaceLHA::setupSampling() { |
---|
3749 | |
---|
3750 | // Find which strategy Les Houches events are produced with. |
---|
3751 | strategy = lhaUpPtr->strategy(); |
---|
3752 | stratAbs = abs(strategy); |
---|
3753 | if (strategy == 0 || stratAbs > 4) { |
---|
3754 | ostringstream stratCode; |
---|
3755 | stratCode << strategy; |
---|
3756 | infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: unknown " |
---|
3757 | "Les Houches Accord weighting stategy", stratCode.str()); |
---|
3758 | return false; |
---|
3759 | } |
---|
3760 | |
---|
3761 | // Number of contributing processes. |
---|
3762 | nProc = lhaUpPtr->sizeProc(); |
---|
3763 | |
---|
3764 | // Loop over all processes. Read out maximum and cross section. |
---|
3765 | xMaxAbsSum = 0.; |
---|
3766 | xSecSgnSum = 0.; |
---|
3767 | int idPr; |
---|
3768 | double xMax, xSec, xMaxAbs; |
---|
3769 | for (int iProc = 0 ; iProc < nProc; ++iProc) { |
---|
3770 | idPr = lhaUpPtr->idProcess(iProc); |
---|
3771 | xMax = lhaUpPtr->xMax(iProc); |
---|
3772 | xSec = lhaUpPtr->xSec(iProc); |
---|
3773 | |
---|
3774 | // Check for inconsistencies between strategy and stored values. |
---|
3775 | if ( (strategy == 1 || strategy == 2) && xMax < 0.) { |
---|
3776 | infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: " |
---|
3777 | "negative maximum not allowed"); |
---|
3778 | return false; |
---|
3779 | } |
---|
3780 | if ( ( strategy == 2 || strategy == 3) && xSec < 0.) { |
---|
3781 | infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: " |
---|
3782 | "negative cross section not allowed"); |
---|
3783 | return false; |
---|
3784 | } |
---|
3785 | |
---|
3786 | // Store maximal cross sections for later choice. |
---|
3787 | if (stratAbs == 1) xMaxAbs = abs(xMax); |
---|
3788 | else if (stratAbs < 4) xMaxAbs = abs(xSec); |
---|
3789 | else xMaxAbs = 1.; |
---|
3790 | idProc.push_back( idPr ); |
---|
3791 | xMaxAbsProc.push_back( xMaxAbs ); |
---|
3792 | |
---|
3793 | // Find sum and convert to mb. |
---|
3794 | xMaxAbsSum += xMaxAbs; |
---|
3795 | xSecSgnSum += xSec; |
---|
3796 | } |
---|
3797 | sigmaMx = xMaxAbsSum * CONVERTPB2MB; |
---|
3798 | sigmaSgn = xSecSgnSum * CONVERTPB2MB; |
---|
3799 | |
---|
3800 | // Done. |
---|
3801 | return true; |
---|
3802 | |
---|
3803 | } |
---|
3804 | |
---|
3805 | //-------------------------------------------------------------------------- |
---|
3806 | |
---|
3807 | // Construct the next process, by interface to Les Houches class. |
---|
3808 | |
---|
3809 | bool PhaseSpaceLHA::trialKin( bool, bool repeatSame ) { |
---|
3810 | |
---|
3811 | // Must select process type in some cases. |
---|
3812 | int idProcNow = 0; |
---|
3813 | if (repeatSame) idProcNow = idProcSave; |
---|
3814 | else if (stratAbs <= 2) { |
---|
3815 | double xMaxAbsRndm = xMaxAbsSum * rndmPtr->flat(); |
---|
3816 | int iProc = -1; |
---|
3817 | do xMaxAbsRndm -= xMaxAbsProc[++iProc]; |
---|
3818 | while (xMaxAbsRndm > 0. && iProc < nProc - 1); |
---|
3819 | idProcNow = idProc[iProc]; |
---|
3820 | } |
---|
3821 | |
---|
3822 | // Generate Les Houches event. Return if fail (= end of file). |
---|
3823 | bool physical = lhaUpPtr->setEvent(idProcNow); |
---|
3824 | if (!physical) return false; |
---|
3825 | |
---|
3826 | // Find which process was generated. |
---|
3827 | int idPr = lhaUpPtr->idProcess(); |
---|
3828 | int iProc = 0; |
---|
3829 | for (int iP = 0; iP < int(idProc.size()); ++iP) |
---|
3830 | if (idProc[iP] == idPr) iProc = iP; |
---|
3831 | idProcSave = idPr; |
---|
3832 | |
---|
3833 | // Extract cross section and rescale according to strategy. |
---|
3834 | double wtPr = lhaUpPtr->weight(); |
---|
3835 | if (stratAbs == 1) sigmaNw = wtPr * CONVERTPB2MB |
---|
3836 | * xMaxAbsSum / xMaxAbsProc[iProc]; |
---|
3837 | else if (stratAbs == 2) sigmaNw = (wtPr / abs(lhaUpPtr->xMax(iProc))) |
---|
3838 | * sigmaMx; |
---|
3839 | else if (strategy == 3) sigmaNw = sigmaMx; |
---|
3840 | else if (strategy == -3 && wtPr > 0.) sigmaNw = sigmaMx; |
---|
3841 | else if (strategy == -3) sigmaNw = -sigmaMx; |
---|
3842 | else if (stratAbs == 4) sigmaNw = wtPr * CONVERTPB2MB; |
---|
3843 | |
---|
3844 | // Set x scales. |
---|
3845 | x1H = lhaUpPtr->x1(); |
---|
3846 | x2H = lhaUpPtr->x2(); |
---|
3847 | |
---|
3848 | // Done. |
---|
3849 | return true; |
---|
3850 | |
---|
3851 | } |
---|
3852 | |
---|
3853 | //========================================================================== |
---|
3854 | |
---|
3855 | } // end namespace Pythia8 |
---|