PYTHIA8/pythia8145/src/PhaseSpace.cxx

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