using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / PhaseSpace.cxx
CommitLineData
5ad4eb21 1// PhaseSpace.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 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
11namespace 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.
24const int PhaseSpace::NMAXTRY = 2;
25
26// Number of three-body trials in phase space optimization.
27const int PhaseSpace::NTRY3BODY = 20;
28
29// Maximum cross section increase, just in case true maximum not found.
30const double PhaseSpace::SAFETYMARGIN = 1.05;
31
32// Small number to avoid division by zero.
33const double PhaseSpace::TINY = 1e-20;
34
35// Fraction of total weight that is shared evenly between all shapes.
36const double PhaseSpace::EVENFRAC = 0.4;
37
38// Two cross sections with a small relative error are assumed same.
39const double PhaseSpace::SAMESIGMA = 1e-6;
40
41// Do not include resonances peaked too far outside allowed mass region.
42const double PhaseSpace::WIDTHMARGIN = 20.;
43
44// Special optimization treatment when two resonances at almost same mass.
45const double PhaseSpace::SAMEMASS = 0.01;
46
47// Minimum phase space left when kinematics constraints are combined.
48const double PhaseSpace::MASSMARGIN = 0.01;
49
50// When using Breit-Wigners in 2 -> 2 raise maximum weight estimate.
51const double PhaseSpace::EXTRABWWTMAX = 1.25;
52
53// Size of Breit-Wigner threshold region, for mass selection biasing.
54const double PhaseSpace::THRESHOLDSIZE = 3.;
55
56// Step size in optimal-mass search, for mass selection biasing.
57const double PhaseSpace::THRESHOLDSTEP = 0.2;
58
59// Minimal rapidity range for allowed open range (in 2 -> 3).
60const 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.
64const double PhaseSpace::LEPTONXMIN = 1e-10;
65const double PhaseSpace::LEPTONXMAX = 1. - 1e-10;
66const double PhaseSpace::LEPTONXLOGMIN = log(1e-10);
67const double PhaseSpace::LEPTONXLOGMAX = log(1. - 1e-10);
68const double PhaseSpace::LEPTONTAUMIN = 2e-10;
69
70// Safety to avoid division with unreasonably small value for z selection.
71const double PhaseSpace::SHATMINZ = 1.;
72
73// Regularization for small pT2min in z = cos(theta) selection.
74const double PhaseSpace::PT2RATMINZ = 0.0001;
75
76// These numbers are hardwired empirical parameters,
77// intended to speed up the M-generator.
78const double PhaseSpace::WTCORRECTION[11] = { 1., 1., 1.,
79 2., 5., 15., 60., 250., 1250., 7000., 50000. };
80
81//*********
82
83// Perform simple initialization and store pointers.
84
85void PhaseSpace::init(SigmaProcess* sigmaProcessPtrIn, Info* infoPtrIn,
86 BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
87 SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn) {
88
89 // Store input pointers for future use.
90 sigmaProcessPtr = sigmaProcessPtrIn;
91 infoPtr = infoPtrIn;
92 beamAPtr = beamAPtrIn;
93 beamBPtr = beamBPtrIn;
94 sigmaTotPtr = sigmaTotPtrIn;
95 userHooksPtr = userHooksPtrIn;
96
97 // Some commonly used beam information.
98 idA = beamAPtr->id();
99 idB = beamBPtr->id();
100 mA = beamAPtr->m();
101 mB = beamBPtr->m();
102 eCM = infoPtr->eCM();
103 s = eCM * eCM;
104
105 // Flag if lepton beams, and if non-resolved ones.
106 hasLeptonBeams = ( beamAPtr->isLepton() || beamBPtr->isLepton() );
107 hasPointLeptons = ( hasLeptonBeams
108 && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) );
109
110 // Standard phase space cuts.
111 mHatGlobalMin = Settings::parm("PhaseSpace:mHatMin");
112 mHatGlobalMax = Settings::parm("PhaseSpace:mHatMax");
113 pTHatGlobalMin = Settings::parm("PhaseSpace:pTHatMin");
114 pTHatGlobalMax = Settings::parm("PhaseSpace:pTHatMax");
115 pTHatMinDiverge = Settings::parm("PhaseSpace:pTHatMinDiverge");
116
117 // When to use Breit-Wigners.
118 useBreitWigners = Settings::flag("PhaseSpace:useBreitWigners");
119 minWidthBreitWigners = Settings::parm("PhaseSpace:minWidthBreitWigners");
120
121 // Whether generation is with variable energy.
122 doEnergySpread = Settings::flag("Beams:allowMomentumSpread");
123
124 // Print flag for maximization information.
125 showSearch = Settings::flag("PhaseSpace:showSearch");
126 showViolation = Settings::flag("PhaseSpace:showViolation");
127
128 // Know whether a Z0 is pure Z0 or admixed with gamma*.
129 gmZmodeGlobal = Settings::mode("WeakZ0:gmZmode");
130
131 // Default event-specific kinematics properties.
132 x1H = 1.;
133 x2H = 1.;
134 m3 = 0.;
135 m4 = 0.;
136 m5 = 0.;
137 s3 = m3 * m3;
138 s4 = m4 * m4;
139 s5 = m5 * m5;
140 mHat = eCM;
141 sH = s;
142 tH = 0.;
143 uH = 0.;
144 pTH = 0.;
145 theta = 0.;
146 phi = 0.;
147 runBW3H = 1.;
148 runBW4H = 1.;
149 runBW5H = 1.;
150
151 // Default cross section information.
152 sigmaNw = 0.;
153 sigmaMx = 0.;
154 sigmaNeg = 0.;
155 newSigmaMx = false;
156
157 // Flag if user should be allow to reweight cross section.
158 canModifySigma = (userHooksPtr > 0)
159 ? userHooksPtr->canModifySigma() : false;
160
161}
162
163//*********
164
165// Allow for nonisotropic decays when ME's available.
166
167void PhaseSpace::decayKinematics( Event& process) {
168
169 // Identify sets of sister partons.
170 int iResEnd = 4;
171 for (int iResBeg = 5; iResBeg < process.size(); ++iResBeg) {
172 if (iResBeg <= iResEnd) continue;
173 iResEnd = iResBeg;
174 while ( iResEnd < process.size() - 1
175 && process[iResEnd + 1].mother1() == process[iResBeg].mother1()
176 && process[iResEnd + 1].mother2() == process[iResBeg].mother2() )
177 ++iResEnd;
178
179 // Check that at least one of them is a resonance.
180 bool hasRes = false;
181 for (int iRes = iResBeg; iRes <= iResEnd; ++iRes)
182 if ( !process[iRes].isFinal() ) hasRes = true;
183 if ( !hasRes ) continue;
184
185 // Evaluate matrix element and decide whether to keep kinematics.
186 double decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd);
187 if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
188 "Kinematics: negative angular weight");
189 if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
190 "Kinematics: angular weight above unity");
191 while (decWt < Rndm::flat() ) {
192
193 // Find resonances for which to redo decay angles.
194 for (int iRes = iResBeg; iRes < process.size(); ++iRes) {
195 if ( process[iRes].isFinal() ) continue;
196 int iResMother = iRes;
197 while (iResMother > iResEnd)
198 iResMother = process[iResMother].mother1();
199 if (iResMother < iResBeg) continue;
200
201 // Do decay of this mother isotropically in phase space.
202 decayKinematicsStep( process, iRes);
203
204 // End loop over resonance decay chains.
205 }
206
207 // Ready to allow new test of matrix element.
208 decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd);
209 if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
210 "Kinematics: negative angular weight");
211 if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
212 "Kinematics: angular weight above unity");
213 }
214
215 // End loop over sets of sister resonances/partons.
216 }
217
218}
219
220//*********
221
222// Reselect decay products momenta isotropically in phase space.
223// Does not redo secondary vertex position!
224
225void PhaseSpace::decayKinematicsStep( Event& process, int iRes) {
226
227 // Multiplicity and mother mass and four-momentum.
228 int i1 = process[iRes].daughter1();
229 int mult = process[iRes].daughter2() + 1 - i1;
230 double m0 = process[iRes].m();
231 Vec4 pRes = process[iRes].p();
232
233 // Description of two-body decays as simple special case.
234 if (mult == 2) {
235
236 // Products and product masses.
237 int i2 = i1 + 1;
238 double m1t = process[i1].m();
239 double m2t = process[i2].m();
240
241 // Energies and absolute momentum in the rest frame.
242 double e1 = 0.5 * (m0*m0 + m1t*m1t - m2t*m2t) / m0;
243 double e2 = 0.5 * (m0*m0 + m2t*m2t - m1t*m1t) / m0;
244 double p12 = 0.5 * sqrtpos( (m0 - m1t - m2t) * (m0 + m1t + m2t)
245 * (m0 + m1t - m2t) * (m0 - m1t + m2t) ) / m0;
246
247 // Pick isotropic angles to give three-momentum.
248 double cosTheta = 2. * Rndm::flat() - 1.;
249 double sinTheta = sqrt(1. - cosTheta*cosTheta);
250 double phi12 = 2. * M_PI * Rndm::flat();
251 double pX = p12 * sinTheta * cos(phi12);
252 double pY = p12 * sinTheta * sin(phi12);
253 double pZ = p12 * cosTheta;
254
255 // Fill four-momenta in mother rest frame and then boost to lab frame.
256 Vec4 p1( pX, pY, pZ, e1);
257 Vec4 p2( -pX, -pY, -pZ, e2);
258 p1.bst( pRes );
259 p2.bst( pRes );
260
261 // Done for two-body decay.
262 process[i1].p( p1 );
263 process[i2].p( p2 );
264 return;
265 }
266
267 // Description of three-body decays as semi-simple special case.
268 if (mult == 3) {
269
270 // Products and product masses.
271 int i2 = i1 + 1;
272 int i3 = i2 + 1;
273 double m1t = process[i1].m();
274 double m2t = process[i2].m();
275 double m3t = process[i3].m();
276 double mDiff = m0 - (m1t + m2t + m3t);
277
278 // Kinematical limits for 2+3 mass. Maximum phase-space weight.
279 double m23Min = m2t + m3t;
280 double m23Max = m0 - m1t;
281 double p1Max = 0.5 * sqrtpos( (m0 - m1t - m23Min)
282 * (m0 + m1t + m23Min) * (m0 + m1t - m23Min)
283 * (m0 - m1t + m23Min) ) / m0;
284 double p23Max = 0.5 * sqrtpos( (m23Max - m2t - m3t)
285 * (m23Max + m2t + m3t) * (m23Max + m2t - m3t)
286 * (m23Max - m2t + m3t) ) / m23Max;
287 double wtPSmax = 0.5 * p1Max * p23Max;
288
289 // Pick an intermediate mass m23 flat in the allowed range.
290 double wtPS, m23, p1Abs, p23Abs;
291 do {
292 m23 = m23Min + Rndm::flat() * mDiff;
293
294 // Translate into relative momenta and find phase-space weight.
295 p1Abs = 0.5 * sqrtpos( (m0 - m1t - m23) * (m0 + m1t + m23)
296 * (m0 + m1t - m23) * (m0 - m1t + m23) ) / m0;
297 p23Abs = 0.5 * sqrtpos( (m23 - m2t - m3t) * (m23 + m2t + m3t)
298 * (m23 + m2t - m3t) * (m23 - m2t + m3t) ) / m23;
299 wtPS = p1Abs * p23Abs;
300
301 // If rejected, try again with new invariant masses.
302 } while ( wtPS < Rndm::flat() * wtPSmax );
303
304 // Set up m23 -> m2 + m3 isotropic in its rest frame.
305 double cosTheta = 2. * Rndm::flat() - 1.;
306 double sinTheta = sqrt(1. - cosTheta*cosTheta);
307 double phi23 = 2. * M_PI * Rndm::flat();
308 double pX = p23Abs * sinTheta * cos(phi23);
309 double pY = p23Abs * sinTheta * sin(phi23);
310 double pZ = p23Abs * cosTheta;
311 double e2 = sqrt( m2t*m2t + p23Abs*p23Abs);
312 double e3 = sqrt( m3t*m3t + p23Abs*p23Abs);
313 Vec4 p2( pX, pY, pZ, e2);
314 Vec4 p3( -pX, -pY, -pZ, e3);
315
316 // Set up 0 -> 1 + 23 isotropic in its rest frame.
317 cosTheta = 2. * Rndm::flat() - 1.;
318 sinTheta = sqrt(1. - cosTheta*cosTheta);
319 phi23 = 2. * M_PI * Rndm::flat();
320 pX = p1Abs * sinTheta * cos(phi23);
321 pY = p1Abs * sinTheta * sin(phi23);
322 pZ = p1Abs * cosTheta;
323 double e1 = sqrt( m1t*m1t + p1Abs*p1Abs);
324 double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
325 Vec4 p1( pX, pY, pZ, e1);
326
327 // Boost 2 + 3 to the 0 rest frame and then boost to lab frame.
328 Vec4 p23( -pX, -pY, -pZ, e23);
329 p2.bst( p23 );
330 p3.bst( p23 );
331 p1.bst( pRes );
332 p2.bst( pRes );
333 p3.bst( pRes );
334
335 // Done for three-body decay.
336 process[i1].p( p1 );
337 process[i2].p( p2 );
338 process[i3].p( p3 );
339 return;
340 }
341
342 // Do a multibody decay using the M-generator algorithm.
343
344 // Set up masses and four-momenta in a vector, with mother in slot 0.
345 vector<double> mProd;
346 mProd.push_back( m0);
347 for (int i = i1; i <= process[iRes].daughter2(); ++i)
348 mProd.push_back( process[i].m() );
349 vector<Vec4> pProd;
350 pProd.push_back( pRes);
351
352 // Sum of daughter masses.
353 double mSum = mProd[1];
354 for (int i = 2; i <= mult; ++i) mSum += mProd[i];
355 double mDiff = m0 - mSum;
356
357 // Begin setup of intermediate invariant masses.
358 vector<double> mInv;
359 for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
360
361 // Calculate the maximum weight in the decay.
362 double wtPSmax = 1. / WTCORRECTION[mult];
363 double mMaxWT = mDiff + mProd[mult];
364 double mMinWT = 0.;
365 for (int i = mult - 1; i > 0; --i) {
366 mMaxWT += mProd[i];
367 mMinWT += mProd[i+1];
368 double mNow = mProd[i];
369 wtPSmax *= 0.5 * sqrtpos( (mMaxWT - mMinWT - mNow)
370 * (mMaxWT + mMinWT + mNow) * (mMaxWT + mMinWT - mNow)
371 * (mMaxWT - mMinWT + mNow) ) / mMaxWT;
372 }
373
374 // Begin loop to find the set of intermediate invariant masses.
375 vector<double> rndmOrd;
376 double wtPS;
377 do {
378 wtPS = 1.;
379
380 // Find and order random numbers in descending order.
381 rndmOrd.resize(0);
382 rndmOrd.push_back(1.);
383 for (int i = 1; i < mult - 1; ++i) {
384 double rndm = Rndm::flat();
385 rndmOrd.push_back(rndm);
386 for (int j = i - 1; j > 0; --j) {
387 if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
388 else break;
389 }
390 }
391 rndmOrd.push_back(0.);
392
393 // Translate into intermediate masses and find weight.
394 for (int i = mult - 1; i > 0; --i) {
395 mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
396 wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
397 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
398 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
399 }
400
401 // If rejected, try again with new invariant masses.
402 } while ( wtPS < Rndm::flat() * wtPSmax );
403
404 // Perform two-particle decays in the respective rest frame.
405 vector<Vec4> pInv;
406 pInv.resize(mult + 1);
407 for (int i = 1; i < mult; ++i) {
408 double p12 = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
409 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
410 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
411
412 // Isotropic angles give three-momentum.
413 double cosTheta = 2. * Rndm::flat() - 1.;
414 double sinTheta = sqrt(1. - cosTheta*cosTheta);
415 double phiLoc = 2. * M_PI * Rndm::flat();
416 double pX = p12 * sinTheta * cos(phiLoc);
417 double pY = p12 * sinTheta * sin(phiLoc);
418 double pZ = p12 * cosTheta;
419
420 // Calculate energies, fill four-momenta.
421 double eHad = sqrt( mProd[i]*mProd[i] + p12*p12);
422 double eInv = sqrt( mInv[i+1]*mInv[i+1] + p12*p12);
423 pProd.push_back( Vec4( pX, pY, pZ, eHad) );
424 pInv[i+1].p( -pX, -pY, -pZ, eInv);
425 }
426 pProd.push_back( pInv[mult] );
427
428 // Boost decay products to the mother rest frame and on to lab frame.
429 pInv[1] = pProd[0];
430 for (int iFrame = mult - 1; iFrame > 0; --iFrame)
431 for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);
432
433 // Done for multibody decay.
434 for (int i = 1; i < mult; ++i)
435 process[i1 + i - 1].p( pProd[i] );
436 return;
437
438}
439
440//*********
441
442// Determine how 3-body phase space should be sampled.
443
444void PhaseSpace::setup3Body() {
445
446 // Check for massive t-channel propagator particles.
447 int idTchan1 = abs( sigmaProcessPtr->idTchan1() );
448 int idTchan2 = abs( sigmaProcessPtr->idTchan2() );
449 mTchan1 = (idTchan1 == 0) ? 0. : ParticleDataTable::m0(idTchan1);
450 mTchan2 = (idTchan2 == 0) ? 0. : ParticleDataTable::m0(idTchan2);
451 sTchan1 = mTchan1 * mTchan1;
452 sTchan2 = mTchan2 * mTchan2;
453
454 // Find coefficients of different pT2 selection terms. Mirror choice.
455 frac3Pow1 = sigmaProcessPtr->tChanFracPow1();
456 frac3Pow2 = sigmaProcessPtr->tChanFracPow2();
457 frac3Flat = 1. - frac3Pow1 - frac3Pow2;
458 useMirrorWeight = sigmaProcessPtr->useMirrorWeight();
459
460}
461
462//*********
463
464// Determine how phase space should be sampled.
465
466bool PhaseSpace::setupSampling123(bool is2, bool is3, ostream& os) {
467
468 // Optional printout.
469 if (showSearch) os << "\n PYTHIA Optimization printout for "
470 << sigmaProcessPtr->name() << "\n \n" << scientific << setprecision(3);
471
472 // Check that open range in tau (+ set tauMin, tauMax).
473 if (!limitTau(is2, is3)) return false;
474
475 // Reset coefficients and matrices of equation system to solve.
476 int binTau[8], binY[8], binZ[8];
477 double vecTau[8], matTau[8][8], vecY[8], matY[8][8], vecZ[8], matZ[8][8];
478 for (int i = 0; i < 8; ++i) {
479 tauCoef[i] = 0.;
480 yCoef[i] = 0.;
481 zCoef[i] = 0.;
482 binTau[i] = 0;
483 binY[i] = 0;
484 binZ[i] = 0;
485 vecTau[i] = 0.;
486 vecY[i] = 0.;
487 vecZ[i] = 0.;
488 for (int j = 0; j < 8; ++j) {
489 matTau[i][j] = 0.;
490 matY[i][j] = 0.;
491 matZ[i][j] = 0.;
492 }
493 }
494 sigmaMx = 0.;
495 sigmaNeg = 0.;
496
497 // Number of used coefficients/points for each dimension: tau, y, c.
498 nTau = (hasPointLeptons) ? 1 : 2;
499 nY = (hasPointLeptons) ? 1 : 3;
500 nZ = (is2) ? 5 : 1;
501
502 // Identify if any resonances contribute in s-channel.
503 idResA = sigmaProcessPtr->resonanceA();
504 if (idResA != 0) {
505 mResA = ParticleDataTable::m0(idResA);
506 GammaResA = ParticleDataTable::mWidth(idResA);
507 if (mHatMin > mResA + WIDTHMARGIN * GammaResA || (mHatMax > 0.
508 && mHatMax < mResA - WIDTHMARGIN * GammaResA) ) idResA = 0;
509 }
510 idResB = sigmaProcessPtr->resonanceB();
511 if (idResB != 0) {
512 mResB = ParticleDataTable::m0(idResB);
513 GammaResB = ParticleDataTable::mWidth(idResB);
514 if (mHatMin > mResB + WIDTHMARGIN * GammaResB || (mHatMax > 0.
515 && mHatMax < mResB - WIDTHMARGIN * GammaResB) ) idResB = 0;
516 }
517 if (idResA == 0 && idResB != 0) {
518 idResA = idResB;
519 mResA = mResB;
520 GammaResA = GammaResB;
521 idResB = 0;
522 }
523
524 // More sampling in tau if resonances in s-channel.
525 if (idResA !=0 && !hasPointLeptons) {
526 nTau += 2;
527 tauResA = mResA * mResA / s;
528 widResA = mResA * GammaResA / s;
529 }
530 if (idResB != 0 && !hasPointLeptons) {
531 nTau += 2;
532 tauResB = mResB * mResB / s;
533 widResB = mResB * GammaResB / s;
534 }
535
536 // More sampling in tau and y if incoming lepton beams.
537 if (hasLeptonBeams && !hasPointLeptons) {
538 ++nTau;
539 nY += 2;
540 }
541
542 // Special case when both resonances have same mass.
543 sameResMass = false;
544 if (idResB != 0 && abs(mResA - mResB) < SAMEMASS * (GammaResA + GammaResB))
545 sameResMass = true;
546
547 // Default z value and weight required for 2 -> 1. Number of dimensions.
548 z = 0.;
549 wtZ = 1.;
550 int nVar = (is2) ? 3 : 2;
551
552 // Initial values, to be modified later.
553 tauCoef[0] = 1.;
554 yCoef[0] = 0.5;
555 yCoef[1] = 0.5;
556 zCoef[0] = 1.;
557
558 // Step through grid in tau. Set limits on y and z generation.
559 for (int iTau = 0; iTau < nTau; ++iTau) {
560 double posTau = 0.5;
561 if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6;
562 selectTau( iTau, posTau, is2);
563 if (!limitY()) continue;
564 if (is2 && !limitZ()) continue;
565
566 // Step through grids in y and z.
567 for (int iY = 0; iY < nY; ++iY) {
568 selectY( iY, 0.5);
569 for (int iZ = 0; iZ < nZ; ++iZ) {
570 if (is2) selectZ( iZ, 0.5);
571 double sigmaTmp = 0.;
572
573 // 2 -> 1: calculate cross section, weighted by phase-space volume.
574 if (!is2 && !is3) {
575 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
576 sigmaTmp = sigmaProcessPtr->sigmaPDF();
577 sigmaTmp *= wtTau * wtY;
578
579 // 2 -> 2: calculate cross section, weighted by phase-space volume
580 // and Breit-Wigners for masses
581 } else if (is2) {
582 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
583 runBW3H, runBW4H);
584 sigmaTmp = sigmaProcessPtr->sigmaPDF();
585 sigmaTmp *= wtTau * wtY * wtZ * wtBW;
586
587 // 2 -> 3: repeat internal 3-body phase space several times and
588 // keep maximal cross section, weighted by phase-space volume
589 // and Breit-Wigners for masses
590 } else if (is3) {
591 for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
592 if (!select3Body()) continue;
593 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
594 m3, m4, m5, runBW3H, runBW4H, runBW5H);
595 double sigmaTry = sigmaProcessPtr->sigmaPDF();
596 sigmaTry *= wtTau * wtY * wt3Body * wtBW;
597 if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
598 }
599 }
600
601 // Allow possibility for user to modify cross section. (3body??)
602 if (canModifySigma) sigmaTmp
603 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
604
605 // Check if current maximum exceeded.
606 if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp;
607
608 // Optional printout. Protect against negative cross sections.
609 if (showSearch) os << " tau =" << setw(11) << tau << " y ="
610 << setw(11) << y << " z =" << setw(11) << z
611 << " sigma =" << setw(11) << sigmaTmp << "\n";
612 if (sigmaTmp < 0.) sigmaTmp = 0.;
613
614 // Sum up tau cross-section pieces in points used.
615 if (!hasPointLeptons) {
616 binTau[iTau] += 1;
617 vecTau[iTau] += sigmaTmp;
618 matTau[iTau][0] += 1. / intTau0;
619 matTau[iTau][1] += (1. / intTau1) / tau;
620 if (idResA != 0) {
621 matTau[iTau][2] += (1. / intTau2) / (tau + tauResA);
622 matTau[iTau][3] += (1. / intTau3)
623 * tau / ( pow2(tau - tauResA) + pow2(widResA) );
624 }
625 if (idResB != 0) {
626 matTau[iTau][4] += (1. / intTau4) / (tau + tauResB);
627 matTau[iTau][5] += (1. / intTau5)
628 * tau / ( pow2(tau - tauResB) + pow2(widResB) );
629 }
630 if (hasLeptonBeams) matTau[iTau][nTau - 1] += (1. / intTau6)
631 * tau / max( LEPTONTAUMIN, 1. - tau);
632 }
633
634 // Sum up y cross-section pieces in points used.
635 if (!hasPointLeptons) {
636 binY[iY] += 1;
637 vecY[iY] += sigmaTmp;
638 matY[iY][0] += (yMax / intY01) * (y + yMax);
639 matY[iY][1] += (yMax / intY01) * (yMax - y);
640 matY[iY][2] += (yMax / intY2) / cosh(y);
641 if (hasLeptonBeams) {
642 matY[iY][3] += (yMax / intY34)
643 / max( LEPTONXMIN, 1. - exp( y - yMax) );
644 matY[iY][4] += (yMax / intY34)
645 / max( LEPTONXMIN, 1. - exp(-y - yMax) );
646 }
647 }
648
649 // Integrals over z expressions at tauMax, to be used below.
650 if (is2) {
651 double p2AbsMax = 0.25 * (pow2(tauMax * s - s3 - s4)
652 - 4. * s3 * s4) / (tauMax * s);
653 double zMaxMax = sqrtpos( 1. - pT2HatMin / p2AbsMax );
654 double zPosMaxMax = max(ratio34, unity34 + zMaxMax);
655 double zNegMaxMax = max(ratio34, unity34 - zMaxMax);
656 double intZ0Max = 2. * zMaxMax;
657 double intZ12Max = log( zPosMaxMax / zNegMaxMax);
658 double intZ34Max = 1. / zNegMaxMax - 1. / zPosMaxMax;
659
660 // Sum up z cross-section pieces in points used.
661 binZ[iZ] += 1;
662 vecZ[iZ] += sigmaTmp;
663 matZ[iZ][0] += 1.;
664 matZ[iZ][1] += (intZ0Max / intZ12Max) / zNeg;
665 matZ[iZ][2] += (intZ0Max / intZ12Max) / zPos;
666 matZ[iZ][3] += (intZ0Max / intZ34Max) / pow2(zNeg);
667 matZ[iZ][4] += (intZ0Max / intZ34Max) / pow2(zPos);
668 }
669
670 // End of loops over phase space points.
671 }
672 }
673 }
674
675 // Fail if no non-vanishing cross sections.
676 if (sigmaMx <= 0.) {
677 sigmaMx = 0.;
678 return false;
679 }
680
681 // Solve respective equation system for better phase space coefficients.
682 if (!hasPointLeptons) solveSys( nTau, binTau, vecTau, matTau, tauCoef);
683 if (!hasPointLeptons) solveSys( nY, binY, vecY, matY, yCoef);
684 if (is2) solveSys( nZ, binZ, vecZ, matZ, zCoef);
685 if (showSearch) os << "\n";
686
687 // Provide cumulative sum of coefficients.
688 tauCoefSum[0] = tauCoef[0];
689 yCoefSum[0] = yCoef[0];
690 zCoefSum[0] = zCoef[0];
691 for (int i = 1; i < 8; ++ i) {
692 tauCoefSum[i] = tauCoefSum[i - 1] + tauCoef[i];
693 yCoefSum[i] = yCoefSum[i - 1] + yCoef[i];
694 zCoefSum[i] = zCoefSum[i - 1] + zCoef[i];
695 }
696 // The last element should be > 1 to be on safe side in selection below.
697 tauCoefSum[nTau - 1] = 2.;
698 yCoefSum[nY - 1] = 2.;
699 zCoefSum[nZ - 1] = 2.;
700
701
702 // Begin find two most promising maxima among same points as before.
703 int iMaxTau[NMAXTRY + 2], iMaxY[NMAXTRY + 2], iMaxZ[NMAXTRY + 2];
704 double sigMax[NMAXTRY + 2];
705 int nMax = 0;
706
707 // Scan same grid as before in tau, y, z.
708 for (int iTau = 0; iTau < nTau; ++iTau) {
709 double posTau = 0.5;
710 if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6;
711 selectTau( iTau, posTau, is2);
712 if (!limitY()) continue;
713 if (is2 && !limitZ()) continue;
714 for (int iY = 0; iY < nY; ++iY) {
715 selectY( iY, 0.5);
716 for (int iZ = 0; iZ < nZ; ++iZ) {
717 if (is2) selectZ( iZ, 0.5);
718 double sigmaTmp = 0.;
719
720 // 2 -> 1: calculate cross section, weighted by phase-space volume.
721 if (!is2 && !is3) {
722 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
723 sigmaTmp = sigmaProcessPtr->sigmaPDF();
724 sigmaTmp *= wtTau * wtY;
725
726 // 2 -> 2: calculate cross section, weighted by phase-space volume
727 // and Breit-Wigners for masses
728 } else if (is2) {
729 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
730 runBW3H, runBW4H);
731 sigmaTmp = sigmaProcessPtr->sigmaPDF();
732 sigmaTmp *= wtTau * wtY * wtZ * wtBW;
733
734 // 2 -> 3: repeat internal 3-body phase space several times and
735 // keep maximal cross section, weighted by phase-space volume
736 // and Breit-Wigners for masses
737 } else if (is3) {
738 for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
739 if (!select3Body()) continue;
740 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
741 m3, m4, m5, runBW3H, runBW4H, runBW5H);
742 double sigmaTry = sigmaProcessPtr->sigmaPDF();
743 sigmaTry *= wtTau * wtY * wt3Body * wtBW;
744 if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
745 }
746 }
747
748 // Allow possibility for user to modify cross section. (3body??)
749 if (canModifySigma) sigmaTmp
750 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
751
752 // Optional printout. Protect against negative cross section.
753 if (showSearch) os << " tau =" << setw(11) << tau << " y ="
754 << setw(11) << y << " z =" << setw(11) << z
755 << " sigma =" << setw(11) << sigmaTmp << "\n";
756 if (sigmaTmp < 0.) sigmaTmp = 0.;
757
758 // Check that point is not simply mirror of already found one.
759 bool mirrorPoint = false;
760 for (int iMove = 0; iMove < nMax; ++iMove)
761 if (abs(sigmaTmp - sigMax[iMove]) < SAMESIGMA
762 * (sigmaTmp + sigMax[iMove])) mirrorPoint = true;
763
764 // Add to or insert in maximum list. Only first two count.
765 if (!mirrorPoint) {
766 int iInsert = 0;
767 for (int iMove = nMax - 1; iMove >= -1; --iMove) {
768 iInsert = iMove + 1;
769 if (iInsert == 0 || sigmaTmp < sigMax[iMove]) break;
770 iMaxTau[iMove + 1] = iMaxTau[iMove];
771 iMaxY[iMove + 1] = iMaxY[iMove];
772 iMaxZ[iMove + 1] = iMaxZ[iMove];
773 sigMax[iMove + 1] = sigMax[iMove];
774 }
775 iMaxTau[iInsert] = iTau;
776 iMaxY[iInsert] = iY;
777 iMaxZ[iInsert] = iZ;
778 sigMax[iInsert] = sigmaTmp;
779 if (nMax < NMAXTRY) ++nMax;
780 }
781
782 // Found two most promising maxima.
783 }
784 }
785 }
786 if (showSearch) os << "\n";
787
788 // Read out starting position for search.
789 sigmaMx = sigMax[0];
790 int beginVar = (hasPointLeptons) ? 2 : 0;
791 for (int iMax = 0; iMax < nMax; ++iMax) {
792 int iTau = iMaxTau[iMax];
793 int iY = iMaxY[iMax];
794 int iZ = iMaxZ[iMax];
795 double tauVal = 0.5;
796 double yVal = 0.5;
797 double zVal = 0.5;
798 int iGrid;
799 double varVal, varNew, deltaVar, marginVar, sigGrid[3];
800
801 // Starting point and step size in parameter space.
802 for (int iRepeat = 0; iRepeat < 2; ++iRepeat) {
803 // Run through (possibly a subset of) tau, y and z.
804 for (int iVar = beginVar; iVar < nVar; ++iVar) {
805 if (iVar == 0) varVal = tauVal;
806 else if (iVar == 1) varVal = yVal;
807 else varVal = zVal;
808 deltaVar = (iRepeat == 0) ? 0.1
809 : max( 0.01, min( 0.05, min( varVal - 0.02, 0.98 - varVal) ) );
810 marginVar = (iRepeat == 0) ? 0.02 : 0.002;
811 int moveStart = (iRepeat == 0 && iVar == 0) ? 0 : 1;
812 for (int move = moveStart; move < 9; ++move) {
813
814 // Define new parameter-space point by step in one dimension.
815 if (move == 0) {
816 iGrid = 1;
817 varNew = varVal;
818 } else if (move == 1) {
819 iGrid = 2;
820 varNew = varVal + deltaVar;
821 } else if (move == 2) {
822 iGrid = 0;
823 varNew = varVal - deltaVar;
824 } else if (sigGrid[2] >= max( sigGrid[0], sigGrid[1])
825 && varVal + 2. * deltaVar < 1. - marginVar) {
826 varVal += deltaVar;
827 sigGrid[0] = sigGrid[1];
828 sigGrid[1] = sigGrid[2];
829 iGrid = 2;
830 varNew = varVal + deltaVar;
831 } else if (sigGrid[0] >= max( sigGrid[1], sigGrid[2])
832 && varVal - 2. * deltaVar > marginVar) {
833 varVal -= deltaVar;
834 sigGrid[2] = sigGrid[1];
835 sigGrid[1] = sigGrid[0];
836 iGrid = 0;
837 varNew = varVal - deltaVar;
838 } else if (sigGrid[2] >= sigGrid[0]) {
839 deltaVar *= 0.5;
840 varVal += deltaVar;
841 sigGrid[0] = sigGrid[1];
842 iGrid = 1;
843 varNew = varVal;
844 } else {
845 deltaVar *= 0.5;
846 varVal -= deltaVar;
847 sigGrid[2] = sigGrid[1];
848 iGrid = 1;
849 varNew = varVal;
850 }
851
852 // Convert to relevant variables and find derived new limits.
853 bool insideLimits = true;
854 if (iVar == 0) {
855 tauVal = varNew;
856 selectTau( iTau, tauVal, is2);
857 if (!limitY()) insideLimits = false;
858 if (is2 && !limitZ()) insideLimits = false;
859 if (insideLimits) {
860 selectY( iY, yVal);
861 if (is2) selectZ( iZ, zVal);
862 }
863 } else if (iVar == 1) {
864 yVal = varNew;
865 selectY( iY, yVal);
866 } else if (iVar == 2) {
867 zVal = varNew;
868 selectZ( iZ, zVal);
869 }
870
871 // Evaluate cross-section.
872 double sigmaTmp = 0.;
873 if (insideLimits) {
874
875 // 2 -> 1: calculate cross section, weighted by phase-space volume.
876 if (!is2 && !is3) {
877 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
878 sigmaTmp = sigmaProcessPtr->sigmaPDF();
879 sigmaTmp *= wtTau * wtY;
880
881 // 2 -> 2: calculate cross section, weighted by phase-space volume
882 // and Breit-Wigners for masses
883 } else if (is2) {
884 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
885 runBW3H, runBW4H);
886 sigmaTmp = sigmaProcessPtr->sigmaPDF();
887 sigmaTmp *= wtTau * wtY * wtZ * wtBW;
888
889 // 2 -> 3: repeat internal 3-body phase space several times and
890 // keep maximal cross section, weighted by phase-space volume
891 // and Breit-Wigners for masses
892 } else if (is3) {
893 for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
894 if (!select3Body()) continue;
895 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
896 m3, m4, m5, runBW3H, runBW4H, runBW5H);
897 double sigmaTry = sigmaProcessPtr->sigmaPDF();
898 sigmaTry *= wtTau * wtY * wt3Body * wtBW;
899 if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
900 }
901 }
902
903 // Allow possibility for user to modify cross section.
904 if (canModifySigma) sigmaTmp
905 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
906
907 // Optional printout. Protect against negative cross section.
908 if (showSearch) os << " tau =" << setw(11) << tau << " y ="
909 << setw(11) << y << " z =" << setw(11) << z
910 << " sigma =" << setw(11) << sigmaTmp << "\n";
911 if (sigmaTmp < 0.) sigmaTmp = 0.;
912 }
913
914 // Save new maximum. Final maximum.
915 sigGrid[iGrid] = sigmaTmp;
916 if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp;
917 }
918 }
919 }
920 }
921 sigmaMx *= SAFETYMARGIN;
922
923 // Optional printout.
924 if (showSearch) os << "\n Final maximum = " << setw(11) << sigmaMx << endl;
925
926 // Done.
927 return true;
928}
929
930//*********
931
932// Select a trial kinematics phase space point.
933// Note: by In is meant the integral over the quantity multiplying
934// coefficient cn. The sum of cn is normalized to unity.
935
936bool PhaseSpace::trialKin123(bool is2, bool is3, bool inEvent, ostream& os) {
937
938 // Allow for possibility that energy varies from event to event.
939 if (doEnergySpread) {
940 eCM = infoPtr->eCM();
941 s = eCM * eCM;
942
943 // Find shifted tauRes values.
944 if (idResA !=0 && !hasPointLeptons) {
945 tauResA = mResA * mResA / s;
946 widResA = mResA * GammaResA / s;
947 }
948 if (idResB != 0 && !hasPointLeptons) {
949 tauResB = mResB * mResB / s;
950 widResB = mResB * GammaResB / s;
951 }
952 }
953
954 // Choose tau according to h1(tau)/tau, where
955 // h1(tau) = c0/I0 + (c1/I1) * 1/tau
956 // + (c2/I2) / (tau + tauResA)
957 // + (c3/I3) * tau / ((tau - tauResA)^2 + widResA^2)
958 // + (c4/I4) / (tau + tauResB)
959 // + (c5/I5) * tau / ((tau - tauResB)^2 + widResB^2)
960 // + (c6/I6) * tau / (1 - tau).
961 if (!limitTau(is2, is3)) return false;
962 int iTau = 0;
963 if (!hasPointLeptons) {
964 double rTau = Rndm::flat();
965 while (rTau > tauCoefSum[iTau]) ++iTau;
966 }
967 selectTau( iTau, Rndm::flat(), is2);
968
969 // Choose y according to h2(y), where
970 // h2(y) = (c0/I0) * (y-ymin) + (c1/I1) * (ymax-y)
971 // + (c2/I2) * 1/cosh(y) + (c3/I3) * 1 / (1 - exp(y-ymax))
972 // + (c4/I4) * 1 / (1 - exp(ymin-y)).
973 if (!limitY()) return false;
974 int iY = 0;
975 if (!hasPointLeptons) {
976 double rY = Rndm::flat();
977 while (rY > yCoefSum[iY]) ++iY;
978 }
979 selectY( iY, Rndm::flat());
980
981 // Choose z = cos(thetaHat) according to h3(z), where
982 // h3(z) = c0/I0 + (c1/I1) * 1/(A - z) + (c2/I2) * 1/(A + z)
983 // + (c3/I3) * 1/(A - z)^2 + (c4/I4) * 1/(A + z)^2,
984 // where A = 1 + 2*(m3*m4/sH)^2 (= 1 for massless products).
985 if (is2) {
986 if (!limitZ()) return false;
987 int iZ = 0;
988 double rZ = Rndm::flat();
989 while (rZ > zCoefSum[iZ]) ++iZ;
990 selectZ( iZ, Rndm::flat());
991 }
992
993 // 2 -> 1: calculate cross section, weighted by phase-space volume.
994 if (!is2 && !is3) {
995 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
996 sigmaNw = sigmaProcessPtr->sigmaPDF();
997 sigmaNw *= wtTau * wtY;
998
999 // 2 -> 2: calculate cross section, weighted by phase-space volume
1000 // and Breit-Wigners for masses
1001 } else if (is2) {
1002 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, runBW3H, runBW4H);
1003 sigmaNw = sigmaProcessPtr->sigmaPDF();
1004 sigmaNw *= wtTau * wtY * wtZ * wtBW;
1005
1006 // 2 -> 3: also sample internal 3-body phase, weighted by
1007 // 2 -> 1 phase-space volume and Breit-Wigners for masses
1008 } else if (is3) {
1009 if (!select3Body()) sigmaNw = 0.;
1010 else {
1011 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
1012 m3, m4, m5, runBW3H, runBW4H, runBW5H);
1013 sigmaNw = sigmaProcessPtr->sigmaPDF();
1014 sigmaNw *= wtTau * wtY * wt3Body * wtBW;
1015 }
1016 }
1017
1018 // Allow possibility for user to modify cross section.
1019 if (canModifySigma) sigmaNw
1020 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, inEvent);
1021
1022 // Check if maximum violated.
1023 newSigmaMx = false;
1024 if (sigmaNw > sigmaMx) {
1025 infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin: "
1026 "maximum for cross section violated");
1027 double violFact = SAFETYMARGIN * sigmaNw / sigmaMx;
1028 sigmaMx = SAFETYMARGIN * sigmaNw;
1029 newSigmaMx = true;
1030
1031 // Optional printout of (all) violations.
1032 if (showViolation) {
1033 if (violFact < 9.99) os << fixed;
1034 else os << scientific;
1035 os << " PYTHIA Maximum for " << sigmaProcessPtr->name()
1036 << " increased by factor " << setprecision(3) << violFact
1037 << " to " << scientific << sigmaMx << endl;
1038 }
1039 }
1040
1041 // Check if negative cross section.
1042 if (sigmaNw < sigmaNeg) {
1043 infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin:"
1044 " negative cross section set 0", "for " + sigmaProcessPtr->name() );
1045 sigmaNeg = sigmaNw;
1046
1047 // Optional printout of (all) violations.
1048 if (showViolation) os << " PYTHIA Negative minimum for "
1049 << sigmaProcessPtr->name() << " changed to " << scientific
1050 << setprecision(3) << sigmaNeg << endl;
1051 }
1052 if (sigmaNw < 0.) sigmaNw = 0.;
1053
1054 // Done.
1055 return true;
1056}
1057
1058//*********
1059
1060// Find range of allowed tau values.
1061
1062bool PhaseSpace::limitTau(bool is2, bool is3) {
1063
1064 // Trivial reply for unresolved lepton beams.
1065 if (hasPointLeptons) {
1066 tauMin = 1.;
1067 tauMax = 1.;
1068 return true;
1069 }
1070
1071 // Requirements from allowed mHat range.
1072 tauMin = sHatMin / s;
1073 tauMax = (mHatMax < mHatMin) ? 1. : min( 1., sHatMax / s);
1074
1075 // Requirements from allowed pT range and masses.
1076 if (is2 || is3) {
1077 double mT3Min = sqrt(s3 + pT2HatMin);
1078 double mT4Min = sqrt(s4 + pT2HatMin);
1079 double mT5Min = (is3) ? sqrt(s5 + pT2HatMin) : 0.;
1080 tauMin = max( tauMin, pow2(mT3Min + mT4Min + mT5Min) / s);
1081 }
1082
1083 // Check that there is an open range.
1084 return (tauMax > tauMin);
1085}
1086
1087//*********
1088
1089// Find range of allowed y values.
1090
1091bool PhaseSpace::limitY() {
1092
1093 // Trivial reply for unresolved lepton beams.
1094 if (hasPointLeptons) {
1095 yMax = 1.;
1096 return true;
1097 }
1098
1099 // Requirements from selected tau value.
1100 yMax = -0.5 * log(tau);
1101
1102 // For lepton beams requirements from cutoff for f_e^e.
1103 double yMaxMargin = (hasLeptonBeams) ? yMax + LEPTONXLOGMAX : yMax;
1104
1105 // Check that there is an open range.
1106 return (yMaxMargin > 0.);
1107}
1108
1109//*********
1110
1111// Find range of allowed z = cos(theta) values.
1112
1113bool PhaseSpace::limitZ() {
1114
1115 // Default limits.
1116 zMin = 0.;
1117 zMax = 1.;
1118
1119 // Requirements from pTHat limits.
1120 zMax = sqrtpos( 1. - pT2HatMin / p2Abs );
1121 if (pTHatMax > pTHatMin) zMin = sqrtpos( 1. - pT2HatMax / p2Abs );
1122
1123 // Check that there is an open range.
1124 return (zMax > zMin);
1125}
1126
1127//*********
1128
1129// Select tau according to a choice of shapes.
1130
1131void PhaseSpace::selectTau(int iTau, double tauVal, bool is2) {
1132
1133 // Trivial reply for unresolved lepton beams.
1134 if (hasPointLeptons) {
1135 tau = 1.;
1136 wtTau = 1.;
1137 sH = s;
1138 mHat = sqrt(sH);
1139 if (is2) {
1140 p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1141 pAbs = sqrtpos( p2Abs );
1142 }
1143 return;
1144 }
1145
1146 // Contributions from s-channel resonances.
1147 double tRatA = 0.;
1148 double aLowA = 0.;
1149 double aUppA = 0.;
1150 if (idResA !=0) {
1151 tRatA = ((tauResA + tauMax) / (tauResA + tauMin)) * (tauMin / tauMax);
1152 aLowA = atan( (tauMin - tauResA) / widResA);
1153 aUppA = atan( (tauMax - tauResA) / widResA);
1154 }
1155 double tRatB = 0.;
1156 double aLowB = 0.;
1157 double aUppB = 0.;
1158 if (idResB != 0) {
1159 tRatB = ((tauResB + tauMax) / (tauResB + tauMin)) * (tauMin / tauMax);
1160 aLowB = atan( (tauMin - tauResB) / widResB);
1161 aUppB = atan( (tauMax - tauResB) / widResB);
1162 }
1163
1164 // Contributions from 1 / (1 - tau) for lepton beams.
1165 double aLowT = 0.;
1166 double aUppT = 0.;
1167 if (hasLeptonBeams) {
1168 aLowT = log( max( LEPTONTAUMIN, 1. - tauMin) );
1169 aUppT = log( max( LEPTONTAUMIN, 1. - tauMax) );
1170 intTau6 = aLowT - aUppT;
1171 }
1172
1173 // Select according to 1/tau or 1/tau^2.
1174 if (iTau == 0) tau = tauMin * pow( tauMax / tauMin, tauVal);
1175 else if (iTau == 1) tau = tauMax * tauMin
1176 / (tauMin + (tauMax - tauMin) * tauVal);
1177
1178 // Select according to 1 / (1 - tau) for lepton beams.
1179 else if (hasLeptonBeams && iTau == nTau - 1)
1180 tau = 1. - exp( aUppT + intTau6 * tauVal );
1181
1182 // Select according to 1 / (tau * (tau + tauRes)) or
1183 // 1 / ((tau - tauRes)^2 + widRes^2) for resonances A and B.
1184 else if (iTau == 2) tau = tauResA * tauMin
1185 / ((tauResA + tauMin) * pow( tRatA, tauVal) - tauMin);
1186 else if (iTau == 3) tau = tauResA + widResA
1187 * tan( aLowA + (aUppA - aLowA) * tauVal);
1188 else if (iTau == 4) tau = tauResB * tauMin
1189 / ((tauResB + tauMin) * pow( tRatB, tauVal) - tauMin);
1190 else if (iTau == 5) tau = tauResB + widResB
1191 * tan( aLowB + (aUppB - aLowB) * tauVal);
1192
1193 // Phase space weight in tau.
1194 intTau0 = log( tauMax / tauMin);
1195 intTau1 = (tauMax - tauMin) / (tauMax * tauMin);
1196 double invWtTau = (tauCoef[0] / intTau0) + (tauCoef[1] / intTau1) / tau;
1197 if (idResA != 0) {
1198 intTau2 = -log(tRatA) / tauResA;
1199 intTau3 = (aUppA - aLowA) / widResA;
1200 invWtTau += (tauCoef[2] / intTau2) / (tau + tauResA)
1201 + (tauCoef[3] / intTau3) * tau / ( pow2(tau - tauResA) + pow2(widResA) );
1202 }
1203 if (idResB != 0) {
1204 intTau4 = -log(tRatB) / tauResB;
1205 intTau5 = (aUppB - aLowB) / widResB;
1206 invWtTau += (tauCoef[4] / intTau4) / (tau + tauResB)
1207 + (tauCoef[5] / intTau5) * tau / ( pow2(tau - tauResB) + pow2(widResB) );
1208 }
1209 if (hasLeptonBeams)
1210 invWtTau += (tauCoef[nTau - 1] / intTau6)
1211 * tau / max( LEPTONTAUMIN, 1. - tau);
1212 wtTau = 1. / invWtTau;
1213
1214 // Calculate sHat and absolute momentum of outgoing partons.
1215 sH = tau * s;
1216 mHat = sqrt(sH);
1217 if (is2) {
1218 p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1219 pAbs = sqrtpos( p2Abs );
1220 }
1221
1222}
1223
1224//*********
1225
1226// Select y according to a choice of shapes.
1227
1228void PhaseSpace::selectY(int iY, double yVal) {
1229
1230 // Trivial reply for unresolved lepton beams.
1231 if (hasPointLeptons) {
1232 y = 0.;
1233 wtY = 1.;
1234 x1H = 1.;
1235 x2H = 1.;
1236 return;
1237 }
1238
1239 // Standard expressions used below.
1240 double atanMax = atan( exp(yMax) );
1241 double atanMin = atan( exp(-yMax) );
1242 double aUppY = (hasLeptonBeams)
1243 ? log( max( LEPTONXMIN, LEPTONXMAX / tau - 1. ) ) : 0.;
1244 double aLowY = LEPTONXLOGMIN;
1245
1246 // y - y_min or mirrored y_max - y.
1247 if (iY <= 1) y = yMax * (2. * sqrt(yVal) - 1.);
1248
1249 // 1 / cosh(y).
1250 else if (iY == 2)
1251 y = log( tan( atanMin + (atanMax - atanMin) * yVal ) );
1252
1253 // 1 / (1 - exp(y - y_max)) or mirrored 1 / (1 - exp(y_min - y)).
1254 else y = yMax - log( 1. + exp(aLowY + (aUppY - aLowY) * yVal) );
1255
1256 // Mirror two cases.
1257 if (iY == 1 || iY == 4) y = -y;
1258
1259 // Phase space integral in y.
1260 intY01 = 0.5 * pow2(2. * yMax);
1261 intY2 = 2. * (atanMax - atanMin);
1262 intY34 = aUppY - aLowY;
1263 double invWtY = (yCoef[0] / intY01) * (y + yMax)
1264 + (yCoef[1] / intY01) * (yMax - y) + (yCoef[2] / intY2) / cosh(y);
1265 if (hasLeptonBeams) invWtY
1266 += (yCoef[3] / intY34) / max( LEPTONXMIN, 1. - exp( y - yMax) )
1267 + (yCoef[4] / intY34) / max( LEPTONXMIN, 1. - exp(-y - yMax) );
1268 wtY = 1. / invWtY;
1269
1270 // Calculate x1 and x2.
1271 x1H = sqrt(tau) * exp(y);
1272 x2H = sqrt(tau) * exp(-y);
1273}
1274
1275//*********
1276
1277// Select z = cos(theta) according to a choice of shapes.
1278// The selection is split in the positive- and negative-z regions,
1279// since a pTmax cut can remove the region around z = 0.
1280
1281void PhaseSpace::selectZ(int iZ, double zVal) {
1282
1283 // Mass-dependent dampening of pT -> 0 limit.
1284 ratio34 = max(TINY, 2. * s3 * s4 / pow2(sH));
1285 unity34 = 1. + ratio34;
1286 double ratiopT2 = 2. * pT2HatMin / max( SHATMINZ, sH);
1287 if (ratiopT2 < PT2RATMINZ) ratio34 = max( ratio34, ratiopT2);
1288
1289 // Common expressions in z limits.
1290 double zPosMax = max(ratio34, unity34 + zMax);
1291 double zNegMax = max(ratio34, unity34 - zMax);
1292 double zPosMin = max(ratio34, unity34 + zMin);
1293 double zNegMin = max(ratio34, unity34 - zMin);
1294
1295 // Flat in z.
1296 if (iZ == 0) {
1297 if (zVal < 0.5) z = -(zMax + (zMin - zMax) * 2. * zVal);
1298 else z = zMin + (zMax - zMin) * (2. * zVal - 1.);
1299
1300 // 1 / (unity34 - z).
1301 } else if (iZ == 1) {
1302 double areaNeg = log(zPosMax / zPosMin);
1303 double areaPos = log(zNegMin / zNegMax);
1304 double area = areaNeg + areaPos;
1305 if (zVal * area < areaNeg) {
1306 double zValMod = zVal * area / areaNeg;
1307 z = unity34 - zPosMax * pow(zPosMin / zPosMax, zValMod);
1308 } else {
1309 double zValMod = (zVal * area - areaNeg)/ areaPos;
1310 z = unity34 - zNegMin * pow(zNegMax / zNegMin, zValMod);
1311 }
1312
1313 // 1 / (unity34 + z).
1314 } else if (iZ == 2) {
1315 double areaNeg = log(zNegMin / zNegMax);
1316 double areaPos = log(zPosMax / zPosMin);
1317 double area = areaNeg + areaPos;
1318 if (zVal * area < areaNeg) {
1319 double zValMod = zVal * area / areaNeg;
1320 z = zNegMax * pow(zNegMin / zNegMax, zValMod) - unity34;
1321 } else {
1322 double zValMod = (zVal * area - areaNeg)/ areaPos;
1323 z = zPosMin * pow(zPosMax / zPosMin, zValMod) - unity34;
1324 }
1325
1326 // 1 / (unity34 - z)^2.
1327 } else if (iZ == 3) {
1328 double areaNeg = 1. / zPosMin - 1. / zPosMax;
1329 double areaPos = 1. / zNegMax - 1. / zNegMin;
1330 double area = areaNeg + areaPos;
1331 if (zVal * area < areaNeg) {
1332 double zValMod = zVal * area / areaNeg;
1333 z = unity34 - 1. / (1./zPosMax + areaNeg * zValMod);
1334 } else {
1335 double zValMod = (zVal * area - areaNeg)/ areaPos;
1336 z = unity34 - 1. / (1./zNegMin + areaPos * zValMod);
1337 }
1338
1339 // 1 / (unity34 + z)^2.
1340 } else if (iZ == 4) {
1341 double areaNeg = 1. / zNegMax - 1. / zNegMin;
1342 double areaPos = 1. / zPosMin - 1. / zPosMax;
1343 double area = areaNeg + areaPos;
1344 if (zVal * area < areaNeg) {
1345 double zValMod = zVal * area / areaNeg;
1346 z = 1. / (1./zNegMax - areaNeg * zValMod) - unity34;
1347 } else {
1348 double zValMod = (zVal * area - areaNeg)/ areaPos;
1349 z = 1. / (1./zPosMin - areaPos * zValMod) - unity34;
1350 }
1351 }
1352
1353 // Safety check for roundoff errors. Combinations with z.
1354 if (z < 0.) z = min(-zMin, max(-zMax, z));
1355 else z = min(zMax, max(zMin, z));
1356 zNeg = max(ratio34, unity34 - z);
1357 zPos = max(ratio34, unity34 + z);
1358
1359 // Phase space integral in z.
1360 double intZ0 = 2. * (zMax - zMin);
1361 double intZ12 = log( (zPosMax * zNegMin) / (zPosMin * zNegMax) );
1362 double intZ34 = 1. / zPosMin - 1. / zPosMax + 1. / zNegMax
1363 - 1. / zNegMin;
1364 wtZ = mHat * pAbs / ( (zCoef[0] / intZ0) + (zCoef[1] / intZ12) / zNeg
1365 + (zCoef[2] / intZ12) / zPos + (zCoef[3] / intZ34) / pow2(zNeg)
1366 + (zCoef[4] / intZ34) / pow2(zPos) );
1367
1368 // Calculate tHat and uHat. Also gives pTHat.
1369 double sH34 = -0.5 * (sH - s3 - s4);
1370 tH = sH34 + mHat * pAbs * z;
1371 uH = sH34 - mHat * pAbs * z;
1372 pTH = sqrtpos( (tH * uH - s3 * s4) / sH);
1373
1374}
1375
1376//*********
1377
1378// Select three-body phase space according to a cylindrically based form
1379// that can be chosen to favour low pT based on the form of propagators.
1380
1381bool PhaseSpace::select3Body() {
1382
1383 // Upper and lower limits of pT choice for 4 and 5.
1384 double m35S = pow2(m3 + m5);
1385 double pT4Smax = 0.25 * ( pow2(sH - s4 - m35S) - 4. * s4 * m35S ) / sH;
1386 if (pTHatMax > pTHatMin) pT4Smax = min( pT2HatMax, pT4Smax);
1387 double pT4Smin = pT2HatMin;
1388 double m34S = pow2(m3 + m4);
1389 double pT5Smax = 0.25 * ( pow2(sH - s5 - m34S) - 4. * s5 * m34S ) / sH;
1390 if (pTHatMax > pTHatMin) pT5Smax = min( pT2HatMax, pT5Smax);
1391 double pT5Smin = pT2HatMin;
1392
1393 // Check that pT ranges not closed.
1394 if ( pT4Smax < pow2(pTHatMin + MASSMARGIN) ) return false;
1395 if ( pT5Smax < pow2(pTHatMin + MASSMARGIN) ) return false;
1396
1397 // Select pT4S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2.
1398 double pTSmaxProp = pT4Smax + sTchan1;
1399 double pTSminProp = pT4Smin + sTchan1;
1400 double pTSratProp = pTSmaxProp / pTSminProp;
1401 double pTSdiff = pT4Smax - pT4Smin;
1402 double rShape = Rndm::flat();
1403 double pT4S = 0.;
1404 if (rShape < frac3Flat) pT4S = pT4Smin + Rndm::flat() * pTSdiff;
1405 else if (rShape < frac3Flat + frac3Pow1) pT4S = max( pT2HatMin,
1406 pTSminProp * pow( pTSratProp, Rndm::flat() ) - sTchan1 );
1407 else pT4S = max( pT2HatMin, pTSminProp * pTSmaxProp
1408 / (pTSminProp + Rndm::flat()* pTSdiff) - sTchan1 );
1409 double wt4 = pTSdiff / ( frac3Flat
1410 + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT4S + sTchan1))
1411 + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT4S + sTchan1) );
1412
1413 // Select pT5S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2.
1414 pTSmaxProp = pT5Smax + sTchan2;
1415 pTSminProp = pT5Smin + sTchan2;
1416 pTSratProp = pTSmaxProp / pTSminProp;
1417 pTSdiff = pT5Smax - pT5Smin;
1418 rShape = Rndm::flat();
1419 double pT5S = 0.;
1420 if (rShape < frac3Flat) pT5S = pT5Smin + Rndm::flat() * pTSdiff;
1421 else if (rShape < frac3Flat + frac3Pow1) pT5S = max( pT2HatMin,
1422 pTSminProp * pow( pTSratProp, Rndm::flat() ) - sTchan2 );
1423 else pT5S = max( pT2HatMin, pTSminProp * pTSmaxProp
1424 / (pTSminProp + Rndm::flat()* pTSdiff) - sTchan2 );
1425 double wt5 = pTSdiff / ( frac3Flat
1426 + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT5S + sTchan2))
1427 + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT5S + sTchan2) );
1428
1429 // Select azimuthal angles and check that third pT in range.
1430 double phi4 = 2. * M_PI * Rndm::flat();
1431 double phi5 = 2. * M_PI * Rndm::flat();
1432 double pT3S = max( 0., pT4S + pT5S + 2. * sqrt(pT4S * pT5S)
1433 * cos(phi4 - phi5) );
1434 if ( pT3S < pT2HatMin || (pTHatMax > pTHatMin && pT3S > pT2HatMax) )
1435 return false;
1436
1437 // Calculate transverse masses and check that phase space not closed.
1438 double sT3 = s3 + pT3S;
1439 double sT4 = s4 + pT4S;
1440 double sT5 = s5 + pT5S;
1441 double mT3 = sqrt(sT3);
1442 double mT4 = sqrt(sT4);
1443 double mT5 = sqrt(sT5);
1444 if ( mT3 + mT4 + mT5 + MASSMARGIN > mHat ) return false;
1445
1446 // Select rapidity for particle 3 and check that phase space not closed.
1447 double m45S = pow2(mT4 + mT5);
1448 double y3max = log( ( sH + sT3 - m45S + sqrtpos( pow2(sH - sT3 - m45S)
1449 - 4 * sT3 * m45S ) ) / (2. * mHat * mT3) );
1450 if (y3max < YRANGEMARGIN) return false;
1451 double y3 = (2. * Rndm::flat() - 1.) * (1. - YRANGEMARGIN) * y3max;
1452 double pz3 = mT3 * sinh(y3);
1453 double e3 = mT3 * cosh(y3);
1454
1455 // Find momentum transfers in the two mirror solutions (in 4-5 frame).
1456 double pz45 = -pz3;
1457 double e45 = mHat - e3;
1458 double sT45 = e45 * e45 - pz45 * pz45;
1459 double lam45 = sqrtpos( pow2(sT45 - sT4 - sT5) - 4. * sT4 * sT5 );
1460 if (lam45 < YRANGEMARGIN * sH) return false;
1461 double lam4e = sT45 + sT4 - sT5;
1462 double lam5e = sT45 + sT5 - sT4;
1463 double tFac = -0.5 * mHat / sT45;
1464 double t1Pos = tFac * (e45 - pz45) * (lam4e - lam45);
1465 double t1Neg = tFac * (e45 - pz45) * (lam4e + lam45);
1466 double t2Pos = tFac * (e45 + pz45) * (lam5e - lam45);
1467 double t2Neg = tFac * (e45 + pz45) * (lam5e + lam45);
1468
1469 // Construct relative mirror weights and make choice.
1470 double wtPosUnnorm = 1.;
1471 double wtNegUnnorm = 1.;
1472 if (useMirrorWeight) {
1473 wtPosUnnorm = 1./ pow2( (t1Pos - sTchan1) * (t2Pos - sTchan2) );
1474 wtNegUnnorm = 1./ pow2( (t1Neg - sTchan1) * (t2Neg - sTchan2) );
1475 }
1476 double wtPos = wtPosUnnorm / (wtPosUnnorm + wtNegUnnorm);
1477 double wtNeg = wtNegUnnorm / (wtPosUnnorm + wtNegUnnorm);
1478 double epsilon = (Rndm::flat() < wtPos) ? 1. : -1.;
1479
1480 // Construct four-vectors in rest frame of subprocess.
1481 double px4 = sqrt(pT4S) * cos(phi4);
1482 double py4 = sqrt(pT4S) * sin(phi4);
1483 double px5 = sqrt(pT5S) * cos(phi5);
1484 double py5 = sqrt(pT5S) * sin(phi5);
1485 double pz4 = 0.5 * (pz45 * lam4e + epsilon * e45 * lam45) / sT45;
1486 double pz5 = pz45 - pz4;
1487 double e4 = sqrt(sT4 + pz4 * pz4);
1488 double e5 = sqrt(sT5 + pz5 * pz5);
1489 p3cm = Vec4( -(px4 + px5), -(py4 + py5), pz3, e3);
1490 p4cm = Vec4( px4, py4, pz4, e4);
1491 p5cm = Vec4( px5, py5, pz5, e5);
1492
1493 // Total weight to associate with kinematics choice.
1494 wt3Body = wt4 * wt5 * (2. * y3max) / (128. * pow3(M_PI) * lam45);
1495 wt3Body *= (epsilon > 0.) ? 1. / wtPos : 1. / wtNeg;
1496
1497 // Cross section of form |M|^2/(2 sHat) dPS_3 so need 1/(2 sHat).
1498 wt3Body /= (2. * sH);
1499
1500 // Done.
1501 return true;
1502
1503}
1504
1505//*********
1506
1507// Solve linear equation system for better phase space coefficients.
1508
1509void PhaseSpace::solveSys( int n, int bin[8], double vec[8],
1510 double mat[8][8], double coef[8], ostream& os) {
1511
1512 // Optional printout.
1513 if (showSearch) {
1514 os << "\n Equation system: " << setw(5) << bin[0];
1515 for (int j = 0; j < n; ++j) os << setw(12) << mat[0][j];
1516 os << setw(12) << vec[0] << "\n";
1517 for (int i = 1; i < n; ++i) {
1518 os << " " << setw(5) << bin[i];
1519 for (int j = 0; j < n; ++j) os << setw(12) << mat[i][j];
1520 os << setw(12) << vec[i] << "\n";
1521 }
1522 }
1523
1524 // Local variables.
1525 double vecNor[8], coefTmp[8];
1526 for (int i = 0; i < n; ++i) coefTmp[i] = 0.;
1527
1528 // Check if equation system solvable.
1529 bool canSolve = true;
1530 for (int i = 0; i < n; ++i) if (bin[i] == 0) canSolve = false;
1531 double vecSum = 0.;
1532 for (int i = 0; i < n; ++i) vecSum += vec[i];
1533 if (abs(vecSum) < TINY) canSolve = false;
1534
1535 // Solve to find relative importance of cross-section pieces.
1536 if (canSolve) {
1537 for (int i = 0; i < n; ++i) vecNor[i] = max( 0.1, vec[i] / vecSum);
1538 for (int k = 0; k < n - 1; ++k) {
1539 for (int i = k + 1; i < n; ++i) {
1540 if (abs(mat[k][k]) < TINY) {canSolve = false; break;}
1541 double ratio = mat[i][k] / mat[k][k];
1542 vec[i] -= ratio * vec[k];
1543 for (int j = k; j < n; ++j) mat[i][j] -= ratio * mat[k][j];
1544 }
1545 if (!canSolve) break;
1546 }
1547 if (canSolve) {
1548 for (int k = n - 1; k >= 0; --k) {
1549 for (int j = k + 1; j < n; ++j) vec[k] -= mat[k][j] * coefTmp[j];
1550 coefTmp[k] = vec[k] / mat[k][k];
1551 }
1552 }
1553 }
1554
1555 // Share evenly if failure.
1556 if (!canSolve) for (int i = 0; i < n; ++i) {
1557 coefTmp[i] = 1.;
1558 vecNor[i] = 0.1;
1559 if (vecSum > TINY) vecNor[i] = max(0.1, vec[i] / vecSum);
1560 }
1561
1562 // Normalize coefficients, with piece shared democratically.
1563 double coefSum = 0.;
1564 vecSum = 0.;
1565 for (int i = 0; i < n; ++i) {
1566 coefTmp[i] = max( 0., coefTmp[i]);
1567 coefSum += coefTmp[i];
1568 vecSum += vecNor[i];
1569 }
1570 if (coefSum > 0.) for (int i = 0; i < n; ++i) coef[i] = EVENFRAC / n
1571 + (1. - EVENFRAC) * 0.5 * (coefTmp[i] / coefSum + vecNor[i] / vecSum);
1572 else for (int i = 0; i < n; ++i) coef[i] = 1. / n;
1573
1574 // Optional printout.
1575 if (showSearch) {
1576 os << " Solution: ";
1577 for (int i = 0; i < n; ++i) os << setw(12) << coef[i];
1578 os << "\n";
1579 }
1580}
1581
1582//*********
1583
1584// Setup mass selection for one resonance at a time - part 1.
1585
1586void PhaseSpace::setupMass1(int iM) {
1587
1588 // Identity for mass seletion; is 0 also for light quarks (not yet selected).
1589 if (iM == 3) idMass[iM] = abs(sigmaProcessPtr->id3Mass());
1590 if (iM == 4) idMass[iM] = abs(sigmaProcessPtr->id4Mass());
1591 if (iM == 5) idMass[iM] = abs(sigmaProcessPtr->id5Mass());
1592
1593 // Masses and widths of resonances.
1594 if (idMass[iM] == 0) {
1595 mPeak[iM] = 0.;
1596 mWidth[iM] = 0.;
1597 mMin[iM] = 0.;
1598 mMax[iM] = 0.;
1599 } else {
1600 mPeak[iM] = ParticleDataTable::m0(idMass[iM]);
1601 mWidth[iM] = ParticleDataTable::mWidth(idMass[iM]);
1602 mMin[iM] = ParticleDataTable::mMin(idMass[iM]);
1603 mMax[iM] = ParticleDataTable::mMax(idMass[iM]);
1604 // gmZmode == 1 means pure photon propagator; set at lower mass limit.
1605 if (idMass[iM] == 23 && gmZmode == 1) mPeak[iM] = mMin[iM];
1606 }
1607
1608 // Mass and width combinations for Breit-Wigners.
1609 sPeak[iM] = mPeak[iM] * mPeak[iM];
1610 useBW[iM] = useBreitWigners && (mWidth[iM] > minWidthBreitWigners);
1611 if (!useBW[iM]) mWidth[iM] = 0.;
1612 mw[iM] = mPeak[iM] * mWidth[iM];
1613 wmRat[iM] = (idMass[iM] == 0 || mPeak[iM] == 0.)
1614 ? 0. : mWidth[iM] / mPeak[iM];
1615
1616 // Simple Breit-Wigner range, upper edge to be corrected subsequently.
1617 if (useBW[iM]) {
1618 mLower[iM] = mMin[iM];
1619 mUpper[iM] = mHatMax;
1620 }
1621
1622}
1623
1624//*********
1625
1626// Setup mass selection for one resonance at a time - part 2.
1627
1628void PhaseSpace::setupMass2(int iM, double distToThresh) {
1629
1630 // Store reduced Breit-Wigner range.
1631 if (mMax[iM] > mMin[iM]) mUpper[iM] = min( mUpper[iM], mMax[iM]);
1632 sLower[iM] = mLower[iM] * mLower[iM];
1633 sUpper[iM] = mUpper[iM] * mUpper[iM];
1634
1635 // Prepare to select m3 by BW + flat + 1/s_3.
1636 // Determine relative coefficients by allowed mass range.
1637 if (distToThresh > THRESHOLDSIZE) {
1638 fracFlat[iM] = 0.1;
1639 fracInv[iM] = 0.1;
1640 } else if (distToThresh > - THRESHOLDSIZE) {
1641 fracFlat[iM] = 0.25 - 0.15 * distToThresh / THRESHOLDSIZE;
1642 fracInv [iM] = 0.15 - 0.05 * distToThresh / THRESHOLDSIZE;
1643 } else {
1644 fracFlat[iM] = 0.4;
1645 fracInv[iM] = 0.2;
1646 }
1647
1648 // For gamma*/Z0: increase 1/s_i part and introduce 1/s_i^2 part.
1649 fracInv2[iM] = 0.;
1650 if (idMass[iM] == 23 && gmZmode == 0) {
1651 fracFlat[iM] *= 0.5;
1652 fracInv[iM] = 0.5 * fracInv[iM] + 0.25;
1653 fracInv2[iM] = 0.25;
1654 } else if (idMass[iM] == 23 && gmZmode == 1) {
1655 fracFlat[iM] = 0.1;
1656 fracInv[iM] = 0.4;
1657 fracInv2[iM] = 0.4;
1658 }
1659
1660 // Normalization integrals for the respective contribution.
1661 atanLower[iM] = atan( (sLower[iM] - sPeak[iM])/ mw[iM] );
1662 atanUpper[iM] = atan( (sUpper[iM] - sPeak[iM])/ mw[iM] );
1663 intBW[iM] = atanUpper[iM] - atanLower[iM];
1664 intFlat[iM] = sUpper[iM] - sLower[iM];
1665 intInv[iM] = log( sUpper[iM] / sLower[iM] );
1666 intInv2[iM] = 1./sLower[iM] - 1./sUpper[iM];
1667
1668}
1669
1670//*********
1671
1672// Select Breit-Wigner-distributed or fixed masses.
1673
1674void PhaseSpace::trialMass(int iM) {
1675
1676 // References to masses to be set.
1677 double& mSet = (iM == 3) ? m3 : ( (iM == 4) ? m4 : m5 );
1678 double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 );
1679
1680 // Distribution for m_i is BW + flat + 1/s_i + 1/s_i^2.
1681 if (useBW[iM]) {
1682 double pickForm = Rndm::flat();
1683 if (pickForm > fracFlat[iM] + fracInv[iM] + fracInv2[iM])
1684 sSet = sPeak[iM] + mw[iM] * tan( atanLower[iM]
1685 + Rndm::flat() * intBW[iM] );
1686 else if (pickForm > fracInv[iM] + fracInv2[iM])
1687 sSet = sLower[iM] + Rndm::flat() * (sUpper[iM] - sLower[iM]);
1688 else if (pickForm > fracInv2[iM])
1689 sSet = sLower[iM] * pow( sUpper[iM] / sLower[iM], Rndm::flat() );
1690 else sSet = sLower[iM] * sUpper[iM]
1691 / (sLower[iM] + Rndm::flat() * (sUpper[iM] - sLower[iM]));
1692 mSet = sqrt(sSet);
1693
1694 // Else m_i is fixed at peak value.
1695 } else {
1696 mSet = mPeak[iM];
1697 sSet = sPeak[iM];
1698 }
1699
1700}
1701
1702//*********
1703
1704// Naively a fixed-width Breit-Wigner is used to pick the mass.
1705// Here come the correction factors for
1706// (i) preselection according to BW + flat in s_i + 1/s_i + 1/s_i^2,
1707// (ii) reduced allowed mass range,
1708// (iii) running width, i.e. m0*Gamma0 -> s*Gamma0/m0.
1709// In the end, the weighted distribution is a running-width BW.
1710
1711double PhaseSpace::weightMass(int iM) {
1712
1713 // Reference to mass and to Breit-Wigner weight to be set.
1714 double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 );
1715 double& runBWH = (iM == 3) ? runBW3H : ( (iM == 4) ? runBW4H : runBW5H );
1716
1717 // Default weight if no Breit-Wigner.
1718 runBWH = 1.;
1719 if (!useBW[iM]) return 1.;
1720
1721 // Weight of generated distribution.
1722 double genBW = (1. - fracFlat[iM] - fracInv[iM] - fracInv2[iM])
1723 * mw[iM] / ( (pow2(sSet - sPeak[iM]) + pow2(mw[iM])) * intBW[iM])
1724 + fracFlat[iM] / intFlat[iM] + fracInv[iM] / (sSet * intInv[iM])
1725 + fracInv2[iM] / (sSet*sSet * intInv2[iM]);
1726
1727 // Weight of distribution with running width in Breit-Wigner.
1728 double mwRun = sSet * wmRat[iM];
1729 runBWH = mwRun / (pow2(sSet - sPeak[iM]) + pow2(mwRun)) / M_PI;
1730
1731 // Done.
1732 return (runBWH / genBW);
1733
1734}
1735
1736//**************************************************************************
1737
1738// PhaseSpace2to1tauy class.
1739// 2 -> 1 kinematics for normal subprocesses.
1740
1741//*********
1742
1743// Set limits for resonance mass selection.
1744
1745bool PhaseSpace2to1tauy::setupMass() {
1746
1747 // Treat Z0 as such or as gamma*/Z0
1748 gmZmode = gmZmodeGlobal;
1749 int gmZmodeProc = sigmaProcessPtr->gmZmode();
1750 if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
1751
1752 // Mass limits for current resonance.
1753 int idRes = abs(sigmaProcessPtr->resonanceA());
1754 int idTmp = abs(sigmaProcessPtr->resonanceB());
1755 if (idTmp > 0) idRes = idTmp;
1756 double mResMin = (idRes == 0) ? 0. : ParticleDataTable::mMin(idRes);
1757 double mResMax = (idRes == 0) ? 0. : ParticleDataTable::mMax(idRes);
1758
1759 // Compare with global mass limits and pick tighter of them.
1760 mHatMin = max( mResMin, mHatGlobalMin);
1761 sHatMin = mHatMin*mHatMin;
1762 mHatMax = eCM;
1763 if (mResMax > mResMin) mHatMax = min( mHatMax, mResMax);
1764 if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( mHatMax, mHatGlobalMax);
1765 sHatMax = mHatMax*mHatMax;
1766
1767 // Default Breit-Wigner weight.
1768 wtBW = 1.;
1769
1770 // Fail if mass window (almost) closed.
1771 return (mHatMax > mHatMin + MASSMARGIN);
1772
1773}
1774
1775//*********
1776
1777// Construct the four-vector kinematics from the trial values.
1778
1779bool PhaseSpace2to1tauy::finalKin() {
1780
1781 // Particle masses; incoming always on mass shell.
1782 mH[1] = 0.;
1783 mH[2] = 0.;
1784 mH[3] = mHat;
1785
1786 // Incoming partons along beam axes. Outgoing has sum of momenta.
1787 pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
1788 pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
1789 pH[3] = pH[1] + pH[2];
1790
1791 // Done.
1792 return true;
1793}
1794
1795//**************************************************************************
1796
1797// PhaseSpace2to2tauyz class.
1798// 2 -> 2 kinematics for normal subprocesses.
1799
1800//*********
1801
1802// Set up for fixed or Breit-Wigner mass selection.
1803
1804bool PhaseSpace2to2tauyz::setupMasses() {
1805
1806 // Treat Z0 as such or as gamma*/Z0
1807 gmZmode = gmZmodeGlobal;
1808 int gmZmodeProc = sigmaProcessPtr->gmZmode();
1809 if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
1810
1811 // Set sHat limits - based on global limits only.
1812 mHatMin = mHatGlobalMin;
1813 sHatMin = mHatMin*mHatMin;
1814 mHatMax = eCM;
1815 if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax);
1816 sHatMax = mHatMax*mHatMax;
1817
1818 // Masses and widths of resonances.
1819 setupMass1(3);
1820 setupMass1(4);
1821
1822 // Reduced mass range when two massive particles.
1823 if (useBW[3]) mUpper[3] -= (useBW[4]) ? mMin[4] : mPeak[4];
1824 if (useBW[4]) mUpper[4] -= (useBW[3]) ? mMin[3] : mPeak[3];
1825
1826 // If closed phase space then unallowed process.
1827 bool physical = true;
1828 if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false;
1829 if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false;
1830 if (!useBW[3] && !useBW[4] && mHatMax < mPeak[3] + mPeak[4] + MASSMARGIN)
1831 physical = false;
1832 if (!physical) return false;
1833
1834 // If either particle is massless then need extra pTHat cut.
1835 pTHatMin = pTHatGlobalMin;
1836 if (mPeak[3] < pTHatMinDiverge || mPeak[4] < pTHatMinDiverge)
1837 pTHatMin = max( pTHatMin, pTHatMinDiverge);
1838 pT2HatMin = pTHatMin * pTHatMin;
1839 pTHatMax = pTHatGlobalMax;
1840 pT2HatMax = pTHatMax * pTHatMax;
1841
1842 // Prepare to select m3 by BW + flat + 1/s_3.
1843 if (useBW[3]) {
1844 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[3]
1845 / (pow2(mWidth[3]) + pow2(mWidth[4]));
1846 double distToThreshB = (mHatMax - mPeak[3] - mMin[4]) / mWidth[3];
1847 double distToThresh = min( distToThreshA, distToThreshB);
1848 setupMass2(3, distToThresh);
1849 }
1850
1851 // Prepare to select m4 by BW + flat + 1/s_4.
1852 if (useBW[4]) {
1853 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[4]
1854 / (pow2(mWidth[3]) + pow2(mWidth[4]));
1855 double distToThreshB = (mHatMax - mMin[3] - mPeak[4]) / mWidth[4];
1856 double distToThresh = min( distToThreshA, distToThreshB);
1857 setupMass2(4, distToThresh);
1858 }
1859
1860 // Initialization masses. Special cases when constrained phase space.
1861 m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3];
1862 m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4];
1863 if (m3 + m4 + THRESHOLDSIZE * (mWidth[3] + mWidth[4]) + MASSMARGIN
1864 > mHatMax) {
1865 if (useBW[3] && useBW[4]) physical = constrainedM3M4();
1866 else if (useBW[3]) physical = constrainedM3();
1867 else if (useBW[4]) physical = constrainedM4();
1868 }
1869 s3 = m3*m3;
1870 s4 = m4*m4;
1871
1872 // Correct selected mass-spectrum to running-width Breit-Wigner.
1873 // Extra safety margin for maximum search.
1874 wtBW = 1.;
1875 if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX;
1876 if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX;
1877
1878 // Done.
1879 return physical;
1880
1881}
1882
1883
1884//*********
1885
1886// Select Breit-Wigner-distributed or fixed masses.
1887
1888bool PhaseSpace2to2tauyz::trialMasses() {
1889
1890 // By default vanishing cross section.
1891 sigmaNw = 0.;
1892 wtBW = 1.;
1893
1894 // Pick m3 and m4 independently.
1895 trialMass(3);
1896 trialMass(4);
1897
1898 // If outside phase space then reject event.
1899 if (m3 + m4 + MASSMARGIN > mHatMax) return false;
1900
1901 // Correct selected mass-spectrum to running-width Breit-Wigner.
1902 if (useBW[3]) wtBW *= weightMass(3);
1903 if (useBW[4]) wtBW *= weightMass(4);
1904
1905 // Done.
1906 return true;
1907}
1908
1909//*********
1910
1911// Construct the four-vector kinematics from the trial values.
1912
1913bool PhaseSpace2to2tauyz::finalKin() {
1914
1915 // Assign masses to particles assumed massless in matrix elements.
1916 int id3 = sigmaProcessPtr->id(3);
1917 int id4 = sigmaProcessPtr->id(4);
1918 if (idMass[3] == 0) { m3 = ParticleDataTable::m0(id3); s3 = m3*m3; }
1919 if (idMass[4] == 0) { m4 = ParticleDataTable::m0(id4); s4 = m4*m4; }
1920
1921 // Sometimes swap tHat <-> uHat to reflect chosen final-state order.
1922 if (sigmaProcessPtr->swappedTU()) {
1923 swap(tH, uH);
1924 z = -z;
1925 }
1926
1927 // Check that phase space still open after new mass assignment.
1928 if (m3 + m4 + MASSMARGIN > mHat) {
1929 infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::finalKin: "
1930 "failed after mass assignment");
1931 return false;
1932 }
1933 p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1934 pAbs = sqrtpos( p2Abs );
1935
1936 // Particle masses; incoming always on mass shell.
1937 mH[1] = 0.;
1938 mH[2] = 0.;
1939 mH[3] = m3;
1940 mH[4] = m4;
1941
1942 // Incoming partons along beam axes.
1943 pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
1944 pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
1945
1946 // Outgoing partons initially in collision CM frame along beam axes.
1947 pH[3] = Vec4( 0., 0., pAbs, 0.5 * (sH + s3 - s4) / mHat);
1948 pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (sH + s4 - s3) / mHat);
1949
1950 // Then rotate and boost them to overall CM frame
1951 theta = acos(z);
1952 phi = 2. * M_PI * Rndm::flat();
1953 betaZ = (x1H - x2H)/(x1H + x2H);
1954 pH[3].rot( theta, phi);
1955 pH[4].rot( theta, phi);
1956 pH[3].bst( 0., 0., betaZ);
1957 pH[4].bst( 0., 0., betaZ);
1958 pTH = pAbs * sin(theta);
1959
1960 // Done.
1961 return true;
1962}
1963
1964//*********
1965
1966// Special choice of m3 and m4 when mHatMax push them off mass shell.
1967// Vary x in expression m3 + m4 = mHatMax - x * (Gamma3 + Gamma4).
1968// For each x try to put either 3 or 4 as close to mass shell as possible.
1969// Maximize BW_3 * BW_4 * beta_34, where latter approximate phase space.
1970
1971bool PhaseSpace2to2tauyz::constrainedM3M4() {
1972
1973 // Initial values.
1974 bool foundNonZero = false;
1975 double wtMassMax = 0.;
1976 double m3WtMax = 0.;
1977 double m4WtMax = 0.;
1978 double xMax = (mHatMax - mLower[3] - mLower[4]) / (mWidth[3] + mWidth[4]);
1979 double xStep = THRESHOLDSTEP * min(1., xMax);
1980 double xNow = 0.;
1981 double wtMassXbin, wtMassMaxOld, m34, mT34Min, wtMassNow,
1982 wtBW3Now, wtBW4Now, beta34Now;
1983
1984 // Step through increasing x values.
1985 do {
1986 xNow += xStep;
1987 wtMassXbin = 0.;
1988 wtMassMaxOld = wtMassMax;
1989 m34 = mHatMax - xNow * (mWidth[3] + mWidth[4]);
1990
1991 // Study point where m3 as close as possible to on-shell.
1992 m3 = min( mUpper[3], m34 - mLower[4]);
1993 if (m3 > mPeak[3]) m3 = max( mLower[3], mPeak[3]);
1994 m4 = m34 - m3;
1995 if (m4 < mLower[4]) {m4 = mLower[4]; m3 = m34 - m4;}
1996
1997 // Check that inside phase space limit set by pTmin.
1998 mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin);
1999 if (mT34Min < mHatMax) {
2000
2001 // Breit-Wigners and beta factor give total weight.
2002 wtMassNow = 0.;
2003 if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4]
2004 && m4 < mUpper[4]) {
2005 wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2006 wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2007 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2008 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2009 wtMassNow = wtBW3Now * wtBW4Now * beta34Now;
2010 }
2011
2012 // Store new maximum, if any.
2013 if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow;
2014 if (wtMassNow > wtMassMax) {
2015 foundNonZero = true;
2016 wtMassMax = wtMassNow;
2017 m3WtMax = m3;
2018 m4WtMax = m4;
2019 }
2020 }
2021
2022 // Study point where m4 as close as possible to on-shell.
2023 m4 = min( mUpper[4], m34 - mLower[3]);
2024 if (m4 > mPeak[4]) m4 = max( mLower[4], mPeak[4]);
2025 m3 = m34 - m4;
2026 if (m3 < mLower[3]) {m3 = mLower[3]; m4 = m34 - m3;}
2027
2028 // Check that inside phase space limit set by pTmin.
2029 mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin);
2030 if (mT34Min < mHatMax) {
2031
2032 // Breit-Wigners and beta factor give total weight.
2033 wtMassNow = 0.;
2034 if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4]
2035 && m4 < mUpper[4]) {
2036 wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2037 wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2038 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2039 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2040 wtMassNow = wtBW3Now * wtBW4Now * beta34Now;
2041 }
2042
2043 // Store new maximum, if any.
2044 if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow;
2045 if (wtMassNow > wtMassMax) {
2046 foundNonZero = true;
2047 wtMassMax = wtMassNow;
2048 m3WtMax = m3;
2049 m4WtMax = m4;
2050 }
2051 }
2052
2053 // Continue stepping if increasing trend and more x range available.
2054 } while ( (!foundNonZero || wtMassXbin > wtMassMaxOld)
2055 && xNow < xMax - xStep);
2056
2057 // Restore best values for subsequent maximization. Return.
2058 m3 = m3WtMax;
2059 m4 = m4WtMax;
2060 return foundNonZero;
2061
2062}
2063
2064//*********
2065
2066// Special choice of m3 when mHatMax pushes it off mass shell.
2067// Vary x in expression m3 = mHatMax - m4 - x * Gamma3.
2068// Maximize BW_3 * beta_34, where latter approximate phase space.
2069
2070bool PhaseSpace2to2tauyz::constrainedM3() {
2071
2072 // Initial values.
2073 bool foundNonZero = false;
2074 double wtMassMax = 0.;
2075 double m3WtMax = 0.;
2076 double mT4Min = sqrt(m4*m4 + pT2HatMin);
2077 double xMax = (mHatMax - mLower[3] - m4) / mWidth[3];
2078 double xStep = THRESHOLDSTEP * min(1., xMax);
2079 double xNow = 0.;
2080 double wtMassNow, mT34Min, wtBW3Now, beta34Now;
2081
2082 // Step through increasing x values; gives m3 unambiguously.
2083 do {
2084 xNow += xStep;
2085 wtMassNow = 0.;
2086 m3 = mHatMax - m4 - xNow * mWidth[3];
2087
2088 // Check that inside phase space limit set by pTmin.
2089 mT34Min = sqrt(m3*m3 + pT2HatMin) + mT4Min;
2090 if (mT34Min < mHatMax) {
2091
2092 // Breit-Wigner and beta factor give total weight.
2093 wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2094 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2095 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2096 wtMassNow = wtBW3Now * beta34Now;
2097
2098 // Store new maximum, if any.
2099 if (wtMassNow > wtMassMax) {
2100 foundNonZero = true;
2101 wtMassMax = wtMassNow;
2102 m3WtMax = m3;
2103 }
2104 }
2105
2106 // Continue stepping if increasing trend and more x range available.
2107 } while ( (!foundNonZero || wtMassNow > wtMassMax)
2108 && xNow < xMax - xStep);
2109
2110 // Restore best value for subsequent maximization. Return.
2111 m3 = m3WtMax;
2112 return foundNonZero;
2113
2114}
2115
2116//*********
2117
2118// Special choice of m4 when mHatMax pushes it off mass shell.
2119// Vary x in expression m4 = mHatMax - m3 - x * Gamma4.
2120// Maximize BW_4 * beta_34, where latter approximate phase space.
2121
2122bool PhaseSpace2to2tauyz::constrainedM4() {
2123
2124 // Initial values.
2125 bool foundNonZero = false;
2126 double wtMassMax = 0.;
2127 double m4WtMax = 0.;
2128 double mT3Min = sqrt(m3*m3 + pT2HatMin);
2129 double xMax = (mHatMax - mLower[4] - m3) / mWidth[4];
2130 double xStep = THRESHOLDSTEP * min(1., xMax);
2131 double xNow = 0.;
2132 double wtMassNow, mT34Min, wtBW4Now, beta34Now;
2133
2134 // Step through increasing x values; gives m4 unambiguously.
2135 do {
2136 xNow += xStep;
2137 wtMassNow = 0.;
2138 m4 = mHatMax - m3 - xNow * mWidth[4];
2139
2140 // Check that inside phase space limit set by pTmin.
2141 mT34Min = mT3Min + sqrt(m4*m4 + pT2HatMin);
2142 if (mT34Min < mHatMax) {
2143
2144 // Breit-Wigner and beta factor give total weight.
2145 wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2146 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2147 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2148 wtMassNow = wtBW4Now * beta34Now;
2149
2150 // Store new maximum, if any.
2151 if (wtMassNow > wtMassMax) {
2152 foundNonZero = true;
2153 wtMassMax = wtMassNow;
2154 m4WtMax = m4;
2155 }
2156 }
2157
2158 // Continue stepping if increasing trend and more x range available.
2159 } while ( (!foundNonZero || wtMassNow > wtMassMax)
2160 && xNow < xMax - xStep);
2161
2162 // Restore best value for subsequent maximization.
2163 m4 = m4WtMax;
2164 return foundNonZero;
2165
2166}
2167
2168//**************************************************************************
2169
2170// PhaseSpace2to2elastic class.
2171// 2 -> 2 kinematics set up for elastic scattering.
2172
2173//*********
2174
2175// Constants: could be changed here if desired, but normally should not.
2176// These are of technical nature, as described for each.
2177
2178// Maximum positive/negative argument for exponentiation.
2179const double PhaseSpace2to2elastic::EXPMAX = 50.;
2180
2181// Conversion coefficients = 1/(16pi) * (mb <-> GeV^2).
2182const double PhaseSpace2to2elastic::CONVERTEL = 0.0510925;
2183
2184//*********
2185
2186// Form of phase space sampling already fixed, so no optimization.
2187// However, need to read out relevant parameters from SigmaTotal.
2188
2189bool PhaseSpace2to2elastic::setupSampling() {
2190
2191 // Find maximum = value of cross section.
2192 sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2193 sigmaMx = sigmaNw;
2194
2195 // Squared masses of particles.
2196 s1 = mA * mA;
2197 s2 = mB * mB;
2198
2199 // Elastic slope.
2200 bSlope = sigmaTotPtr->bSlopeEl();
2201
2202 // Determine maximum possible t range.
2203 lambda12S = pow2(s - s1 - s2) - 4. * s1 * s2 ;
2204 tLow = - lambda12S / s;
2205 tUpp = 0;
2206
2207 // Production model with Coulomb corrections need more parameters.
2208 useCoulomb = Settings::flag("SigmaTotal:setOwn")
2209 && Settings::flag("SigmaElastic:setOwn");
2210 if (useCoulomb) {
2211 sigmaTot = sigmaTotPtr->sigmaTot();
2212 rho = Settings::parm("SigmaElastic:rho");
2213 lambda = Settings::parm("SigmaElastic:lambda");
2214 tAbsMin = Settings::parm("SigmaElastic:tAbsMin");
2215 phaseCst = Settings::parm("SigmaElastic:phaseConst");
2216 alphaEM0 = Settings::parm("StandardModel:alphaEM0");
2217
2218 // Relative rate of nuclear and Coulombic parts in trials.
2219 tUpp = -tAbsMin;
2220 sigmaNuc = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho) / bSlope
2221 * exp(-bSlope * tAbsMin);
2222 sigmaCou = (useCoulomb) ?
2223 pow2(alphaEM0) / (4. * CONVERTEL * tAbsMin) : 0.;
2224 signCou = (idA == idB) ? 1. : -1.;
2225
2226 // Dummy values.
2227 } else {
2228 sigmaNuc = sigmaNw;
2229 sigmaCou = 0.;
2230 }
2231
2232 // Calculate coefficient of generation.
2233 tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2234
2235 // Initialize alphaEM generation.
2236 alphaEM.init( Settings::mode("SigmaProcess:alphaEMorder") );
2237
2238 return true;
2239
2240}
2241
2242//*********
2243
2244// Select a trial kinematics phase space point. Perform full
2245// Monte Carlo acceptance/rejection at this stage.
2246
2247bool PhaseSpace2to2elastic::trialKin( bool, bool ) {
2248
2249 // Select t according to exp(bSlope*t).
2250 if (!useCoulomb || sigmaNuc > Rndm::flat() * (sigmaNuc + sigmaCou))
2251 tH = tUpp + log(1. + tAux * Rndm::flat()) / bSlope;
2252
2253 // Select t according to 1/t^2.
2254 else tH = tLow * tUpp / (tUpp + Rndm::flat() * (tLow - tUpp));
2255
2256 // Correction factor for ratio full/simulated.
2257 if (useCoulomb) {
2258 double sigmaN = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho)
2259 * exp(bSlope * tH);
2260 double alpEM = alphaEM.alphaEM(-tH);
2261 double sigmaC = pow2(alpEM) / (4. * CONVERTEL * tH*tH);
2262 double sigmaGen = 2. * (sigmaN + sigmaC);
2263 double form2 = pow4(lambda/(lambda - tH));
2264 double phase = signCou * alpEM
2265 * (-phaseCst - log(-0.5 * bSlope * tH));
2266 double sigmaCor = sigmaN + pow2(form2) * sigmaC
2267 - signCou * alpEM * sigmaTot * (form2 / (-tH))
2268 * exp(0.5 * bSlope * tH) * (rho * cos(phase) + sin(phase));
2269 sigmaNw = sigmaMx * sigmaCor / sigmaGen;
2270 }
2271
2272 // Careful reconstruction of scattering angle.
2273 double tRat = s * tH / lambda12S;
2274 double cosTheta = min(1., max(-1., 1. + 2. * tRat ) );
2275 double sinTheta = 2. * sqrtpos( -tRat * (1. + tRat) );
2276 theta = asin( min(1., sinTheta));
2277 if (cosTheta < 0.) theta = M_PI - theta;
2278
2279 return true;
2280
2281}
2282
2283//*********
2284
2285// Construct the four-vector kinematics from the trial values.
2286
2287bool PhaseSpace2to2elastic::finalKin() {
2288
2289 // Particle masses.
2290 mH[1] = mA;
2291 mH[2] = mB;
2292 mH[3] = m3;
2293 mH[4] = m4;
2294
2295 // Incoming particles along beam axes.
2296 pAbs = 0.5 * sqrtpos(lambda12S) / eCM;
2297 pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2298 pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2299
2300 // Outgoing particles initially along beam axes.
2301 pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2302 pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2303
2304 // Then rotate them
2305 phi = 2. * M_PI * Rndm::flat();
2306 pH[3].rot( theta, phi);
2307 pH[4].rot( theta, phi);
2308
2309 // Set some further info for completeness.
2310 x1H = 1.;
2311 x2H = 1.;
2312 sH = s;
2313 uH = 2. * (s1 + s2) - sH - tH;
2314 mHat = eCM;
2315 p2Abs = pAbs * pAbs;
2316 betaZ = 0.;
2317 pTH = pAbs * sin(theta);
2318
2319 // Done.
2320 return true;
2321
2322}
2323
2324//**************************************************************************
2325
2326// PhaseSpace2to2diffractive class.
2327// 2 -> 2 kinematics set up for diffractive scattering.
2328
2329//*********
2330
2331// Constants: could be changed here if desired, but normally should not.
2332// These are of technical nature, as described for each.
2333
2334// Number of tries to find acceptable (m^2, t) set.
2335const int PhaseSpace2to2diffractive::NTRY = 500;
2336
2337// Maximum positive/negative argument for exponentiation.
2338const double PhaseSpace2to2diffractive::EXPMAX = 50.;
2339
2340// Safety margin so sum of diffractive masses not too close to eCM.
2341const double PhaseSpace2to2diffractive::DIFFMASSMAX = 1e-8;
2342
2343//*********
2344
2345// Form of phase space sampling already fixed, so no optimization.
2346// However, need to read out relevant parameters from SigmaTotal.
2347
2348bool PhaseSpace2to2diffractive::setupSampling() {
2349
2350 // Find maximum = value of cross section.
2351 sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2352 sigmaMx = sigmaNw;
2353
2354 // Masses of particles and minimal masses of diffractive states.
2355 m3ElDiff = (isDiffA) ? sigmaTotPtr->mMinXB() : mA;
2356 m4ElDiff = (isDiffB) ? sigmaTotPtr->mMinAX() : mB;
2357 s1 = mA * mA;
2358 s2 = mB * mB;
2359 s3 = pow2( m3ElDiff);
2360 s4 = pow2( m4ElDiff);
2361
2362 // Parameters of low-mass-resonance diffractive enhancement.
2363 cRes = sigmaTotPtr->cRes();
2364 sResXB = pow2( sigmaTotPtr->mResXB());
2365 sResAX = pow2( sigmaTotPtr->mResAX());
2366 sProton = sigmaTotPtr->sProton();
2367
2368 // Lower limit diffractive slope.
2369 if (!isDiffB) bMin = sigmaTotPtr->bMinSlopeXB();
2370 else if (!isDiffA) bMin = sigmaTotPtr->bMinSlopeAX();
2371 else bMin = sigmaTotPtr->bMinSlopeXX();
2372
2373 // Determine maximum possible t range and coefficient of generation.
2374 lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
2375 lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2376 double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2377 double tempB = lambda12 * lambda34 / s;
2378 double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2379 * (s1 * s4 - s2 * s3) / s;
2380 tLow = -0.5 * (tempA + tempB);
2381 tUpp = tempC / tLow;
2382 tAux = exp( max(-EXPMAX, bMin * (tLow - tUpp)) ) - 1.;
2383
2384 return true;
2385
2386}
2387
2388//*********
2389
2390// Select a trial kinematics phase space point. Perform full
2391// Monte Carlo acceptance/rejection at this stage.
2392
2393bool PhaseSpace2to2diffractive::trialKin( bool, bool ) {
2394
2395 // Loop over attempts to set up masses and t consistently.
2396 for (int loop = 0; ; ++loop) {
2397 if (loop == NTRY) {
2398 infoPtr->errorMsg("Error in PhaseSpace2to2diffractive::trialKin: "
2399 " quit after repeated tries");
2400 return false;
2401 }
2402
2403 // Select diffractive mass/masses according to dm^2/m^2.
2404 m3 = (isDiffA) ? m3ElDiff * pow( max(mA, eCM - m4ElDiff) / m3ElDiff,
2405 Rndm::flat()) : m3ElDiff;
2406 m4 = (isDiffB) ? m4ElDiff * pow( max(mB, eCM - m3ElDiff) / m4ElDiff,
2407 Rndm::flat()) : m4ElDiff;
2408 s3 = m3 * m3;
2409 s4 = m4 * m4;
2410
2411 // Additional mass factors, including resonance enhancement.
2412 if (m3 + m4 >= eCM) continue;
2413 if (isDiffA && !isDiffB) {
2414 double facXB = (1. - s3 / s)
2415 * (1. + cRes * sResXB / (sResXB + s3));
2416 if (facXB < Rndm::flat() * (1. + cRes)) continue;
2417 } else if (isDiffB && !isDiffA) {
2418 double facAX = (1. - s4 / s)
2419 * (1. + cRes * sResAX / (sResAX + s4));
2420 if (facAX < Rndm::flat() * (1. + cRes)) continue;
2421 } else {
2422 double facXX = (1. - pow2(m3 + m4) / s)
2423 * (s * sProton / (s * sProton + s3 * s4))
2424 * (1. + cRes * sResXB / (sResXB + s3))
2425 * (1. + cRes * sResAX / (sResAX + s4));
2426 if (facXX < Rndm::flat() * pow2(1. + cRes)) continue;
2427 }
2428
2429 // Select t according to exp(bMin*t) and correct to right slope.
2430 tH = tUpp + log(1. + tAux * Rndm::flat()) / bMin;
2431 double bDiff = 0.;
2432 if (isDiffA && !isDiffB) bDiff = sigmaTotPtr->bSlopeXB(s3) - bMin;
2433 else if (!isDiffA) bDiff = sigmaTotPtr->bSlopeAX(s4) - bMin;
2434 else bDiff = sigmaTotPtr->bSlopeXX(s3, s4) - bMin;
2435 bDiff = max(0., bDiff);
2436 if (exp( max(-50., bDiff * (tH - tUpp)) ) < Rndm::flat()) continue;
2437
2438 // Check whether m^2 and t choices are consistent.
2439 lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2440 double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2441 double tempB = lambda12 * lambda34 / s;
2442 if (tempB < DIFFMASSMAX) continue;
2443 double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2444 * (s1 * s4 - s2 * s3) / s;
2445 double tLowNow = -0.5 * (tempA + tempB);
2446 double tUppNow = tempC / tLowNow;
2447 if (tH < tLowNow || tH > tUppNow) continue;
2448
2449 // Careful reconstruction of scattering angle.
2450 double cosTheta = min(1., max(-1., (tempA + 2. * tH) / tempB));
2451 double sinTheta = 2. * sqrtpos( -(tempC + tempA * tH + tH * tH) )
2452 / tempB;
2453 theta = asin( min(1., sinTheta));
2454 if (cosTheta < 0.) theta = M_PI - theta;
2455
2456 // Found acceptable kinematics, so no more looping.
2457 break;
2458 }
2459
2460 return true;
2461
2462}
2463
2464//*********
2465
2466// Construct the four-vector kinematics from the trial values.
2467
2468bool PhaseSpace2to2diffractive::finalKin() {
2469
2470 // Particle masses; incoming always on mass shell.
2471 mH[1] = mA;
2472 mH[2] = mB;
2473 mH[3] = m3;
2474 mH[4] = m4;
2475
2476 // Incoming particles along beam axes.
2477 pAbs = 0.5 * lambda12 / eCM;
2478 pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2479 pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2480
2481 // Outgoing particles initially along beam axes.
2482 pAbs = 0.5 * lambda34 / eCM;
2483 pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s3 - s4) / eCM);
2484 pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s4 - s3) / eCM);
2485
2486 // Then rotate them
2487 phi = 2. * M_PI * Rndm::flat();
2488 pH[3].rot( theta, phi);
2489 pH[4].rot( theta, phi);
2490
2491 // Set some further info for completeness.
2492 x1H = 1.;
2493 x2H = 1.;
2494 sH = s;
2495 uH = s1 + s2 + s3 + s4 - sH - tH;
2496 mHat = eCM;
2497 p2Abs = pAbs * pAbs;
2498 betaZ = 0.;
2499 pTH = pAbs * sin(theta);
2500
2501 // Done.
2502 return true;
2503
2504}
2505
2506//**************************************************************************
2507
2508// PhaseSpace2to3tauycyl class.
2509// 2 -> 3 kinematics for normal subprocesses.
2510
2511//*********
2512
2513// Constants: could be changed here if desired, but normally should not.
2514// These are of technical nature, as described for each.
2515
2516// Number of Newton-Raphson iterations of kinematics when masses introduced.
2517const int PhaseSpace2to3tauycyl::NITERNR = 5;
2518
2519//*********
2520
2521// Set up for fixed or Breit-Wigner mass selection.
2522
2523bool PhaseSpace2to3tauycyl::setupMasses() {
2524
2525 // Treat Z0 as such or as gamma*/Z0
2526 gmZmode = gmZmodeGlobal;
2527 int gmZmodeProc = sigmaProcessPtr->gmZmode();
2528 if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
2529
2530 // Set sHat limits - based on global limits only.
2531 mHatMin = mHatGlobalMin;
2532 sHatMin = mHatMin*mHatMin;
2533 mHatMax = eCM;
2534 if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax);
2535 sHatMax = mHatMax*mHatMax;
2536
2537 // Masses and widths of resonances.
2538 setupMass1(3);
2539 setupMass1(4);
2540 setupMass1(5);
2541
2542 // Reduced mass range - do not make it as fancy as in two-body case.
2543 if (useBW[3]) mUpper[3] -= (mPeak[4] + mPeak[5]);
2544 if (useBW[4]) mUpper[4] -= (mPeak[3] + mPeak[5]);
2545 if (useBW[5]) mUpper[5] -= (mPeak[3] + mPeak[4]);
2546
2547 // If closed phase space then unallowed process.
2548 bool physical = true;
2549 if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false;
2550 if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false;
2551 if (useBW[5] && mUpper[5] < mLower[5] + MASSMARGIN) physical = false;
2552 if (!useBW[3] && !useBW[4] && !useBW[5] && mHatMax < mPeak[3]
2553 + mPeak[4] + mPeak[5] + MASSMARGIN) physical = false;
2554 if (!physical) return false;
2555
2556 // No extra pT precautions in massless limit - assumed fixed by ME's.
2557 pTHatMin = pTHatGlobalMin;
2558 pT2HatMin = pTHatMin * pTHatMin;
2559 pTHatMax = pTHatGlobalMax;
2560 pT2HatMax = pTHatMax * pTHatMax;
2561
2562 // Prepare to select m3 by BW + flat + 1/s_3.
2563 if (useBW[3]) {
2564 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
2565 * mWidth[3] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
2566 double distToThreshB = (mHatMax - mPeak[3] - mMin[4] - mMin[5])
2567 / mWidth[3];
2568 double distToThresh = min( distToThreshA, distToThreshB);
2569 setupMass2(3, distToThresh);
2570 }
2571
2572 // Prepare to select m4 by BW + flat + 1/s_3.
2573 if (useBW[4]) {
2574 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
2575 * mWidth[4] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
2576 double distToThreshB = (mHatMax - mPeak[4] - mMin[3] - mMin[5])
2577 / mWidth[4];
2578 double distToThresh = min( distToThreshA, distToThreshB);
2579 setupMass2(4, distToThresh);
2580 }
2581
2582 // Prepare to select m5 by BW + flat + 1/s_3.
2583 if (useBW[5]) {
2584 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
2585 * mWidth[5] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
2586 double distToThreshB = (mHatMax - mPeak[5] - mMin[3] - mMin[4])
2587 / mWidth[5];
2588 double distToThresh = min( distToThreshA, distToThreshB);
2589 setupMass2(5, distToThresh);
2590 }
2591
2592 // Initialization masses. For now give up when constrained phase space.
2593 m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3];
2594 m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4];
2595 m5 = (useBW[5]) ? min(mPeak[5], mUpper[5]) : mPeak[5];
2596 if (m3 + m4 + m5 + MASSMARGIN > mHatMax) physical = false;
2597 s3 = m3*m3;
2598 s4 = m4*m4;
2599 s5 = m5*m5;
2600
2601 // Correct selected mass-spectrum to running-width Breit-Wigner.
2602 // Extra safety margin for maximum search.
2603 wtBW = 1.;
2604 if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX;
2605 if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX;
2606 if (useBW[5]) wtBW *= weightMass(5) * EXTRABWWTMAX;
2607
2608 // Done.
2609 return physical;
2610
2611}
2612
2613//*********
2614
2615// Select Breit-Wigner-distributed or fixed masses.
2616
2617bool PhaseSpace2to3tauycyl::trialMasses() {
2618
2619 // By default vanishing cross section.
2620 sigmaNw = 0.;
2621 wtBW = 1.;
2622
2623 // Pick m3, m4 and m5 independently.
2624 trialMass(3);
2625 trialMass(4);
2626 trialMass(5);
2627
2628 // If outside phase space then reject event.
2629 if (m3 + m4 + m5 + MASSMARGIN > mHatMax) return false;
2630
2631 // Correct selected mass-spectrum to running-width Breit-Wigner.
2632 if (useBW[3]) wtBW *= weightMass(3);
2633 if (useBW[4]) wtBW *= weightMass(4);
2634 if (useBW[5]) wtBW *= weightMass(5);
2635
2636 // Done.
2637 return true;
2638
2639}
2640
2641//*********
2642
2643// Construct the four-vector kinematics from the trial values.
2644
2645bool PhaseSpace2to3tauycyl::finalKin() {
2646
2647 // Assign masses to particles assumed massless in matrix elements.
2648 int id3 = sigmaProcessPtr->id(3);
2649 int id4 = sigmaProcessPtr->id(4);
2650 int id5 = sigmaProcessPtr->id(5);
2651 if (idMass[3] == 0) { m3 = ParticleDataTable::m0(id3); s3 = m3*m3; }
2652 if (idMass[4] == 0) { m4 = ParticleDataTable::m0(id4); s4 = m4*m4; }
2653 if (idMass[5] == 0) { m5 = ParticleDataTable::m0(id5); s5 = m5*m5; }
2654
2655 // Check that phase space still open after new mass assignment.
2656 if (m3 + m4 + m5 + MASSMARGIN > mHat) {
2657 infoPtr->errorMsg("Warning in PhaseSpace2to3tauycyl::finalKin: "
2658 "failed after mass assignment");
2659 return false;
2660 }
2661
2662 // Particle masses; incoming always on mass shell.
2663 mH[1] = 0.;
2664 mH[2] = 0.;
2665 mH[3] = m3;
2666 mH[4] = m4;
2667 mH[5] = m5;
2668
2669 // Incoming partons along beam axes.
2670 pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
2671 pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
2672
2673 // Begin three-momentum rescaling to compensate for masses.
2674 if (idMass[3] == 0 || idMass[4] == 0 || idMass[5] == 0) {
2675 double p3S = p3cm.pAbs2();
2676 double p4S = p4cm.pAbs2();
2677 double p5S = p5cm.pAbs2();
2678 double fac = 1.;
2679 double e3, e4, e5, value, deriv;
2680
2681 // Iterate rescaling solution five times, using Newton-Raphson.
2682 for (int i = 0; i < NITERNR; ++i) {
2683 e3 = sqrt(s3 + fac * p3S);
2684 e4 = sqrt(s4 + fac * p4S);
2685 e5 = sqrt(s5 + fac * p5S);
2686 value = e3 + e4 + e5 - mHat;
2687 deriv = 0.5 * (p3S / e3 + p4S / e4 + p5S / e5);
2688 fac -= value / deriv;
2689 }
2690
2691 // Rescale momenta appropriately.
2692 double facRoot = sqrt(fac);
2693 p3cm.rescale3( facRoot );
2694 p4cm.rescale3( facRoot );
2695 p5cm.rescale3( facRoot );
2696 p3cm.e( sqrt(s3 + fac * p3S) );
2697 p4cm.e( sqrt(s4 + fac * p4S) );
2698 p5cm.e( sqrt(s5 + fac * p5S) );
2699 }
2700
2701 // Outgoing partons initially in collision CM frame along beam axes.
2702 pH[3] = p3cm;
2703 pH[4] = p4cm;
2704 pH[5] = p5cm;
2705
2706 // Then boost them to overall CM frame
2707 betaZ = (x1H - x2H)/(x1H + x2H);
2708 pH[3].rot( theta, phi);
2709 pH[4].rot( theta, phi);
2710 pH[3].bst( 0., 0., betaZ);
2711 pH[4].bst( 0., 0., betaZ);
2712 pH[5].bst( 0., 0., betaZ);
2713
2714 // Store average pT of three final particles for documentation.
2715 pTH = (p3cm.pT() + p4cm.pT() + p5cm.pT()) / 3.;
2716
2717 // Done.
2718 return true;
2719}
2720
2721//**************************************************************************
2722
2723// The PhaseSpaceLHA class.
2724// A derived class for Les Houches events.
2725// Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
2726
2727//*********
2728
2729// Constants: could be changed here if desired, but normally should not.
2730// These are of technical nature, as described for each.
2731
2732// LHA convention with cross section in pb forces conversion to mb.
2733const double PhaseSpaceLHA::CONVERTPB2MB = 1e-9;
2734
2735//*********
2736
2737// Find maximal cross section for comparison with internal processes.
2738
2739bool PhaseSpaceLHA::setupSampling() {
2740
2741 // Find which strategy Les Houches events are produced with.
2742 strategy = lhaUpPtr->strategy();
2743 stratAbs = abs(strategy);
2744 if (strategy == 0 || stratAbs > 4) {
2745 ostringstream stratCode;
2746 stratCode << strategy;
2747 infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: unknown "
2748 "Les Houches Accord weighting stategy", stratCode.str());
2749 return false;
2750 }
2751
2752 // Number of contributing processes.
2753 nProc = lhaUpPtr->sizeProc();
2754
2755 // Loop over all processes. Read out maximum and cross section.
2756 xMaxAbsSum = 0.;
2757 xSecSgnSum = 0.;
2758 int idPr;
2759 double xMax, xSec, xMaxAbs;
2760 for (int iProc = 0 ; iProc < nProc; ++iProc) {
2761 idPr = lhaUpPtr->idProcess(iProc);
2762 xMax = lhaUpPtr->xMax(iProc);
2763 xSec = lhaUpPtr->xSec(iProc);
2764
2765 // Check for inconsistencies between strategy and stored values.
2766 if ( (strategy == 1 || strategy == 2) && xMax < 0.) {
2767 infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: "
2768 "negative maximum not allowed");
2769 return false;
2770 }
2771 if ( ( strategy == 2 || strategy == 3) && xSec < 0.) {
2772 infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: "
2773 "negative cross section not allowed");
2774 return false;
2775 }
2776
2777 // Store maximal cross sections for later choice.
2778 if (stratAbs == 1) xMaxAbs = abs(xMax);
2779 else if (stratAbs < 4) xMaxAbs = abs(xSec);
2780 else xMaxAbs = 1.;
2781 idProc.push_back( idPr );
2782 xMaxAbsProc.push_back( xMaxAbs );
2783
2784 // Find sum and convert to mb.
2785 xMaxAbsSum += xMaxAbs;
2786 xSecSgnSum += xSec;
2787 }
2788 sigmaMx = xMaxAbsSum * CONVERTPB2MB;
2789 sigmaSgn = xSecSgnSum * CONVERTPB2MB;
2790
2791 // Done.
2792 return true;
2793
2794}
2795
2796//*********
2797
2798// Construct the next process, by interface to Les Houches class.
2799
2800bool PhaseSpaceLHA::trialKin( bool, bool repeatSame ) {
2801
2802 // Must select process type in some cases.
2803 int idProcNow = 0;
2804 if (repeatSame) idProcNow = idProcSave;
2805 else if (stratAbs <= 2) {
2806 double xMaxAbsRndm = xMaxAbsSum * Rndm::flat();
2807 int iProc = -1;
2808 do xMaxAbsRndm -= xMaxAbsProc[++iProc];
2809 while (xMaxAbsRndm > 0. && iProc < nProc - 1);
2810 idProcNow = idProc[iProc];
2811 }
2812
2813 // Generate Les Houches event. Return if fail (= end of file).
2814 bool physical = lhaUpPtr->setEvent(idProcNow);
2815 if (!physical) return false;
2816
2817 // Find which process was generated.
2818 int idPr = lhaUpPtr->idProcess();
2819 int iProc = 0;
2820 for (int iP = 0; iP < int(idProc.size()); ++iP)
2821 if (idProc[iP] == idPr) iProc = iP;
2822 idProcSave = idPr;
2823
2824 // Extract cross section and rescale according to strategy.
2825 double wtPr = lhaUpPtr->weight();
2826 if (stratAbs == 1) sigmaNw = wtPr * CONVERTPB2MB
2827 * xMaxAbsSum / xMaxAbsProc[iProc];
2828 else if (stratAbs == 2) sigmaNw = (wtPr / abs(lhaUpPtr->xMax(iProc)))
2829 * sigmaMx;
2830 else if (strategy == 3) sigmaNw = sigmaMx;
2831 else if (strategy == -3 && wtPr > 0.) sigmaNw = sigmaMx;
2832 else if (strategy == -3) sigmaNw = -sigmaMx;
2833 else if (stratAbs == 4) sigmaNw = wtPr * CONVERTPB2MB;
2834
2835 // Done.
2836 return true;
2837
2838}
2839
2840//**************************************************************************
2841
2842} // end namespace Pythia8