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