PYTHIA8/pythia8170/src/PhaseSpace.cxx

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