1 // PhaseSpace.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2008 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 // Function definitions (not found in the header) for the
7 // PhaseSpace and PhaseSpace2to2tauyz classes.
9 #include "PhaseSpace.h"
13 //**************************************************************************
15 // The PhaseSpace class.
16 // Base class for phase space generators.
20 // Constants: could be changed here if desired, but normally should not.
21 // These are of technical nature, as described for each.
23 // Number of trial maxima around which maximum search is performed.
24 const int PhaseSpace::NMAXTRY = 2;
26 // Number of three-body trials in phase space optimization.
27 const int PhaseSpace::NTRY3BODY = 20;
29 // Maximum cross section increase, just in case true maximum not found.
30 const double PhaseSpace::SAFETYMARGIN = 1.05;
32 // Small number to avoid division by zero.
33 const double PhaseSpace::TINY = 1e-20;
35 // Fraction of total weight that is shared evenly between all shapes.
36 const double PhaseSpace::EVENFRAC = 0.4;
38 // Two cross sections with a small relative error are assumed same.
39 const double PhaseSpace::SAMESIGMA = 1e-6;
41 // Do not include resonances peaked too far outside allowed mass region.
42 const double PhaseSpace::WIDTHMARGIN = 20.;
44 // Special optimization treatment when two resonances at almost same mass.
45 const double PhaseSpace::SAMEMASS = 0.01;
47 // Minimum phase space left when kinematics constraints are combined.
48 const double PhaseSpace::MASSMARGIN = 0.01;
50 // When using Breit-Wigners in 2 -> 2 raise maximum weight estimate.
51 const double PhaseSpace::EXTRABWWTMAX = 1.25;
53 // Size of Breit-Wigner threshold region, for mass selection biasing.
54 const double PhaseSpace::THRESHOLDSIZE = 3.;
56 // Step size in optimal-mass search, for mass selection biasing.
57 const double PhaseSpace::THRESHOLDSTEP = 0.2;
59 // Minimal rapidity range for allowed open range (in 2 -> 3).
60 const double PhaseSpace::YRANGEMARGIN = 1e-6;
62 // Cutoff for f_e^e at x < 1 - 10^{-10} to be used in phase space selection.
63 // Note: the ...MIN quantities come from 1 - x_max or 1 - tau_max.
64 const double PhaseSpace::LEPTONXMIN = 1e-10;
65 const double PhaseSpace::LEPTONXMAX = 1. - 1e-10;
66 const double PhaseSpace::LEPTONXLOGMIN = log(1e-10);
67 const double PhaseSpace::LEPTONXLOGMAX = log(1. - 1e-10);
68 const double PhaseSpace::LEPTONTAUMIN = 2e-10;
70 // Safety to avoid division with unreasonably small value for z selection.
71 const double PhaseSpace::SHATMINZ = 1.;
73 // Regularization for small pT2min in z = cos(theta) selection.
74 const double PhaseSpace::PT2RATMINZ = 0.0001;
76 // These numbers are hardwired empirical parameters,
77 // intended to speed up the M-generator.
78 const double PhaseSpace::WTCORRECTION[11] = { 1., 1., 1.,
79 2., 5., 15., 60., 250., 1250., 7000., 50000. };
83 // Perform simple initialization and store pointers.
85 void PhaseSpace::init(SigmaProcess* sigmaProcessPtrIn, Info* infoPtrIn,
86 BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
87 SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn) {
89 // Store input pointers for future use.
90 sigmaProcessPtr = sigmaProcessPtrIn;
92 beamAPtr = beamAPtrIn;
93 beamBPtr = beamBPtrIn;
94 sigmaTotPtr = sigmaTotPtrIn;
95 userHooksPtr = userHooksPtrIn;
97 // Some commonly used beam information.
102 eCM = infoPtr->eCM();
105 // Flag if lepton beams, and if non-resolved ones.
106 hasLeptonBeams = ( beamAPtr->isLepton() || beamBPtr->isLepton() );
107 hasPointLeptons = ( hasLeptonBeams
108 && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) );
110 // Standard phase space cuts.
111 mHatGlobalMin = Settings::parm("PhaseSpace:mHatMin");
112 mHatGlobalMax = Settings::parm("PhaseSpace:mHatMax");
113 pTHatGlobalMin = Settings::parm("PhaseSpace:pTHatMin");
114 pTHatGlobalMax = Settings::parm("PhaseSpace:pTHatMax");
115 pTHatMinDiverge = Settings::parm("PhaseSpace:pTHatMinDiverge");
117 // When to use Breit-Wigners.
118 useBreitWigners = Settings::flag("PhaseSpace:useBreitWigners");
119 minWidthBreitWigners = Settings::parm("PhaseSpace:minWidthBreitWigners");
121 // Whether generation is with variable energy.
122 doEnergySpread = Settings::flag("Beams:allowMomentumSpread");
124 // Print flag for maximization information.
125 showSearch = Settings::flag("PhaseSpace:showSearch");
126 showViolation = Settings::flag("PhaseSpace:showViolation");
128 // Know whether a Z0 is pure Z0 or admixed with gamma*.
129 gmZmodeGlobal = Settings::mode("WeakZ0:gmZmode");
131 // Default event-specific kinematics properties.
151 // Default cross section information.
157 // Flag if user should be allow to reweight cross section.
158 canModifySigma = (userHooksPtr > 0)
159 ? userHooksPtr->canModifySigma() : false;
165 // Allow for nonisotropic decays when ME's available.
167 void PhaseSpace::decayKinematics( Event& process) {
169 // Identify sets of sister partons.
171 for (int iResBeg = 5; iResBeg < process.size(); ++iResBeg) {
172 if (iResBeg <= iResEnd) continue;
174 while ( iResEnd < process.size() - 1
175 && process[iResEnd + 1].mother1() == process[iResBeg].mother1()
176 && process[iResEnd + 1].mother2() == process[iResBeg].mother2() )
179 // Check that at least one of them is a resonance.
181 for (int iRes = iResBeg; iRes <= iResEnd; ++iRes)
182 if ( !process[iRes].isFinal() ) hasRes = true;
183 if ( !hasRes ) continue;
185 // Evaluate matrix element and decide whether to keep kinematics.
186 double decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd);
187 if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
188 "Kinematics: negative angular weight");
189 if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
190 "Kinematics: angular weight above unity");
191 while (decWt < Rndm::flat() ) {
193 // Find resonances for which to redo decay angles.
194 for (int iRes = iResBeg; iRes < process.size(); ++iRes) {
195 if ( process[iRes].isFinal() ) continue;
196 int iResMother = iRes;
197 while (iResMother > iResEnd)
198 iResMother = process[iResMother].mother1();
199 if (iResMother < iResBeg) continue;
201 // Do decay of this mother isotropically in phase space.
202 decayKinematicsStep( process, iRes);
204 // End loop over resonance decay chains.
207 // Ready to allow new test of matrix element.
208 decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd);
209 if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
210 "Kinematics: negative angular weight");
211 if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
212 "Kinematics: angular weight above unity");
215 // End loop over sets of sister resonances/partons.
222 // Reselect decay products momenta isotropically in phase space.
223 // Does not redo secondary vertex position!
225 void PhaseSpace::decayKinematicsStep( Event& process, int iRes) {
227 // Multiplicity and mother mass and four-momentum.
228 int i1 = process[iRes].daughter1();
229 int mult = process[iRes].daughter2() + 1 - i1;
230 double m0 = process[iRes].m();
231 Vec4 pRes = process[iRes].p();
233 // Description of two-body decays as simple special case.
236 // Products and product masses.
238 double m1t = process[i1].m();
239 double m2t = process[i2].m();
241 // Energies and absolute momentum in the rest frame.
242 double e1 = 0.5 * (m0*m0 + m1t*m1t - m2t*m2t) / m0;
243 double e2 = 0.5 * (m0*m0 + m2t*m2t - m1t*m1t) / m0;
244 double p12 = 0.5 * sqrtpos( (m0 - m1t - m2t) * (m0 + m1t + m2t)
245 * (m0 + m1t - m2t) * (m0 - m1t + m2t) ) / m0;
247 // Pick isotropic angles to give three-momentum.
248 double cosTheta = 2. * Rndm::flat() - 1.;
249 double sinTheta = sqrt(1. - cosTheta*cosTheta);
250 double phi12 = 2. * M_PI * Rndm::flat();
251 double pX = p12 * sinTheta * cos(phi12);
252 double pY = p12 * sinTheta * sin(phi12);
253 double pZ = p12 * cosTheta;
255 // Fill four-momenta in mother rest frame and then boost to lab frame.
256 Vec4 p1( pX, pY, pZ, e1);
257 Vec4 p2( -pX, -pY, -pZ, e2);
261 // Done for two-body decay.
267 // Description of three-body decays as semi-simple special case.
270 // Products and product masses.
273 double m1t = process[i1].m();
274 double m2t = process[i2].m();
275 double m3t = process[i3].m();
276 double mDiff = m0 - (m1t + m2t + m3t);
278 // Kinematical limits for 2+3 mass. Maximum phase-space weight.
279 double m23Min = m2t + m3t;
280 double m23Max = m0 - m1t;
281 double p1Max = 0.5 * sqrtpos( (m0 - m1t - m23Min)
282 * (m0 + m1t + m23Min) * (m0 + m1t - m23Min)
283 * (m0 - m1t + m23Min) ) / m0;
284 double p23Max = 0.5 * sqrtpos( (m23Max - m2t - m3t)
285 * (m23Max + m2t + m3t) * (m23Max + m2t - m3t)
286 * (m23Max - m2t + m3t) ) / m23Max;
287 double wtPSmax = 0.5 * p1Max * p23Max;
289 // Pick an intermediate mass m23 flat in the allowed range.
290 double wtPS, m23, p1Abs, p23Abs;
292 m23 = m23Min + Rndm::flat() * mDiff;
294 // Translate into relative momenta and find phase-space weight.
295 p1Abs = 0.5 * sqrtpos( (m0 - m1t - m23) * (m0 + m1t + m23)
296 * (m0 + m1t - m23) * (m0 - m1t + m23) ) / m0;
297 p23Abs = 0.5 * sqrtpos( (m23 - m2t - m3t) * (m23 + m2t + m3t)
298 * (m23 + m2t - m3t) * (m23 - m2t + m3t) ) / m23;
299 wtPS = p1Abs * p23Abs;
301 // If rejected, try again with new invariant masses.
302 } while ( wtPS < Rndm::flat() * wtPSmax );
304 // Set up m23 -> m2 + m3 isotropic in its rest frame.
305 double cosTheta = 2. * Rndm::flat() - 1.;
306 double sinTheta = sqrt(1. - cosTheta*cosTheta);
307 double phi23 = 2. * M_PI * Rndm::flat();
308 double pX = p23Abs * sinTheta * cos(phi23);
309 double pY = p23Abs * sinTheta * sin(phi23);
310 double pZ = p23Abs * cosTheta;
311 double e2 = sqrt( m2t*m2t + p23Abs*p23Abs);
312 double e3 = sqrt( m3t*m3t + p23Abs*p23Abs);
313 Vec4 p2( pX, pY, pZ, e2);
314 Vec4 p3( -pX, -pY, -pZ, e3);
316 // Set up 0 -> 1 + 23 isotropic in its rest frame.
317 cosTheta = 2. * Rndm::flat() - 1.;
318 sinTheta = sqrt(1. - cosTheta*cosTheta);
319 phi23 = 2. * M_PI * Rndm::flat();
320 pX = p1Abs * sinTheta * cos(phi23);
321 pY = p1Abs * sinTheta * sin(phi23);
322 pZ = p1Abs * cosTheta;
323 double e1 = sqrt( m1t*m1t + p1Abs*p1Abs);
324 double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
325 Vec4 p1( pX, pY, pZ, e1);
327 // Boost 2 + 3 to the 0 rest frame and then boost to lab frame.
328 Vec4 p23( -pX, -pY, -pZ, e23);
335 // Done for three-body decay.
342 // Do a multibody decay using the M-generator algorithm.
344 // Set up masses and four-momenta in a vector, with mother in slot 0.
345 vector<double> mProd;
346 mProd.push_back( m0);
347 for (int i = i1; i <= process[iRes].daughter2(); ++i)
348 mProd.push_back( process[i].m() );
350 pProd.push_back( pRes);
352 // Sum of daughter masses.
353 double mSum = mProd[1];
354 for (int i = 2; i <= mult; ++i) mSum += mProd[i];
355 double mDiff = m0 - mSum;
357 // Begin setup of intermediate invariant masses.
359 for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
361 // Calculate the maximum weight in the decay.
362 double wtPSmax = 1. / WTCORRECTION[mult];
363 double mMaxWT = mDiff + mProd[mult];
365 for (int i = mult - 1; i > 0; --i) {
367 mMinWT += mProd[i+1];
368 double mNow = mProd[i];
369 wtPSmax *= 0.5 * sqrtpos( (mMaxWT - mMinWT - mNow)
370 * (mMaxWT + mMinWT + mNow) * (mMaxWT + mMinWT - mNow)
371 * (mMaxWT - mMinWT + mNow) ) / mMaxWT;
374 // Begin loop to find the set of intermediate invariant masses.
375 vector<double> rndmOrd;
380 // Find and order random numbers in descending order.
382 rndmOrd.push_back(1.);
383 for (int i = 1; i < mult - 1; ++i) {
384 double rndm = Rndm::flat();
385 rndmOrd.push_back(rndm);
386 for (int j = i - 1; j > 0; --j) {
387 if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
391 rndmOrd.push_back(0.);
393 // Translate into intermediate masses and find weight.
394 for (int i = mult - 1; i > 0; --i) {
395 mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
396 wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
397 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
398 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
401 // If rejected, try again with new invariant masses.
402 } while ( wtPS < Rndm::flat() * wtPSmax );
404 // Perform two-particle decays in the respective rest frame.
406 pInv.resize(mult + 1);
407 for (int i = 1; i < mult; ++i) {
408 double p12 = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
409 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
410 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
412 // Isotropic angles give three-momentum.
413 double cosTheta = 2. * Rndm::flat() - 1.;
414 double sinTheta = sqrt(1. - cosTheta*cosTheta);
415 double phiLoc = 2. * M_PI * Rndm::flat();
416 double pX = p12 * sinTheta * cos(phiLoc);
417 double pY = p12 * sinTheta * sin(phiLoc);
418 double pZ = p12 * cosTheta;
420 // Calculate energies, fill four-momenta.
421 double eHad = sqrt( mProd[i]*mProd[i] + p12*p12);
422 double eInv = sqrt( mInv[i+1]*mInv[i+1] + p12*p12);
423 pProd.push_back( Vec4( pX, pY, pZ, eHad) );
424 pInv[i+1].p( -pX, -pY, -pZ, eInv);
426 pProd.push_back( pInv[mult] );
428 // Boost decay products to the mother rest frame and on to lab frame.
430 for (int iFrame = mult - 1; iFrame > 0; --iFrame)
431 for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);
433 // Done for multibody decay.
434 for (int i = 1; i < mult; ++i)
435 process[i1 + i - 1].p( pProd[i] );
442 // Determine how 3-body phase space should be sampled.
444 void PhaseSpace::setup3Body() {
446 // Check for massive t-channel propagator particles.
447 int idTchan1 = abs( sigmaProcessPtr->idTchan1() );
448 int idTchan2 = abs( sigmaProcessPtr->idTchan2() );
449 mTchan1 = (idTchan1 == 0) ? 0. : ParticleDataTable::m0(idTchan1);
450 mTchan2 = (idTchan2 == 0) ? 0. : ParticleDataTable::m0(idTchan2);
451 sTchan1 = mTchan1 * mTchan1;
452 sTchan2 = mTchan2 * mTchan2;
454 // Find coefficients of different pT2 selection terms. Mirror choice.
455 frac3Pow1 = sigmaProcessPtr->tChanFracPow1();
456 frac3Pow2 = sigmaProcessPtr->tChanFracPow2();
457 frac3Flat = 1. - frac3Pow1 - frac3Pow2;
458 useMirrorWeight = sigmaProcessPtr->useMirrorWeight();
464 // Determine how phase space should be sampled.
466 bool PhaseSpace::setupSampling123(bool is2, bool is3, ostream& os) {
468 // Optional printout.
469 if (showSearch) os << "\n PYTHIA Optimization printout for "
470 << sigmaProcessPtr->name() << "\n \n" << scientific << setprecision(3);
472 // Check that open range in tau (+ set tauMin, tauMax).
473 if (!limitTau(is2, is3)) return false;
475 // Reset coefficients and matrices of equation system to solve.
476 int binTau[8], binY[8], binZ[8];
477 double vecTau[8], matTau[8][8], vecY[8], matY[8][8], vecZ[8], matZ[8][8];
478 for (int i = 0; i < 8; ++i) {
488 for (int j = 0; j < 8; ++j) {
497 // Number of used coefficients/points for each dimension: tau, y, c.
498 nTau = (hasPointLeptons) ? 1 : 2;
499 nY = (hasPointLeptons) ? 1 : 3;
502 // Identify if any resonances contribute in s-channel.
503 idResA = sigmaProcessPtr->resonanceA();
505 mResA = ParticleDataTable::m0(idResA);
506 GammaResA = ParticleDataTable::mWidth(idResA);
507 if (mHatMin > mResA + WIDTHMARGIN * GammaResA || (mHatMax > 0.
508 && mHatMax < mResA - WIDTHMARGIN * GammaResA) ) idResA = 0;
510 idResB = sigmaProcessPtr->resonanceB();
512 mResB = ParticleDataTable::m0(idResB);
513 GammaResB = ParticleDataTable::mWidth(idResB);
514 if (mHatMin > mResB + WIDTHMARGIN * GammaResB || (mHatMax > 0.
515 && mHatMax < mResB - WIDTHMARGIN * GammaResB) ) idResB = 0;
517 if (idResA == 0 && idResB != 0) {
520 GammaResA = GammaResB;
524 // More sampling in tau if resonances in s-channel.
525 if (idResA !=0 && !hasPointLeptons) {
527 tauResA = mResA * mResA / s;
528 widResA = mResA * GammaResA / s;
530 if (idResB != 0 && !hasPointLeptons) {
532 tauResB = mResB * mResB / s;
533 widResB = mResB * GammaResB / s;
536 // More sampling in tau and y if incoming lepton beams.
537 if (hasLeptonBeams && !hasPointLeptons) {
542 // Special case when both resonances have same mass.
544 if (idResB != 0 && abs(mResA - mResB) < SAMEMASS * (GammaResA + GammaResB))
547 // Default z value and weight required for 2 -> 1. Number of dimensions.
550 int nVar = (is2) ? 3 : 2;
552 // Initial values, to be modified later.
558 // Step through grid in tau. Set limits on y and z generation.
559 for (int iTau = 0; iTau < nTau; ++iTau) {
561 if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6;
562 selectTau( iTau, posTau, is2);
563 if (!limitY()) continue;
564 if (is2 && !limitZ()) continue;
566 // Step through grids in y and z.
567 for (int iY = 0; iY < nY; ++iY) {
569 for (int iZ = 0; iZ < nZ; ++iZ) {
570 if (is2) selectZ( iZ, 0.5);
571 double sigmaTmp = 0.;
573 // 2 -> 1: calculate cross section, weighted by phase-space volume.
575 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
576 sigmaTmp = sigmaProcessPtr->sigmaPDF();
577 sigmaTmp *= wtTau * wtY;
579 // 2 -> 2: calculate cross section, weighted by phase-space volume
580 // and Breit-Wigners for masses
582 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
584 sigmaTmp = sigmaProcessPtr->sigmaPDF();
585 sigmaTmp *= wtTau * wtY * wtZ * wtBW;
587 // 2 -> 3: repeat internal 3-body phase space several times and
588 // keep maximal cross section, weighted by phase-space volume
589 // and Breit-Wigners for masses
591 for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
592 if (!select3Body()) continue;
593 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
594 m3, m4, m5, runBW3H, runBW4H, runBW5H);
595 double sigmaTry = sigmaProcessPtr->sigmaPDF();
596 sigmaTry *= wtTau * wtY * wt3Body * wtBW;
597 if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
601 // Allow possibility for user to modify cross section. (3body??)
602 if (canModifySigma) sigmaTmp
603 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
605 // Check if current maximum exceeded.
606 if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp;
608 // Optional printout. Protect against negative cross sections.
609 if (showSearch) os << " tau =" << setw(11) << tau << " y ="
610 << setw(11) << y << " z =" << setw(11) << z
611 << " sigma =" << setw(11) << sigmaTmp << "\n";
612 if (sigmaTmp < 0.) sigmaTmp = 0.;
614 // Sum up tau cross-section pieces in points used.
615 if (!hasPointLeptons) {
617 vecTau[iTau] += sigmaTmp;
618 matTau[iTau][0] += 1. / intTau0;
619 matTau[iTau][1] += (1. / intTau1) / tau;
621 matTau[iTau][2] += (1. / intTau2) / (tau + tauResA);
622 matTau[iTau][3] += (1. / intTau3)
623 * tau / ( pow2(tau - tauResA) + pow2(widResA) );
626 matTau[iTau][4] += (1. / intTau4) / (tau + tauResB);
627 matTau[iTau][5] += (1. / intTau5)
628 * tau / ( pow2(tau - tauResB) + pow2(widResB) );
630 if (hasLeptonBeams) matTau[iTau][nTau - 1] += (1. / intTau6)
631 * tau / max( LEPTONTAUMIN, 1. - tau);
634 // Sum up y cross-section pieces in points used.
635 if (!hasPointLeptons) {
637 vecY[iY] += sigmaTmp;
638 matY[iY][0] += (yMax / intY01) * (y + yMax);
639 matY[iY][1] += (yMax / intY01) * (yMax - y);
640 matY[iY][2] += (yMax / intY2) / cosh(y);
641 if (hasLeptonBeams) {
642 matY[iY][3] += (yMax / intY34)
643 / max( LEPTONXMIN, 1. - exp( y - yMax) );
644 matY[iY][4] += (yMax / intY34)
645 / max( LEPTONXMIN, 1. - exp(-y - yMax) );
649 // Integrals over z expressions at tauMax, to be used below.
651 double p2AbsMax = 0.25 * (pow2(tauMax * s - s3 - s4)
652 - 4. * s3 * s4) / (tauMax * s);
653 double zMaxMax = sqrtpos( 1. - pT2HatMin / p2AbsMax );
654 double zPosMaxMax = max(ratio34, unity34 + zMaxMax);
655 double zNegMaxMax = max(ratio34, unity34 - zMaxMax);
656 double intZ0Max = 2. * zMaxMax;
657 double intZ12Max = log( zPosMaxMax / zNegMaxMax);
658 double intZ34Max = 1. / zNegMaxMax - 1. / zPosMaxMax;
660 // Sum up z cross-section pieces in points used.
662 vecZ[iZ] += sigmaTmp;
664 matZ[iZ][1] += (intZ0Max / intZ12Max) / zNeg;
665 matZ[iZ][2] += (intZ0Max / intZ12Max) / zPos;
666 matZ[iZ][3] += (intZ0Max / intZ34Max) / pow2(zNeg);
667 matZ[iZ][4] += (intZ0Max / intZ34Max) / pow2(zPos);
670 // End of loops over phase space points.
675 // Fail if no non-vanishing cross sections.
681 // Solve respective equation system for better phase space coefficients.
682 if (!hasPointLeptons) solveSys( nTau, binTau, vecTau, matTau, tauCoef);
683 if (!hasPointLeptons) solveSys( nY, binY, vecY, matY, yCoef);
684 if (is2) solveSys( nZ, binZ, vecZ, matZ, zCoef);
685 if (showSearch) os << "\n";
687 // Provide cumulative sum of coefficients.
688 tauCoefSum[0] = tauCoef[0];
689 yCoefSum[0] = yCoef[0];
690 zCoefSum[0] = zCoef[0];
691 for (int i = 1; i < 8; ++ i) {
692 tauCoefSum[i] = tauCoefSum[i - 1] + tauCoef[i];
693 yCoefSum[i] = yCoefSum[i - 1] + yCoef[i];
694 zCoefSum[i] = zCoefSum[i - 1] + zCoef[i];
696 // The last element should be > 1 to be on safe side in selection below.
697 tauCoefSum[nTau - 1] = 2.;
698 yCoefSum[nY - 1] = 2.;
699 zCoefSum[nZ - 1] = 2.;
702 // Begin find two most promising maxima among same points as before.
703 int iMaxTau[NMAXTRY + 2], iMaxY[NMAXTRY + 2], iMaxZ[NMAXTRY + 2];
704 double sigMax[NMAXTRY + 2];
707 // Scan same grid as before in tau, y, z.
708 for (int iTau = 0; iTau < nTau; ++iTau) {
710 if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6;
711 selectTau( iTau, posTau, is2);
712 if (!limitY()) continue;
713 if (is2 && !limitZ()) continue;
714 for (int iY = 0; iY < nY; ++iY) {
716 for (int iZ = 0; iZ < nZ; ++iZ) {
717 if (is2) selectZ( iZ, 0.5);
718 double sigmaTmp = 0.;
720 // 2 -> 1: calculate cross section, weighted by phase-space volume.
722 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
723 sigmaTmp = sigmaProcessPtr->sigmaPDF();
724 sigmaTmp *= wtTau * wtY;
726 // 2 -> 2: calculate cross section, weighted by phase-space volume
727 // and Breit-Wigners for masses
729 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
731 sigmaTmp = sigmaProcessPtr->sigmaPDF();
732 sigmaTmp *= wtTau * wtY * wtZ * wtBW;
734 // 2 -> 3: repeat internal 3-body phase space several times and
735 // keep maximal cross section, weighted by phase-space volume
736 // and Breit-Wigners for masses
738 for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
739 if (!select3Body()) continue;
740 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
741 m3, m4, m5, runBW3H, runBW4H, runBW5H);
742 double sigmaTry = sigmaProcessPtr->sigmaPDF();
743 sigmaTry *= wtTau * wtY * wt3Body * wtBW;
744 if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
748 // Allow possibility for user to modify cross section. (3body??)
749 if (canModifySigma) sigmaTmp
750 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
752 // Optional printout. Protect against negative cross section.
753 if (showSearch) os << " tau =" << setw(11) << tau << " y ="
754 << setw(11) << y << " z =" << setw(11) << z
755 << " sigma =" << setw(11) << sigmaTmp << "\n";
756 if (sigmaTmp < 0.) sigmaTmp = 0.;
758 // Check that point is not simply mirror of already found one.
759 bool mirrorPoint = false;
760 for (int iMove = 0; iMove < nMax; ++iMove)
761 if (abs(sigmaTmp - sigMax[iMove]) < SAMESIGMA
762 * (sigmaTmp + sigMax[iMove])) mirrorPoint = true;
764 // Add to or insert in maximum list. Only first two count.
767 for (int iMove = nMax - 1; iMove >= -1; --iMove) {
769 if (iInsert == 0 || sigmaTmp < sigMax[iMove]) break;
770 iMaxTau[iMove + 1] = iMaxTau[iMove];
771 iMaxY[iMove + 1] = iMaxY[iMove];
772 iMaxZ[iMove + 1] = iMaxZ[iMove];
773 sigMax[iMove + 1] = sigMax[iMove];
775 iMaxTau[iInsert] = iTau;
778 sigMax[iInsert] = sigmaTmp;
779 if (nMax < NMAXTRY) ++nMax;
782 // Found two most promising maxima.
786 if (showSearch) os << "\n";
788 // Read out starting position for search.
790 int beginVar = (hasPointLeptons) ? 2 : 0;
791 for (int iMax = 0; iMax < nMax; ++iMax) {
792 int iTau = iMaxTau[iMax];
793 int iY = iMaxY[iMax];
794 int iZ = iMaxZ[iMax];
799 double varVal, varNew, deltaVar, marginVar, sigGrid[3];
801 // Starting point and step size in parameter space.
802 for (int iRepeat = 0; iRepeat < 2; ++iRepeat) {
803 // Run through (possibly a subset of) tau, y and z.
804 for (int iVar = beginVar; iVar < nVar; ++iVar) {
805 if (iVar == 0) varVal = tauVal;
806 else if (iVar == 1) varVal = yVal;
808 deltaVar = (iRepeat == 0) ? 0.1
809 : max( 0.01, min( 0.05, min( varVal - 0.02, 0.98 - varVal) ) );
810 marginVar = (iRepeat == 0) ? 0.02 : 0.002;
811 int moveStart = (iRepeat == 0 && iVar == 0) ? 0 : 1;
812 for (int move = moveStart; move < 9; ++move) {
814 // Define new parameter-space point by step in one dimension.
818 } else if (move == 1) {
820 varNew = varVal + deltaVar;
821 } else if (move == 2) {
823 varNew = varVal - deltaVar;
824 } else if (sigGrid[2] >= max( sigGrid[0], sigGrid[1])
825 && varVal + 2. * deltaVar < 1. - marginVar) {
827 sigGrid[0] = sigGrid[1];
828 sigGrid[1] = sigGrid[2];
830 varNew = varVal + deltaVar;
831 } else if (sigGrid[0] >= max( sigGrid[1], sigGrid[2])
832 && varVal - 2. * deltaVar > marginVar) {
834 sigGrid[2] = sigGrid[1];
835 sigGrid[1] = sigGrid[0];
837 varNew = varVal - deltaVar;
838 } else if (sigGrid[2] >= sigGrid[0]) {
841 sigGrid[0] = sigGrid[1];
847 sigGrid[2] = sigGrid[1];
852 // Convert to relevant variables and find derived new limits.
853 bool insideLimits = true;
856 selectTau( iTau, tauVal, is2);
857 if (!limitY()) insideLimits = false;
858 if (is2 && !limitZ()) insideLimits = false;
861 if (is2) selectZ( iZ, zVal);
863 } else if (iVar == 1) {
866 } else if (iVar == 2) {
871 // Evaluate cross-section.
872 double sigmaTmp = 0.;
875 // 2 -> 1: calculate cross section, weighted by phase-space volume.
877 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
878 sigmaTmp = sigmaProcessPtr->sigmaPDF();
879 sigmaTmp *= wtTau * wtY;
881 // 2 -> 2: calculate cross section, weighted by phase-space volume
882 // and Breit-Wigners for masses
884 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
886 sigmaTmp = sigmaProcessPtr->sigmaPDF();
887 sigmaTmp *= wtTau * wtY * wtZ * wtBW;
889 // 2 -> 3: repeat internal 3-body phase space several times and
890 // keep maximal cross section, weighted by phase-space volume
891 // and Breit-Wigners for masses
893 for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
894 if (!select3Body()) continue;
895 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
896 m3, m4, m5, runBW3H, runBW4H, runBW5H);
897 double sigmaTry = sigmaProcessPtr->sigmaPDF();
898 sigmaTry *= wtTau * wtY * wt3Body * wtBW;
899 if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
903 // Allow possibility for user to modify cross section.
904 if (canModifySigma) sigmaTmp
905 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
907 // Optional printout. Protect against negative cross section.
908 if (showSearch) os << " tau =" << setw(11) << tau << " y ="
909 << setw(11) << y << " z =" << setw(11) << z
910 << " sigma =" << setw(11) << sigmaTmp << "\n";
911 if (sigmaTmp < 0.) sigmaTmp = 0.;
914 // Save new maximum. Final maximum.
915 sigGrid[iGrid] = sigmaTmp;
916 if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp;
921 sigmaMx *= SAFETYMARGIN;
923 // Optional printout.
924 if (showSearch) os << "\n Final maximum = " << setw(11) << sigmaMx << endl;
932 // Select a trial kinematics phase space point.
933 // Note: by In is meant the integral over the quantity multiplying
934 // coefficient cn. The sum of cn is normalized to unity.
936 bool PhaseSpace::trialKin123(bool is2, bool is3, bool inEvent, ostream& os) {
938 // Allow for possibility that energy varies from event to event.
939 if (doEnergySpread) {
940 eCM = infoPtr->eCM();
943 // Find shifted tauRes values.
944 if (idResA !=0 && !hasPointLeptons) {
945 tauResA = mResA * mResA / s;
946 widResA = mResA * GammaResA / s;
948 if (idResB != 0 && !hasPointLeptons) {
949 tauResB = mResB * mResB / s;
950 widResB = mResB * GammaResB / s;
954 // Choose tau according to h1(tau)/tau, where
955 // h1(tau) = c0/I0 + (c1/I1) * 1/tau
956 // + (c2/I2) / (tau + tauResA)
957 // + (c3/I3) * tau / ((tau - tauResA)^2 + widResA^2)
958 // + (c4/I4) / (tau + tauResB)
959 // + (c5/I5) * tau / ((tau - tauResB)^2 + widResB^2)
960 // + (c6/I6) * tau / (1 - tau).
961 if (!limitTau(is2, is3)) return false;
963 if (!hasPointLeptons) {
964 double rTau = Rndm::flat();
965 while (rTau > tauCoefSum[iTau]) ++iTau;
967 selectTau( iTau, Rndm::flat(), is2);
969 // Choose y according to h2(y), where
970 // h2(y) = (c0/I0) * (y-ymin) + (c1/I1) * (ymax-y)
971 // + (c2/I2) * 1/cosh(y) + (c3/I3) * 1 / (1 - exp(y-ymax))
972 // + (c4/I4) * 1 / (1 - exp(ymin-y)).
973 if (!limitY()) return false;
975 if (!hasPointLeptons) {
976 double rY = Rndm::flat();
977 while (rY > yCoefSum[iY]) ++iY;
979 selectY( iY, Rndm::flat());
981 // Choose z = cos(thetaHat) according to h3(z), where
982 // h3(z) = c0/I0 + (c1/I1) * 1/(A - z) + (c2/I2) * 1/(A + z)
983 // + (c3/I3) * 1/(A - z)^2 + (c4/I4) * 1/(A + z)^2,
984 // where A = 1 + 2*(m3*m4/sH)^2 (= 1 for massless products).
986 if (!limitZ()) return false;
988 double rZ = Rndm::flat();
989 while (rZ > zCoefSum[iZ]) ++iZ;
990 selectZ( iZ, Rndm::flat());
993 // 2 -> 1: calculate cross section, weighted by phase-space volume.
995 sigmaProcessPtr->set1Kin( x1H, x2H, sH);
996 sigmaNw = sigmaProcessPtr->sigmaPDF();
997 sigmaNw *= wtTau * wtY;
999 // 2 -> 2: calculate cross section, weighted by phase-space volume
1000 // and Breit-Wigners for masses
1002 sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, runBW3H, runBW4H);
1003 sigmaNw = sigmaProcessPtr->sigmaPDF();
1004 sigmaNw *= wtTau * wtY * wtZ * wtBW;
1006 // 2 -> 3: also sample internal 3-body phase, weighted by
1007 // 2 -> 1 phase-space volume and Breit-Wigners for masses
1009 if (!select3Body()) sigmaNw = 0.;
1011 sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
1012 m3, m4, m5, runBW3H, runBW4H, runBW5H);
1013 sigmaNw = sigmaProcessPtr->sigmaPDF();
1014 sigmaNw *= wtTau * wtY * wt3Body * wtBW;
1018 // Allow possibility for user to modify cross section.
1019 if (canModifySigma) sigmaNw
1020 *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, inEvent);
1022 // Check if maximum violated.
1024 if (sigmaNw > sigmaMx) {
1025 infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin: "
1026 "maximum for cross section violated");
1027 double violFact = SAFETYMARGIN * sigmaNw / sigmaMx;
1028 sigmaMx = SAFETYMARGIN * sigmaNw;
1031 // Optional printout of (all) violations.
1032 if (showViolation) {
1033 if (violFact < 9.99) os << fixed;
1034 else os << scientific;
1035 os << " PYTHIA Maximum for " << sigmaProcessPtr->name()
1036 << " increased by factor " << setprecision(3) << violFact
1037 << " to " << scientific << sigmaMx << endl;
1041 // Check if negative cross section.
1042 if (sigmaNw < sigmaNeg) {
1043 infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin:"
1044 " negative cross section set 0", "for " + sigmaProcessPtr->name() );
1047 // Optional printout of (all) violations.
1048 if (showViolation) os << " PYTHIA Negative minimum for "
1049 << sigmaProcessPtr->name() << " changed to " << scientific
1050 << setprecision(3) << sigmaNeg << endl;
1052 if (sigmaNw < 0.) sigmaNw = 0.;
1060 // Find range of allowed tau values.
1062 bool PhaseSpace::limitTau(bool is2, bool is3) {
1064 // Trivial reply for unresolved lepton beams.
1065 if (hasPointLeptons) {
1071 // Requirements from allowed mHat range.
1072 tauMin = sHatMin / s;
1073 tauMax = (mHatMax < mHatMin) ? 1. : min( 1., sHatMax / s);
1075 // Requirements from allowed pT range and masses.
1077 double mT3Min = sqrt(s3 + pT2HatMin);
1078 double mT4Min = sqrt(s4 + pT2HatMin);
1079 double mT5Min = (is3) ? sqrt(s5 + pT2HatMin) : 0.;
1080 tauMin = max( tauMin, pow2(mT3Min + mT4Min + mT5Min) / s);
1083 // Check that there is an open range.
1084 return (tauMax > tauMin);
1089 // Find range of allowed y values.
1091 bool PhaseSpace::limitY() {
1093 // Trivial reply for unresolved lepton beams.
1094 if (hasPointLeptons) {
1099 // Requirements from selected tau value.
1100 yMax = -0.5 * log(tau);
1102 // For lepton beams requirements from cutoff for f_e^e.
1103 double yMaxMargin = (hasLeptonBeams) ? yMax + LEPTONXLOGMAX : yMax;
1105 // Check that there is an open range.
1106 return (yMaxMargin > 0.);
1111 // Find range of allowed z = cos(theta) values.
1113 bool PhaseSpace::limitZ() {
1119 // Requirements from pTHat limits.
1120 zMax = sqrtpos( 1. - pT2HatMin / p2Abs );
1121 if (pTHatMax > pTHatMin) zMin = sqrtpos( 1. - pT2HatMax / p2Abs );
1123 // Check that there is an open range.
1124 return (zMax > zMin);
1129 // Select tau according to a choice of shapes.
1131 void PhaseSpace::selectTau(int iTau, double tauVal, bool is2) {
1133 // Trivial reply for unresolved lepton beams.
1134 if (hasPointLeptons) {
1140 p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1141 pAbs = sqrtpos( p2Abs );
1146 // Contributions from s-channel resonances.
1151 tRatA = ((tauResA + tauMax) / (tauResA + tauMin)) * (tauMin / tauMax);
1152 aLowA = atan( (tauMin - tauResA) / widResA);
1153 aUppA = atan( (tauMax - tauResA) / widResA);
1159 tRatB = ((tauResB + tauMax) / (tauResB + tauMin)) * (tauMin / tauMax);
1160 aLowB = atan( (tauMin - tauResB) / widResB);
1161 aUppB = atan( (tauMax - tauResB) / widResB);
1164 // Contributions from 1 / (1 - tau) for lepton beams.
1167 if (hasLeptonBeams) {
1168 aLowT = log( max( LEPTONTAUMIN, 1. - tauMin) );
1169 aUppT = log( max( LEPTONTAUMIN, 1. - tauMax) );
1170 intTau6 = aLowT - aUppT;
1173 // Select according to 1/tau or 1/tau^2.
1174 if (iTau == 0) tau = tauMin * pow( tauMax / tauMin, tauVal);
1175 else if (iTau == 1) tau = tauMax * tauMin
1176 / (tauMin + (tauMax - tauMin) * tauVal);
1178 // Select according to 1 / (1 - tau) for lepton beams.
1179 else if (hasLeptonBeams && iTau == nTau - 1)
1180 tau = 1. - exp( aUppT + intTau6 * tauVal );
1182 // Select according to 1 / (tau * (tau + tauRes)) or
1183 // 1 / ((tau - tauRes)^2 + widRes^2) for resonances A and B.
1184 else if (iTau == 2) tau = tauResA * tauMin
1185 / ((tauResA + tauMin) * pow( tRatA, tauVal) - tauMin);
1186 else if (iTau == 3) tau = tauResA + widResA
1187 * tan( aLowA + (aUppA - aLowA) * tauVal);
1188 else if (iTau == 4) tau = tauResB * tauMin
1189 / ((tauResB + tauMin) * pow( tRatB, tauVal) - tauMin);
1190 else if (iTau == 5) tau = tauResB + widResB
1191 * tan( aLowB + (aUppB - aLowB) * tauVal);
1193 // Phase space weight in tau.
1194 intTau0 = log( tauMax / tauMin);
1195 intTau1 = (tauMax - tauMin) / (tauMax * tauMin);
1196 double invWtTau = (tauCoef[0] / intTau0) + (tauCoef[1] / intTau1) / tau;
1198 intTau2 = -log(tRatA) / tauResA;
1199 intTau3 = (aUppA - aLowA) / widResA;
1200 invWtTau += (tauCoef[2] / intTau2) / (tau + tauResA)
1201 + (tauCoef[3] / intTau3) * tau / ( pow2(tau - tauResA) + pow2(widResA) );
1204 intTau4 = -log(tRatB) / tauResB;
1205 intTau5 = (aUppB - aLowB) / widResB;
1206 invWtTau += (tauCoef[4] / intTau4) / (tau + tauResB)
1207 + (tauCoef[5] / intTau5) * tau / ( pow2(tau - tauResB) + pow2(widResB) );
1210 invWtTau += (tauCoef[nTau - 1] / intTau6)
1211 * tau / max( LEPTONTAUMIN, 1. - tau);
1212 wtTau = 1. / invWtTau;
1214 // Calculate sHat and absolute momentum of outgoing partons.
1218 p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1219 pAbs = sqrtpos( p2Abs );
1226 // Select y according to a choice of shapes.
1228 void PhaseSpace::selectY(int iY, double yVal) {
1230 // Trivial reply for unresolved lepton beams.
1231 if (hasPointLeptons) {
1239 // Standard expressions used below.
1240 double atanMax = atan( exp(yMax) );
1241 double atanMin = atan( exp(-yMax) );
1242 double aUppY = (hasLeptonBeams)
1243 ? log( max( LEPTONXMIN, LEPTONXMAX / tau - 1. ) ) : 0.;
1244 double aLowY = LEPTONXLOGMIN;
1246 // y - y_min or mirrored y_max - y.
1247 if (iY <= 1) y = yMax * (2. * sqrt(yVal) - 1.);
1251 y = log( tan( atanMin + (atanMax - atanMin) * yVal ) );
1253 // 1 / (1 - exp(y - y_max)) or mirrored 1 / (1 - exp(y_min - y)).
1254 else y = yMax - log( 1. + exp(aLowY + (aUppY - aLowY) * yVal) );
1256 // Mirror two cases.
1257 if (iY == 1 || iY == 4) y = -y;
1259 // Phase space integral in y.
1260 intY01 = 0.5 * pow2(2. * yMax);
1261 intY2 = 2. * (atanMax - atanMin);
1262 intY34 = aUppY - aLowY;
1263 double invWtY = (yCoef[0] / intY01) * (y + yMax)
1264 + (yCoef[1] / intY01) * (yMax - y) + (yCoef[2] / intY2) / cosh(y);
1265 if (hasLeptonBeams) invWtY
1266 += (yCoef[3] / intY34) / max( LEPTONXMIN, 1. - exp( y - yMax) )
1267 + (yCoef[4] / intY34) / max( LEPTONXMIN, 1. - exp(-y - yMax) );
1270 // Calculate x1 and x2.
1271 x1H = sqrt(tau) * exp(y);
1272 x2H = sqrt(tau) * exp(-y);
1277 // Select z = cos(theta) according to a choice of shapes.
1278 // The selection is split in the positive- and negative-z regions,
1279 // since a pTmax cut can remove the region around z = 0.
1281 void PhaseSpace::selectZ(int iZ, double zVal) {
1283 // Mass-dependent dampening of pT -> 0 limit.
1284 ratio34 = max(TINY, 2. * s3 * s4 / pow2(sH));
1285 unity34 = 1. + ratio34;
1286 double ratiopT2 = 2. * pT2HatMin / max( SHATMINZ, sH);
1287 if (ratiopT2 < PT2RATMINZ) ratio34 = max( ratio34, ratiopT2);
1289 // Common expressions in z limits.
1290 double zPosMax = max(ratio34, unity34 + zMax);
1291 double zNegMax = max(ratio34, unity34 - zMax);
1292 double zPosMin = max(ratio34, unity34 + zMin);
1293 double zNegMin = max(ratio34, unity34 - zMin);
1297 if (zVal < 0.5) z = -(zMax + (zMin - zMax) * 2. * zVal);
1298 else z = zMin + (zMax - zMin) * (2. * zVal - 1.);
1300 // 1 / (unity34 - z).
1301 } else if (iZ == 1) {
1302 double areaNeg = log(zPosMax / zPosMin);
1303 double areaPos = log(zNegMin / zNegMax);
1304 double area = areaNeg + areaPos;
1305 if (zVal * area < areaNeg) {
1306 double zValMod = zVal * area / areaNeg;
1307 z = unity34 - zPosMax * pow(zPosMin / zPosMax, zValMod);
1309 double zValMod = (zVal * area - areaNeg)/ areaPos;
1310 z = unity34 - zNegMin * pow(zNegMax / zNegMin, zValMod);
1313 // 1 / (unity34 + z).
1314 } else if (iZ == 2) {
1315 double areaNeg = log(zNegMin / zNegMax);
1316 double areaPos = log(zPosMax / zPosMin);
1317 double area = areaNeg + areaPos;
1318 if (zVal * area < areaNeg) {
1319 double zValMod = zVal * area / areaNeg;
1320 z = zNegMax * pow(zNegMin / zNegMax, zValMod) - unity34;
1322 double zValMod = (zVal * area - areaNeg)/ areaPos;
1323 z = zPosMin * pow(zPosMax / zPosMin, zValMod) - unity34;
1326 // 1 / (unity34 - z)^2.
1327 } else if (iZ == 3) {
1328 double areaNeg = 1. / zPosMin - 1. / zPosMax;
1329 double areaPos = 1. / zNegMax - 1. / zNegMin;
1330 double area = areaNeg + areaPos;
1331 if (zVal * area < areaNeg) {
1332 double zValMod = zVal * area / areaNeg;
1333 z = unity34 - 1. / (1./zPosMax + areaNeg * zValMod);
1335 double zValMod = (zVal * area - areaNeg)/ areaPos;
1336 z = unity34 - 1. / (1./zNegMin + areaPos * zValMod);
1339 // 1 / (unity34 + z)^2.
1340 } else if (iZ == 4) {
1341 double areaNeg = 1. / zNegMax - 1. / zNegMin;
1342 double areaPos = 1. / zPosMin - 1. / zPosMax;
1343 double area = areaNeg + areaPos;
1344 if (zVal * area < areaNeg) {
1345 double zValMod = zVal * area / areaNeg;
1346 z = 1. / (1./zNegMax - areaNeg * zValMod) - unity34;
1348 double zValMod = (zVal * area - areaNeg)/ areaPos;
1349 z = 1. / (1./zPosMin - areaPos * zValMod) - unity34;
1353 // Safety check for roundoff errors. Combinations with z.
1354 if (z < 0.) z = min(-zMin, max(-zMax, z));
1355 else z = min(zMax, max(zMin, z));
1356 zNeg = max(ratio34, unity34 - z);
1357 zPos = max(ratio34, unity34 + z);
1359 // Phase space integral in z.
1360 double intZ0 = 2. * (zMax - zMin);
1361 double intZ12 = log( (zPosMax * zNegMin) / (zPosMin * zNegMax) );
1362 double intZ34 = 1. / zPosMin - 1. / zPosMax + 1. / zNegMax
1364 wtZ = mHat * pAbs / ( (zCoef[0] / intZ0) + (zCoef[1] / intZ12) / zNeg
1365 + (zCoef[2] / intZ12) / zPos + (zCoef[3] / intZ34) / pow2(zNeg)
1366 + (zCoef[4] / intZ34) / pow2(zPos) );
1368 // Calculate tHat and uHat. Also gives pTHat.
1369 double sH34 = -0.5 * (sH - s3 - s4);
1370 tH = sH34 + mHat * pAbs * z;
1371 uH = sH34 - mHat * pAbs * z;
1372 pTH = sqrtpos( (tH * uH - s3 * s4) / sH);
1378 // Select three-body phase space according to a cylindrically based form
1379 // that can be chosen to favour low pT based on the form of propagators.
1381 bool PhaseSpace::select3Body() {
1383 // Upper and lower limits of pT choice for 4 and 5.
1384 double m35S = pow2(m3 + m5);
1385 double pT4Smax = 0.25 * ( pow2(sH - s4 - m35S) - 4. * s4 * m35S ) / sH;
1386 if (pTHatMax > pTHatMin) pT4Smax = min( pT2HatMax, pT4Smax);
1387 double pT4Smin = pT2HatMin;
1388 double m34S = pow2(m3 + m4);
1389 double pT5Smax = 0.25 * ( pow2(sH - s5 - m34S) - 4. * s5 * m34S ) / sH;
1390 if (pTHatMax > pTHatMin) pT5Smax = min( pT2HatMax, pT5Smax);
1391 double pT5Smin = pT2HatMin;
1393 // Check that pT ranges not closed.
1394 if ( pT4Smax < pow2(pTHatMin + MASSMARGIN) ) return false;
1395 if ( pT5Smax < pow2(pTHatMin + MASSMARGIN) ) return false;
1397 // Select pT4S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2.
1398 double pTSmaxProp = pT4Smax + sTchan1;
1399 double pTSminProp = pT4Smin + sTchan1;
1400 double pTSratProp = pTSmaxProp / pTSminProp;
1401 double pTSdiff = pT4Smax - pT4Smin;
1402 double rShape = Rndm::flat();
1404 if (rShape < frac3Flat) pT4S = pT4Smin + Rndm::flat() * pTSdiff;
1405 else if (rShape < frac3Flat + frac3Pow1) pT4S = max( pT2HatMin,
1406 pTSminProp * pow( pTSratProp, Rndm::flat() ) - sTchan1 );
1407 else pT4S = max( pT2HatMin, pTSminProp * pTSmaxProp
1408 / (pTSminProp + Rndm::flat()* pTSdiff) - sTchan1 );
1409 double wt4 = pTSdiff / ( frac3Flat
1410 + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT4S + sTchan1))
1411 + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT4S + sTchan1) );
1413 // Select pT5S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2.
1414 pTSmaxProp = pT5Smax + sTchan2;
1415 pTSminProp = pT5Smin + sTchan2;
1416 pTSratProp = pTSmaxProp / pTSminProp;
1417 pTSdiff = pT5Smax - pT5Smin;
1418 rShape = Rndm::flat();
1420 if (rShape < frac3Flat) pT5S = pT5Smin + Rndm::flat() * pTSdiff;
1421 else if (rShape < frac3Flat + frac3Pow1) pT5S = max( pT2HatMin,
1422 pTSminProp * pow( pTSratProp, Rndm::flat() ) - sTchan2 );
1423 else pT5S = max( pT2HatMin, pTSminProp * pTSmaxProp
1424 / (pTSminProp + Rndm::flat()* pTSdiff) - sTchan2 );
1425 double wt5 = pTSdiff / ( frac3Flat
1426 + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT5S + sTchan2))
1427 + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT5S + sTchan2) );
1429 // Select azimuthal angles and check that third pT in range.
1430 double phi4 = 2. * M_PI * Rndm::flat();
1431 double phi5 = 2. * M_PI * Rndm::flat();
1432 double pT3S = max( 0., pT4S + pT5S + 2. * sqrt(pT4S * pT5S)
1433 * cos(phi4 - phi5) );
1434 if ( pT3S < pT2HatMin || (pTHatMax > pTHatMin && pT3S > pT2HatMax) )
1437 // Calculate transverse masses and check that phase space not closed.
1438 double sT3 = s3 + pT3S;
1439 double sT4 = s4 + pT4S;
1440 double sT5 = s5 + pT5S;
1441 double mT3 = sqrt(sT3);
1442 double mT4 = sqrt(sT4);
1443 double mT5 = sqrt(sT5);
1444 if ( mT3 + mT4 + mT5 + MASSMARGIN > mHat ) return false;
1446 // Select rapidity for particle 3 and check that phase space not closed.
1447 double m45S = pow2(mT4 + mT5);
1448 double y3max = log( ( sH + sT3 - m45S + sqrtpos( pow2(sH - sT3 - m45S)
1449 - 4 * sT3 * m45S ) ) / (2. * mHat * mT3) );
1450 if (y3max < YRANGEMARGIN) return false;
1451 double y3 = (2. * Rndm::flat() - 1.) * (1. - YRANGEMARGIN) * y3max;
1452 double pz3 = mT3 * sinh(y3);
1453 double e3 = mT3 * cosh(y3);
1455 // Find momentum transfers in the two mirror solutions (in 4-5 frame).
1457 double e45 = mHat - e3;
1458 double sT45 = e45 * e45 - pz45 * pz45;
1459 double lam45 = sqrtpos( pow2(sT45 - sT4 - sT5) - 4. * sT4 * sT5 );
1460 if (lam45 < YRANGEMARGIN * sH) return false;
1461 double lam4e = sT45 + sT4 - sT5;
1462 double lam5e = sT45 + sT5 - sT4;
1463 double tFac = -0.5 * mHat / sT45;
1464 double t1Pos = tFac * (e45 - pz45) * (lam4e - lam45);
1465 double t1Neg = tFac * (e45 - pz45) * (lam4e + lam45);
1466 double t2Pos = tFac * (e45 + pz45) * (lam5e - lam45);
1467 double t2Neg = tFac * (e45 + pz45) * (lam5e + lam45);
1469 // Construct relative mirror weights and make choice.
1470 double wtPosUnnorm = 1.;
1471 double wtNegUnnorm = 1.;
1472 if (useMirrorWeight) {
1473 wtPosUnnorm = 1./ pow2( (t1Pos - sTchan1) * (t2Pos - sTchan2) );
1474 wtNegUnnorm = 1./ pow2( (t1Neg - sTchan1) * (t2Neg - sTchan2) );
1476 double wtPos = wtPosUnnorm / (wtPosUnnorm + wtNegUnnorm);
1477 double wtNeg = wtNegUnnorm / (wtPosUnnorm + wtNegUnnorm);
1478 double epsilon = (Rndm::flat() < wtPos) ? 1. : -1.;
1480 // Construct four-vectors in rest frame of subprocess.
1481 double px4 = sqrt(pT4S) * cos(phi4);
1482 double py4 = sqrt(pT4S) * sin(phi4);
1483 double px5 = sqrt(pT5S) * cos(phi5);
1484 double py5 = sqrt(pT5S) * sin(phi5);
1485 double pz4 = 0.5 * (pz45 * lam4e + epsilon * e45 * lam45) / sT45;
1486 double pz5 = pz45 - pz4;
1487 double e4 = sqrt(sT4 + pz4 * pz4);
1488 double e5 = sqrt(sT5 + pz5 * pz5);
1489 p3cm = Vec4( -(px4 + px5), -(py4 + py5), pz3, e3);
1490 p4cm = Vec4( px4, py4, pz4, e4);
1491 p5cm = Vec4( px5, py5, pz5, e5);
1493 // Total weight to associate with kinematics choice.
1494 wt3Body = wt4 * wt5 * (2. * y3max) / (128. * pow3(M_PI) * lam45);
1495 wt3Body *= (epsilon > 0.) ? 1. / wtPos : 1. / wtNeg;
1497 // Cross section of form |M|^2/(2 sHat) dPS_3 so need 1/(2 sHat).
1498 wt3Body /= (2. * sH);
1507 // Solve linear equation system for better phase space coefficients.
1509 void PhaseSpace::solveSys( int n, int bin[8], double vec[8],
1510 double mat[8][8], double coef[8], ostream& os) {
1512 // Optional printout.
1514 os << "\n Equation system: " << setw(5) << bin[0];
1515 for (int j = 0; j < n; ++j) os << setw(12) << mat[0][j];
1516 os << setw(12) << vec[0] << "\n";
1517 for (int i = 1; i < n; ++i) {
1518 os << " " << setw(5) << bin[i];
1519 for (int j = 0; j < n; ++j) os << setw(12) << mat[i][j];
1520 os << setw(12) << vec[i] << "\n";
1525 double vecNor[8], coefTmp[8];
1526 for (int i = 0; i < n; ++i) coefTmp[i] = 0.;
1528 // Check if equation system solvable.
1529 bool canSolve = true;
1530 for (int i = 0; i < n; ++i) if (bin[i] == 0) canSolve = false;
1532 for (int i = 0; i < n; ++i) vecSum += vec[i];
1533 if (abs(vecSum) < TINY) canSolve = false;
1535 // Solve to find relative importance of cross-section pieces.
1537 for (int i = 0; i < n; ++i) vecNor[i] = max( 0.1, vec[i] / vecSum);
1538 for (int k = 0; k < n - 1; ++k) {
1539 for (int i = k + 1; i < n; ++i) {
1540 if (abs(mat[k][k]) < TINY) {canSolve = false; break;}
1541 double ratio = mat[i][k] / mat[k][k];
1542 vec[i] -= ratio * vec[k];
1543 for (int j = k; j < n; ++j) mat[i][j] -= ratio * mat[k][j];
1545 if (!canSolve) break;
1548 for (int k = n - 1; k >= 0; --k) {
1549 for (int j = k + 1; j < n; ++j) vec[k] -= mat[k][j] * coefTmp[j];
1550 coefTmp[k] = vec[k] / mat[k][k];
1555 // Share evenly if failure.
1556 if (!canSolve) for (int i = 0; i < n; ++i) {
1559 if (vecSum > TINY) vecNor[i] = max(0.1, vec[i] / vecSum);
1562 // Normalize coefficients, with piece shared democratically.
1563 double coefSum = 0.;
1565 for (int i = 0; i < n; ++i) {
1566 coefTmp[i] = max( 0., coefTmp[i]);
1567 coefSum += coefTmp[i];
1568 vecSum += vecNor[i];
1570 if (coefSum > 0.) for (int i = 0; i < n; ++i) coef[i] = EVENFRAC / n
1571 + (1. - EVENFRAC) * 0.5 * (coefTmp[i] / coefSum + vecNor[i] / vecSum);
1572 else for (int i = 0; i < n; ++i) coef[i] = 1. / n;
1574 // Optional printout.
1576 os << " Solution: ";
1577 for (int i = 0; i < n; ++i) os << setw(12) << coef[i];
1584 // Setup mass selection for one resonance at a time - part 1.
1586 void PhaseSpace::setupMass1(int iM) {
1588 // Identity for mass seletion; is 0 also for light quarks (not yet selected).
1589 if (iM == 3) idMass[iM] = abs(sigmaProcessPtr->id3Mass());
1590 if (iM == 4) idMass[iM] = abs(sigmaProcessPtr->id4Mass());
1591 if (iM == 5) idMass[iM] = abs(sigmaProcessPtr->id5Mass());
1593 // Masses and widths of resonances.
1594 if (idMass[iM] == 0) {
1600 mPeak[iM] = ParticleDataTable::m0(idMass[iM]);
1601 mWidth[iM] = ParticleDataTable::mWidth(idMass[iM]);
1602 mMin[iM] = ParticleDataTable::mMin(idMass[iM]);
1603 mMax[iM] = ParticleDataTable::mMax(idMass[iM]);
1604 // gmZmode == 1 means pure photon propagator; set at lower mass limit.
1605 if (idMass[iM] == 23 && gmZmode == 1) mPeak[iM] = mMin[iM];
1608 // Mass and width combinations for Breit-Wigners.
1609 sPeak[iM] = mPeak[iM] * mPeak[iM];
1610 useBW[iM] = useBreitWigners && (mWidth[iM] > minWidthBreitWigners);
1611 if (!useBW[iM]) mWidth[iM] = 0.;
1612 mw[iM] = mPeak[iM] * mWidth[iM];
1613 wmRat[iM] = (idMass[iM] == 0 || mPeak[iM] == 0.)
1614 ? 0. : mWidth[iM] / mPeak[iM];
1616 // Simple Breit-Wigner range, upper edge to be corrected subsequently.
1618 mLower[iM] = mMin[iM];
1619 mUpper[iM] = mHatMax;
1626 // Setup mass selection for one resonance at a time - part 2.
1628 void PhaseSpace::setupMass2(int iM, double distToThresh) {
1630 // Store reduced Breit-Wigner range.
1631 if (mMax[iM] > mMin[iM]) mUpper[iM] = min( mUpper[iM], mMax[iM]);
1632 sLower[iM] = mLower[iM] * mLower[iM];
1633 sUpper[iM] = mUpper[iM] * mUpper[iM];
1635 // Prepare to select m3 by BW + flat + 1/s_3.
1636 // Determine relative coefficients by allowed mass range.
1637 if (distToThresh > THRESHOLDSIZE) {
1640 } else if (distToThresh > - THRESHOLDSIZE) {
1641 fracFlat[iM] = 0.25 - 0.15 * distToThresh / THRESHOLDSIZE;
1642 fracInv [iM] = 0.15 - 0.05 * distToThresh / THRESHOLDSIZE;
1648 // For gamma*/Z0: increase 1/s_i part and introduce 1/s_i^2 part.
1650 if (idMass[iM] == 23 && gmZmode == 0) {
1651 fracFlat[iM] *= 0.5;
1652 fracInv[iM] = 0.5 * fracInv[iM] + 0.25;
1653 fracInv2[iM] = 0.25;
1654 } else if (idMass[iM] == 23 && gmZmode == 1) {
1660 // Normalization integrals for the respective contribution.
1661 atanLower[iM] = atan( (sLower[iM] - sPeak[iM])/ mw[iM] );
1662 atanUpper[iM] = atan( (sUpper[iM] - sPeak[iM])/ mw[iM] );
1663 intBW[iM] = atanUpper[iM] - atanLower[iM];
1664 intFlat[iM] = sUpper[iM] - sLower[iM];
1665 intInv[iM] = log( sUpper[iM] / sLower[iM] );
1666 intInv2[iM] = 1./sLower[iM] - 1./sUpper[iM];
1672 // Select Breit-Wigner-distributed or fixed masses.
1674 void PhaseSpace::trialMass(int iM) {
1676 // References to masses to be set.
1677 double& mSet = (iM == 3) ? m3 : ( (iM == 4) ? m4 : m5 );
1678 double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 );
1680 // Distribution for m_i is BW + flat + 1/s_i + 1/s_i^2.
1682 double pickForm = Rndm::flat();
1683 if (pickForm > fracFlat[iM] + fracInv[iM] + fracInv2[iM])
1684 sSet = sPeak[iM] + mw[iM] * tan( atanLower[iM]
1685 + Rndm::flat() * intBW[iM] );
1686 else if (pickForm > fracInv[iM] + fracInv2[iM])
1687 sSet = sLower[iM] + Rndm::flat() * (sUpper[iM] - sLower[iM]);
1688 else if (pickForm > fracInv2[iM])
1689 sSet = sLower[iM] * pow( sUpper[iM] / sLower[iM], Rndm::flat() );
1690 else sSet = sLower[iM] * sUpper[iM]
1691 / (sLower[iM] + Rndm::flat() * (sUpper[iM] - sLower[iM]));
1694 // Else m_i is fixed at peak value.
1704 // Naively a fixed-width Breit-Wigner is used to pick the mass.
1705 // Here come the correction factors for
1706 // (i) preselection according to BW + flat in s_i + 1/s_i + 1/s_i^2,
1707 // (ii) reduced allowed mass range,
1708 // (iii) running width, i.e. m0*Gamma0 -> s*Gamma0/m0.
1709 // In the end, the weighted distribution is a running-width BW.
1711 double PhaseSpace::weightMass(int iM) {
1713 // Reference to mass and to Breit-Wigner weight to be set.
1714 double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 );
1715 double& runBWH = (iM == 3) ? runBW3H : ( (iM == 4) ? runBW4H : runBW5H );
1717 // Default weight if no Breit-Wigner.
1719 if (!useBW[iM]) return 1.;
1721 // Weight of generated distribution.
1722 double genBW = (1. - fracFlat[iM] - fracInv[iM] - fracInv2[iM])
1723 * mw[iM] / ( (pow2(sSet - sPeak[iM]) + pow2(mw[iM])) * intBW[iM])
1724 + fracFlat[iM] / intFlat[iM] + fracInv[iM] / (sSet * intInv[iM])
1725 + fracInv2[iM] / (sSet*sSet * intInv2[iM]);
1727 // Weight of distribution with running width in Breit-Wigner.
1728 double mwRun = sSet * wmRat[iM];
1729 runBWH = mwRun / (pow2(sSet - sPeak[iM]) + pow2(mwRun)) / M_PI;
1732 return (runBWH / genBW);
1736 //**************************************************************************
1738 // PhaseSpace2to1tauy class.
1739 // 2 -> 1 kinematics for normal subprocesses.
1743 // Set limits for resonance mass selection.
1745 bool PhaseSpace2to1tauy::setupMass() {
1747 // Treat Z0 as such or as gamma*/Z0
1748 gmZmode = gmZmodeGlobal;
1749 int gmZmodeProc = sigmaProcessPtr->gmZmode();
1750 if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
1752 // Mass limits for current resonance.
1753 int idRes = abs(sigmaProcessPtr->resonanceA());
1754 int idTmp = abs(sigmaProcessPtr->resonanceB());
1755 if (idTmp > 0) idRes = idTmp;
1756 double mResMin = (idRes == 0) ? 0. : ParticleDataTable::mMin(idRes);
1757 double mResMax = (idRes == 0) ? 0. : ParticleDataTable::mMax(idRes);
1759 // Compare with global mass limits and pick tighter of them.
1760 mHatMin = max( mResMin, mHatGlobalMin);
1761 sHatMin = mHatMin*mHatMin;
1763 if (mResMax > mResMin) mHatMax = min( mHatMax, mResMax);
1764 if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( mHatMax, mHatGlobalMax);
1765 sHatMax = mHatMax*mHatMax;
1767 // Default Breit-Wigner weight.
1770 // Fail if mass window (almost) closed.
1771 return (mHatMax > mHatMin + MASSMARGIN);
1777 // Construct the four-vector kinematics from the trial values.
1779 bool PhaseSpace2to1tauy::finalKin() {
1781 // Particle masses; incoming always on mass shell.
1786 // Incoming partons along beam axes. Outgoing has sum of momenta.
1787 pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
1788 pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
1789 pH[3] = pH[1] + pH[2];
1795 //**************************************************************************
1797 // PhaseSpace2to2tauyz class.
1798 // 2 -> 2 kinematics for normal subprocesses.
1802 // Set up for fixed or Breit-Wigner mass selection.
1804 bool PhaseSpace2to2tauyz::setupMasses() {
1806 // Treat Z0 as such or as gamma*/Z0
1807 gmZmode = gmZmodeGlobal;
1808 int gmZmodeProc = sigmaProcessPtr->gmZmode();
1809 if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
1811 // Set sHat limits - based on global limits only.
1812 mHatMin = mHatGlobalMin;
1813 sHatMin = mHatMin*mHatMin;
1815 if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax);
1816 sHatMax = mHatMax*mHatMax;
1818 // Masses and widths of resonances.
1822 // Reduced mass range when two massive particles.
1823 if (useBW[3]) mUpper[3] -= (useBW[4]) ? mMin[4] : mPeak[4];
1824 if (useBW[4]) mUpper[4] -= (useBW[3]) ? mMin[3] : mPeak[3];
1826 // If closed phase space then unallowed process.
1827 bool physical = true;
1828 if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false;
1829 if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false;
1830 if (!useBW[3] && !useBW[4] && mHatMax < mPeak[3] + mPeak[4] + MASSMARGIN)
1832 if (!physical) return false;
1834 // If either particle is massless then need extra pTHat cut.
1835 pTHatMin = pTHatGlobalMin;
1836 if (mPeak[3] < pTHatMinDiverge || mPeak[4] < pTHatMinDiverge)
1837 pTHatMin = max( pTHatMin, pTHatMinDiverge);
1838 pT2HatMin = pTHatMin * pTHatMin;
1839 pTHatMax = pTHatGlobalMax;
1840 pT2HatMax = pTHatMax * pTHatMax;
1842 // Prepare to select m3 by BW + flat + 1/s_3.
1844 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[3]
1845 / (pow2(mWidth[3]) + pow2(mWidth[4]));
1846 double distToThreshB = (mHatMax - mPeak[3] - mMin[4]) / mWidth[3];
1847 double distToThresh = min( distToThreshA, distToThreshB);
1848 setupMass2(3, distToThresh);
1851 // Prepare to select m4 by BW + flat + 1/s_4.
1853 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[4]
1854 / (pow2(mWidth[3]) + pow2(mWidth[4]));
1855 double distToThreshB = (mHatMax - mMin[3] - mPeak[4]) / mWidth[4];
1856 double distToThresh = min( distToThreshA, distToThreshB);
1857 setupMass2(4, distToThresh);
1860 // Initialization masses. Special cases when constrained phase space.
1861 m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3];
1862 m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4];
1863 if (m3 + m4 + THRESHOLDSIZE * (mWidth[3] + mWidth[4]) + MASSMARGIN
1865 if (useBW[3] && useBW[4]) physical = constrainedM3M4();
1866 else if (useBW[3]) physical = constrainedM3();
1867 else if (useBW[4]) physical = constrainedM4();
1872 // Correct selected mass-spectrum to running-width Breit-Wigner.
1873 // Extra safety margin for maximum search.
1875 if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX;
1876 if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX;
1886 // Select Breit-Wigner-distributed or fixed masses.
1888 bool PhaseSpace2to2tauyz::trialMasses() {
1890 // By default vanishing cross section.
1894 // Pick m3 and m4 independently.
1898 // If outside phase space then reject event.
1899 if (m3 + m4 + MASSMARGIN > mHatMax) return false;
1901 // Correct selected mass-spectrum to running-width Breit-Wigner.
1902 if (useBW[3]) wtBW *= weightMass(3);
1903 if (useBW[4]) wtBW *= weightMass(4);
1911 // Construct the four-vector kinematics from the trial values.
1913 bool PhaseSpace2to2tauyz::finalKin() {
1915 // Assign masses to particles assumed massless in matrix elements.
1916 int id3 = sigmaProcessPtr->id(3);
1917 int id4 = sigmaProcessPtr->id(4);
1918 if (idMass[3] == 0) { m3 = ParticleDataTable::m0(id3); s3 = m3*m3; }
1919 if (idMass[4] == 0) { m4 = ParticleDataTable::m0(id4); s4 = m4*m4; }
1921 // Sometimes swap tHat <-> uHat to reflect chosen final-state order.
1922 if (sigmaProcessPtr->swappedTU()) {
1927 // Check that phase space still open after new mass assignment.
1928 if (m3 + m4 + MASSMARGIN > mHat) {
1929 infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::finalKin: "
1930 "failed after mass assignment");
1933 p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1934 pAbs = sqrtpos( p2Abs );
1936 // Particle masses; incoming always on mass shell.
1942 // Incoming partons along beam axes.
1943 pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
1944 pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
1946 // Outgoing partons initially in collision CM frame along beam axes.
1947 pH[3] = Vec4( 0., 0., pAbs, 0.5 * (sH + s3 - s4) / mHat);
1948 pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (sH + s4 - s3) / mHat);
1950 // Then rotate and boost them to overall CM frame
1952 phi = 2. * M_PI * Rndm::flat();
1953 betaZ = (x1H - x2H)/(x1H + x2H);
1954 pH[3].rot( theta, phi);
1955 pH[4].rot( theta, phi);
1956 pH[3].bst( 0., 0., betaZ);
1957 pH[4].bst( 0., 0., betaZ);
1958 pTH = pAbs * sin(theta);
1966 // Special choice of m3 and m4 when mHatMax push them off mass shell.
1967 // Vary x in expression m3 + m4 = mHatMax - x * (Gamma3 + Gamma4).
1968 // For each x try to put either 3 or 4 as close to mass shell as possible.
1969 // Maximize BW_3 * BW_4 * beta_34, where latter approximate phase space.
1971 bool PhaseSpace2to2tauyz::constrainedM3M4() {
1974 bool foundNonZero = false;
1975 double wtMassMax = 0.;
1976 double m3WtMax = 0.;
1977 double m4WtMax = 0.;
1978 double xMax = (mHatMax - mLower[3] - mLower[4]) / (mWidth[3] + mWidth[4]);
1979 double xStep = THRESHOLDSTEP * min(1., xMax);
1981 double wtMassXbin, wtMassMaxOld, m34, mT34Min, wtMassNow,
1982 wtBW3Now, wtBW4Now, beta34Now;
1984 // Step through increasing x values.
1988 wtMassMaxOld = wtMassMax;
1989 m34 = mHatMax - xNow * (mWidth[3] + mWidth[4]);
1991 // Study point where m3 as close as possible to on-shell.
1992 m3 = min( mUpper[3], m34 - mLower[4]);
1993 if (m3 > mPeak[3]) m3 = max( mLower[3], mPeak[3]);
1995 if (m4 < mLower[4]) {m4 = mLower[4]; m3 = m34 - m4;}
1997 // Check that inside phase space limit set by pTmin.
1998 mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin);
1999 if (mT34Min < mHatMax) {
2001 // Breit-Wigners and beta factor give total weight.
2003 if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4]
2004 && m4 < mUpper[4]) {
2005 wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2006 wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2007 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2008 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2009 wtMassNow = wtBW3Now * wtBW4Now * beta34Now;
2012 // Store new maximum, if any.
2013 if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow;
2014 if (wtMassNow > wtMassMax) {
2015 foundNonZero = true;
2016 wtMassMax = wtMassNow;
2022 // Study point where m4 as close as possible to on-shell.
2023 m4 = min( mUpper[4], m34 - mLower[3]);
2024 if (m4 > mPeak[4]) m4 = max( mLower[4], mPeak[4]);
2026 if (m3 < mLower[3]) {m3 = mLower[3]; m4 = m34 - m3;}
2028 // Check that inside phase space limit set by pTmin.
2029 mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin);
2030 if (mT34Min < mHatMax) {
2032 // Breit-Wigners and beta factor give total weight.
2034 if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4]
2035 && m4 < mUpper[4]) {
2036 wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2037 wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2038 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2039 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2040 wtMassNow = wtBW3Now * wtBW4Now * beta34Now;
2043 // Store new maximum, if any.
2044 if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow;
2045 if (wtMassNow > wtMassMax) {
2046 foundNonZero = true;
2047 wtMassMax = wtMassNow;
2053 // Continue stepping if increasing trend and more x range available.
2054 } while ( (!foundNonZero || wtMassXbin > wtMassMaxOld)
2055 && xNow < xMax - xStep);
2057 // Restore best values for subsequent maximization. Return.
2060 return foundNonZero;
2066 // Special choice of m3 when mHatMax pushes it off mass shell.
2067 // Vary x in expression m3 = mHatMax - m4 - x * Gamma3.
2068 // Maximize BW_3 * beta_34, where latter approximate phase space.
2070 bool PhaseSpace2to2tauyz::constrainedM3() {
2073 bool foundNonZero = false;
2074 double wtMassMax = 0.;
2075 double m3WtMax = 0.;
2076 double mT4Min = sqrt(m4*m4 + pT2HatMin);
2077 double xMax = (mHatMax - mLower[3] - m4) / mWidth[3];
2078 double xStep = THRESHOLDSTEP * min(1., xMax);
2080 double wtMassNow, mT34Min, wtBW3Now, beta34Now;
2082 // Step through increasing x values; gives m3 unambiguously.
2086 m3 = mHatMax - m4 - xNow * mWidth[3];
2088 // Check that inside phase space limit set by pTmin.
2089 mT34Min = sqrt(m3*m3 + pT2HatMin) + mT4Min;
2090 if (mT34Min < mHatMax) {
2092 // Breit-Wigner and beta factor give total weight.
2093 wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2094 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2095 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2096 wtMassNow = wtBW3Now * beta34Now;
2098 // Store new maximum, if any.
2099 if (wtMassNow > wtMassMax) {
2100 foundNonZero = true;
2101 wtMassMax = wtMassNow;
2106 // Continue stepping if increasing trend and more x range available.
2107 } while ( (!foundNonZero || wtMassNow > wtMassMax)
2108 && xNow < xMax - xStep);
2110 // Restore best value for subsequent maximization. Return.
2112 return foundNonZero;
2118 // Special choice of m4 when mHatMax pushes it off mass shell.
2119 // Vary x in expression m4 = mHatMax - m3 - x * Gamma4.
2120 // Maximize BW_4 * beta_34, where latter approximate phase space.
2122 bool PhaseSpace2to2tauyz::constrainedM4() {
2125 bool foundNonZero = false;
2126 double wtMassMax = 0.;
2127 double m4WtMax = 0.;
2128 double mT3Min = sqrt(m3*m3 + pT2HatMin);
2129 double xMax = (mHatMax - mLower[4] - m3) / mWidth[4];
2130 double xStep = THRESHOLDSTEP * min(1., xMax);
2132 double wtMassNow, mT34Min, wtBW4Now, beta34Now;
2134 // Step through increasing x values; gives m4 unambiguously.
2138 m4 = mHatMax - m3 - xNow * mWidth[4];
2140 // Check that inside phase space limit set by pTmin.
2141 mT34Min = mT3Min + sqrt(m4*m4 + pT2HatMin);
2142 if (mT34Min < mHatMax) {
2144 // Breit-Wigner and beta factor give total weight.
2145 wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2146 beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2147 - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2148 wtMassNow = wtBW4Now * beta34Now;
2150 // Store new maximum, if any.
2151 if (wtMassNow > wtMassMax) {
2152 foundNonZero = true;
2153 wtMassMax = wtMassNow;
2158 // Continue stepping if increasing trend and more x range available.
2159 } while ( (!foundNonZero || wtMassNow > wtMassMax)
2160 && xNow < xMax - xStep);
2162 // Restore best value for subsequent maximization.
2164 return foundNonZero;
2168 //**************************************************************************
2170 // PhaseSpace2to2elastic class.
2171 // 2 -> 2 kinematics set up for elastic scattering.
2175 // Constants: could be changed here if desired, but normally should not.
2176 // These are of technical nature, as described for each.
2178 // Maximum positive/negative argument for exponentiation.
2179 const double PhaseSpace2to2elastic::EXPMAX = 50.;
2181 // Conversion coefficients = 1/(16pi) * (mb <-> GeV^2).
2182 const double PhaseSpace2to2elastic::CONVERTEL = 0.0510925;
2186 // Form of phase space sampling already fixed, so no optimization.
2187 // However, need to read out relevant parameters from SigmaTotal.
2189 bool PhaseSpace2to2elastic::setupSampling() {
2191 // Find maximum = value of cross section.
2192 sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2195 // Squared masses of particles.
2200 bSlope = sigmaTotPtr->bSlopeEl();
2202 // Determine maximum possible t range.
2203 lambda12S = pow2(s - s1 - s2) - 4. * s1 * s2 ;
2204 tLow = - lambda12S / s;
2207 // Production model with Coulomb corrections need more parameters.
2208 useCoulomb = Settings::flag("SigmaTotal:setOwn")
2209 && Settings::flag("SigmaElastic:setOwn");
2211 sigmaTot = sigmaTotPtr->sigmaTot();
2212 rho = Settings::parm("SigmaElastic:rho");
2213 lambda = Settings::parm("SigmaElastic:lambda");
2214 tAbsMin = Settings::parm("SigmaElastic:tAbsMin");
2215 phaseCst = Settings::parm("SigmaElastic:phaseConst");
2216 alphaEM0 = Settings::parm("StandardModel:alphaEM0");
2218 // Relative rate of nuclear and Coulombic parts in trials.
2220 sigmaNuc = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho) / bSlope
2221 * exp(-bSlope * tAbsMin);
2222 sigmaCou = (useCoulomb) ?
2223 pow2(alphaEM0) / (4. * CONVERTEL * tAbsMin) : 0.;
2224 signCou = (idA == idB) ? 1. : -1.;
2232 // Calculate coefficient of generation.
2233 tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2235 // Initialize alphaEM generation.
2236 alphaEM.init( Settings::mode("SigmaProcess:alphaEMorder") );
2244 // Select a trial kinematics phase space point. Perform full
2245 // Monte Carlo acceptance/rejection at this stage.
2247 bool PhaseSpace2to2elastic::trialKin( bool, bool ) {
2249 // Select t according to exp(bSlope*t).
2250 if (!useCoulomb || sigmaNuc > Rndm::flat() * (sigmaNuc + sigmaCou))
2251 tH = tUpp + log(1. + tAux * Rndm::flat()) / bSlope;
2253 // Select t according to 1/t^2.
2254 else tH = tLow * tUpp / (tUpp + Rndm::flat() * (tLow - tUpp));
2256 // Correction factor for ratio full/simulated.
2258 double sigmaN = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho)
2260 double alpEM = alphaEM.alphaEM(-tH);
2261 double sigmaC = pow2(alpEM) / (4. * CONVERTEL * tH*tH);
2262 double sigmaGen = 2. * (sigmaN + sigmaC);
2263 double form2 = pow4(lambda/(lambda - tH));
2264 double phase = signCou * alpEM
2265 * (-phaseCst - log(-0.5 * bSlope * tH));
2266 double sigmaCor = sigmaN + pow2(form2) * sigmaC
2267 - signCou * alpEM * sigmaTot * (form2 / (-tH))
2268 * exp(0.5 * bSlope * tH) * (rho * cos(phase) + sin(phase));
2269 sigmaNw = sigmaMx * sigmaCor / sigmaGen;
2272 // Careful reconstruction of scattering angle.
2273 double tRat = s * tH / lambda12S;
2274 double cosTheta = min(1., max(-1., 1. + 2. * tRat ) );
2275 double sinTheta = 2. * sqrtpos( -tRat * (1. + tRat) );
2276 theta = asin( min(1., sinTheta));
2277 if (cosTheta < 0.) theta = M_PI - theta;
2285 // Construct the four-vector kinematics from the trial values.
2287 bool PhaseSpace2to2elastic::finalKin() {
2295 // Incoming particles along beam axes.
2296 pAbs = 0.5 * sqrtpos(lambda12S) / eCM;
2297 pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2298 pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2300 // Outgoing particles initially along beam axes.
2301 pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2302 pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2305 phi = 2. * M_PI * Rndm::flat();
2306 pH[3].rot( theta, phi);
2307 pH[4].rot( theta, phi);
2309 // Set some further info for completeness.
2313 uH = 2. * (s1 + s2) - sH - tH;
2315 p2Abs = pAbs * pAbs;
2317 pTH = pAbs * sin(theta);
2324 //**************************************************************************
2326 // PhaseSpace2to2diffractive class.
2327 // 2 -> 2 kinematics set up for diffractive scattering.
2331 // Constants: could be changed here if desired, but normally should not.
2332 // These are of technical nature, as described for each.
2334 // Number of tries to find acceptable (m^2, t) set.
2335 const int PhaseSpace2to2diffractive::NTRY = 500;
2337 // Maximum positive/negative argument for exponentiation.
2338 const double PhaseSpace2to2diffractive::EXPMAX = 50.;
2340 // Safety margin so sum of diffractive masses not too close to eCM.
2341 const double PhaseSpace2to2diffractive::DIFFMASSMAX = 1e-8;
2345 // Form of phase space sampling already fixed, so no optimization.
2346 // However, need to read out relevant parameters from SigmaTotal.
2348 bool PhaseSpace2to2diffractive::setupSampling() {
2350 // Find maximum = value of cross section.
2351 sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2354 // Masses of particles and minimal masses of diffractive states.
2355 m3ElDiff = (isDiffA) ? sigmaTotPtr->mMinXB() : mA;
2356 m4ElDiff = (isDiffB) ? sigmaTotPtr->mMinAX() : mB;
2359 s3 = pow2( m3ElDiff);
2360 s4 = pow2( m4ElDiff);
2362 // Parameters of low-mass-resonance diffractive enhancement.
2363 cRes = sigmaTotPtr->cRes();
2364 sResXB = pow2( sigmaTotPtr->mResXB());
2365 sResAX = pow2( sigmaTotPtr->mResAX());
2366 sProton = sigmaTotPtr->sProton();
2368 // Lower limit diffractive slope.
2369 if (!isDiffB) bMin = sigmaTotPtr->bMinSlopeXB();
2370 else if (!isDiffA) bMin = sigmaTotPtr->bMinSlopeAX();
2371 else bMin = sigmaTotPtr->bMinSlopeXX();
2373 // Determine maximum possible t range and coefficient of generation.
2374 lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
2375 lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2376 double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2377 double tempB = lambda12 * lambda34 / s;
2378 double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2379 * (s1 * s4 - s2 * s3) / s;
2380 tLow = -0.5 * (tempA + tempB);
2381 tUpp = tempC / tLow;
2382 tAux = exp( max(-EXPMAX, bMin * (tLow - tUpp)) ) - 1.;
2390 // Select a trial kinematics phase space point. Perform full
2391 // Monte Carlo acceptance/rejection at this stage.
2393 bool PhaseSpace2to2diffractive::trialKin( bool, bool ) {
2395 // Loop over attempts to set up masses and t consistently.
2396 for (int loop = 0; ; ++loop) {
2398 infoPtr->errorMsg("Error in PhaseSpace2to2diffractive::trialKin: "
2399 " quit after repeated tries");
2403 // Select diffractive mass/masses according to dm^2/m^2.
2404 m3 = (isDiffA) ? m3ElDiff * pow( max(mA, eCM - m4ElDiff) / m3ElDiff,
2405 Rndm::flat()) : m3ElDiff;
2406 m4 = (isDiffB) ? m4ElDiff * pow( max(mB, eCM - m3ElDiff) / m4ElDiff,
2407 Rndm::flat()) : m4ElDiff;
2411 // Additional mass factors, including resonance enhancement.
2412 if (m3 + m4 >= eCM) continue;
2413 if (isDiffA && !isDiffB) {
2414 double facXB = (1. - s3 / s)
2415 * (1. + cRes * sResXB / (sResXB + s3));
2416 if (facXB < Rndm::flat() * (1. + cRes)) continue;
2417 } else if (isDiffB && !isDiffA) {
2418 double facAX = (1. - s4 / s)
2419 * (1. + cRes * sResAX / (sResAX + s4));
2420 if (facAX < Rndm::flat() * (1. + cRes)) continue;
2422 double facXX = (1. - pow2(m3 + m4) / s)
2423 * (s * sProton / (s * sProton + s3 * s4))
2424 * (1. + cRes * sResXB / (sResXB + s3))
2425 * (1. + cRes * sResAX / (sResAX + s4));
2426 if (facXX < Rndm::flat() * pow2(1. + cRes)) continue;
2429 // Select t according to exp(bMin*t) and correct to right slope.
2430 tH = tUpp + log(1. + tAux * Rndm::flat()) / bMin;
2432 if (isDiffA && !isDiffB) bDiff = sigmaTotPtr->bSlopeXB(s3) - bMin;
2433 else if (!isDiffA) bDiff = sigmaTotPtr->bSlopeAX(s4) - bMin;
2434 else bDiff = sigmaTotPtr->bSlopeXX(s3, s4) - bMin;
2435 bDiff = max(0., bDiff);
2436 if (exp( max(-50., bDiff * (tH - tUpp)) ) < Rndm::flat()) continue;
2438 // Check whether m^2 and t choices are consistent.
2439 lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2440 double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2441 double tempB = lambda12 * lambda34 / s;
2442 if (tempB < DIFFMASSMAX) continue;
2443 double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2444 * (s1 * s4 - s2 * s3) / s;
2445 double tLowNow = -0.5 * (tempA + tempB);
2446 double tUppNow = tempC / tLowNow;
2447 if (tH < tLowNow || tH > tUppNow) continue;
2449 // Careful reconstruction of scattering angle.
2450 double cosTheta = min(1., max(-1., (tempA + 2. * tH) / tempB));
2451 double sinTheta = 2. * sqrtpos( -(tempC + tempA * tH + tH * tH) )
2453 theta = asin( min(1., sinTheta));
2454 if (cosTheta < 0.) theta = M_PI - theta;
2456 // Found acceptable kinematics, so no more looping.
2466 // Construct the four-vector kinematics from the trial values.
2468 bool PhaseSpace2to2diffractive::finalKin() {
2470 // Particle masses; incoming always on mass shell.
2476 // Incoming particles along beam axes.
2477 pAbs = 0.5 * lambda12 / eCM;
2478 pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2479 pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2481 // Outgoing particles initially along beam axes.
2482 pAbs = 0.5 * lambda34 / eCM;
2483 pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s3 - s4) / eCM);
2484 pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s4 - s3) / eCM);
2487 phi = 2. * M_PI * Rndm::flat();
2488 pH[3].rot( theta, phi);
2489 pH[4].rot( theta, phi);
2491 // Set some further info for completeness.
2495 uH = s1 + s2 + s3 + s4 - sH - tH;
2497 p2Abs = pAbs * pAbs;
2499 pTH = pAbs * sin(theta);
2506 //**************************************************************************
2508 // PhaseSpace2to3tauycyl class.
2509 // 2 -> 3 kinematics for normal subprocesses.
2513 // Constants: could be changed here if desired, but normally should not.
2514 // These are of technical nature, as described for each.
2516 // Number of Newton-Raphson iterations of kinematics when masses introduced.
2517 const int PhaseSpace2to3tauycyl::NITERNR = 5;
2521 // Set up for fixed or Breit-Wigner mass selection.
2523 bool PhaseSpace2to3tauycyl::setupMasses() {
2525 // Treat Z0 as such or as gamma*/Z0
2526 gmZmode = gmZmodeGlobal;
2527 int gmZmodeProc = sigmaProcessPtr->gmZmode();
2528 if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
2530 // Set sHat limits - based on global limits only.
2531 mHatMin = mHatGlobalMin;
2532 sHatMin = mHatMin*mHatMin;
2534 if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax);
2535 sHatMax = mHatMax*mHatMax;
2537 // Masses and widths of resonances.
2542 // Reduced mass range - do not make it as fancy as in two-body case.
2543 if (useBW[3]) mUpper[3] -= (mPeak[4] + mPeak[5]);
2544 if (useBW[4]) mUpper[4] -= (mPeak[3] + mPeak[5]);
2545 if (useBW[5]) mUpper[5] -= (mPeak[3] + mPeak[4]);
2547 // If closed phase space then unallowed process.
2548 bool physical = true;
2549 if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false;
2550 if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false;
2551 if (useBW[5] && mUpper[5] < mLower[5] + MASSMARGIN) physical = false;
2552 if (!useBW[3] && !useBW[4] && !useBW[5] && mHatMax < mPeak[3]
2553 + mPeak[4] + mPeak[5] + MASSMARGIN) physical = false;
2554 if (!physical) return false;
2556 // No extra pT precautions in massless limit - assumed fixed by ME's.
2557 pTHatMin = pTHatGlobalMin;
2558 pT2HatMin = pTHatMin * pTHatMin;
2559 pTHatMax = pTHatGlobalMax;
2560 pT2HatMax = pTHatMax * pTHatMax;
2562 // Prepare to select m3 by BW + flat + 1/s_3.
2564 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
2565 * mWidth[3] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
2566 double distToThreshB = (mHatMax - mPeak[3] - mMin[4] - mMin[5])
2568 double distToThresh = min( distToThreshA, distToThreshB);
2569 setupMass2(3, distToThresh);
2572 // Prepare to select m4 by BW + flat + 1/s_3.
2574 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
2575 * mWidth[4] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
2576 double distToThreshB = (mHatMax - mPeak[4] - mMin[3] - mMin[5])
2578 double distToThresh = min( distToThreshA, distToThreshB);
2579 setupMass2(4, distToThresh);
2582 // Prepare to select m5 by BW + flat + 1/s_3.
2584 double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
2585 * mWidth[5] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
2586 double distToThreshB = (mHatMax - mPeak[5] - mMin[3] - mMin[4])
2588 double distToThresh = min( distToThreshA, distToThreshB);
2589 setupMass2(5, distToThresh);
2592 // Initialization masses. For now give up when constrained phase space.
2593 m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3];
2594 m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4];
2595 m5 = (useBW[5]) ? min(mPeak[5], mUpper[5]) : mPeak[5];
2596 if (m3 + m4 + m5 + MASSMARGIN > mHatMax) physical = false;
2601 // Correct selected mass-spectrum to running-width Breit-Wigner.
2602 // Extra safety margin for maximum search.
2604 if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX;
2605 if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX;
2606 if (useBW[5]) wtBW *= weightMass(5) * EXTRABWWTMAX;
2615 // Select Breit-Wigner-distributed or fixed masses.
2617 bool PhaseSpace2to3tauycyl::trialMasses() {
2619 // By default vanishing cross section.
2623 // Pick m3, m4 and m5 independently.
2628 // If outside phase space then reject event.
2629 if (m3 + m4 + m5 + MASSMARGIN > mHatMax) return false;
2631 // Correct selected mass-spectrum to running-width Breit-Wigner.
2632 if (useBW[3]) wtBW *= weightMass(3);
2633 if (useBW[4]) wtBW *= weightMass(4);
2634 if (useBW[5]) wtBW *= weightMass(5);
2643 // Construct the four-vector kinematics from the trial values.
2645 bool PhaseSpace2to3tauycyl::finalKin() {
2647 // Assign masses to particles assumed massless in matrix elements.
2648 int id3 = sigmaProcessPtr->id(3);
2649 int id4 = sigmaProcessPtr->id(4);
2650 int id5 = sigmaProcessPtr->id(5);
2651 if (idMass[3] == 0) { m3 = ParticleDataTable::m0(id3); s3 = m3*m3; }
2652 if (idMass[4] == 0) { m4 = ParticleDataTable::m0(id4); s4 = m4*m4; }
2653 if (idMass[5] == 0) { m5 = ParticleDataTable::m0(id5); s5 = m5*m5; }
2655 // Check that phase space still open after new mass assignment.
2656 if (m3 + m4 + m5 + MASSMARGIN > mHat) {
2657 infoPtr->errorMsg("Warning in PhaseSpace2to3tauycyl::finalKin: "
2658 "failed after mass assignment");
2662 // Particle masses; incoming always on mass shell.
2669 // Incoming partons along beam axes.
2670 pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
2671 pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
2673 // Begin three-momentum rescaling to compensate for masses.
2674 if (idMass[3] == 0 || idMass[4] == 0 || idMass[5] == 0) {
2675 double p3S = p3cm.pAbs2();
2676 double p4S = p4cm.pAbs2();
2677 double p5S = p5cm.pAbs2();
2679 double e3, e4, e5, value, deriv;
2681 // Iterate rescaling solution five times, using Newton-Raphson.
2682 for (int i = 0; i < NITERNR; ++i) {
2683 e3 = sqrt(s3 + fac * p3S);
2684 e4 = sqrt(s4 + fac * p4S);
2685 e5 = sqrt(s5 + fac * p5S);
2686 value = e3 + e4 + e5 - mHat;
2687 deriv = 0.5 * (p3S / e3 + p4S / e4 + p5S / e5);
2688 fac -= value / deriv;
2691 // Rescale momenta appropriately.
2692 double facRoot = sqrt(fac);
2693 p3cm.rescale3( facRoot );
2694 p4cm.rescale3( facRoot );
2695 p5cm.rescale3( facRoot );
2696 p3cm.e( sqrt(s3 + fac * p3S) );
2697 p4cm.e( sqrt(s4 + fac * p4S) );
2698 p5cm.e( sqrt(s5 + fac * p5S) );
2701 // Outgoing partons initially in collision CM frame along beam axes.
2706 // Then boost them to overall CM frame
2707 betaZ = (x1H - x2H)/(x1H + x2H);
2708 pH[3].rot( theta, phi);
2709 pH[4].rot( theta, phi);
2710 pH[3].bst( 0., 0., betaZ);
2711 pH[4].bst( 0., 0., betaZ);
2712 pH[5].bst( 0., 0., betaZ);
2714 // Store average pT of three final particles for documentation.
2715 pTH = (p3cm.pT() + p4cm.pT() + p5cm.pT()) / 3.;
2721 //**************************************************************************
2723 // The PhaseSpaceLHA class.
2724 // A derived class for Les Houches events.
2725 // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
2729 // Constants: could be changed here if desired, but normally should not.
2730 // These are of technical nature, as described for each.
2732 // LHA convention with cross section in pb forces conversion to mb.
2733 const double PhaseSpaceLHA::CONVERTPB2MB = 1e-9;
2737 // Find maximal cross section for comparison with internal processes.
2739 bool PhaseSpaceLHA::setupSampling() {
2741 // Find which strategy Les Houches events are produced with.
2742 strategy = lhaUpPtr->strategy();
2743 stratAbs = abs(strategy);
2744 if (strategy == 0 || stratAbs > 4) {
2745 ostringstream stratCode;
2746 stratCode << strategy;
2747 infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: unknown "
2748 "Les Houches Accord weighting stategy", stratCode.str());
2752 // Number of contributing processes.
2753 nProc = lhaUpPtr->sizeProc();
2755 // Loop over all processes. Read out maximum and cross section.
2759 double xMax, xSec, xMaxAbs;
2760 for (int iProc = 0 ; iProc < nProc; ++iProc) {
2761 idPr = lhaUpPtr->idProcess(iProc);
2762 xMax = lhaUpPtr->xMax(iProc);
2763 xSec = lhaUpPtr->xSec(iProc);
2765 // Check for inconsistencies between strategy and stored values.
2766 if ( (strategy == 1 || strategy == 2) && xMax < 0.) {
2767 infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: "
2768 "negative maximum not allowed");
2771 if ( ( strategy == 2 || strategy == 3) && xSec < 0.) {
2772 infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: "
2773 "negative cross section not allowed");
2777 // Store maximal cross sections for later choice.
2778 if (stratAbs == 1) xMaxAbs = abs(xMax);
2779 else if (stratAbs < 4) xMaxAbs = abs(xSec);
2781 idProc.push_back( idPr );
2782 xMaxAbsProc.push_back( xMaxAbs );
2784 // Find sum and convert to mb.
2785 xMaxAbsSum += xMaxAbs;
2788 sigmaMx = xMaxAbsSum * CONVERTPB2MB;
2789 sigmaSgn = xSecSgnSum * CONVERTPB2MB;
2798 // Construct the next process, by interface to Les Houches class.
2800 bool PhaseSpaceLHA::trialKin( bool, bool repeatSame ) {
2802 // Must select process type in some cases.
2804 if (repeatSame) idProcNow = idProcSave;
2805 else if (stratAbs <= 2) {
2806 double xMaxAbsRndm = xMaxAbsSum * Rndm::flat();
2808 do xMaxAbsRndm -= xMaxAbsProc[++iProc];
2809 while (xMaxAbsRndm > 0. && iProc < nProc - 1);
2810 idProcNow = idProc[iProc];
2813 // Generate Les Houches event. Return if fail (= end of file).
2814 bool physical = lhaUpPtr->setEvent(idProcNow);
2815 if (!physical) return false;
2817 // Find which process was generated.
2818 int idPr = lhaUpPtr->idProcess();
2820 for (int iP = 0; iP < int(idProc.size()); ++iP)
2821 if (idProc[iP] == idPr) iProc = iP;
2824 // Extract cross section and rescale according to strategy.
2825 double wtPr = lhaUpPtr->weight();
2826 if (stratAbs == 1) sigmaNw = wtPr * CONVERTPB2MB
2827 * xMaxAbsSum / xMaxAbsProc[iProc];
2828 else if (stratAbs == 2) sigmaNw = (wtPr / abs(lhaUpPtr->xMax(iProc)))
2830 else if (strategy == 3) sigmaNw = sigmaMx;
2831 else if (strategy == -3 && wtPr > 0.) sigmaNw = sigmaMx;
2832 else if (strategy == -3) sigmaNw = -sigmaMx;
2833 else if (stratAbs == 4) sigmaNw = wtPr * CONVERTPB2MB;
2840 //**************************************************************************
2842 } // end namespace Pythia8