PYTHIA8/pythia8170/src/SigmaProcess.cxx

   1 // SigmaProcess.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2012 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the
   7 // SigmaProcess class, and classes derived from it.
   8
   9 #include "SigmaProcess.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The SigmaProcess class.
  16 // Base class for cross sections.
  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 // Conversion of GeV^{-2} to mb for cross section.
  24 const double SigmaProcess::CONVERT2MB    = 0.389380;
  25
  26 // The sum of outgoing masses must not be too close to the cm energy.
  27 const double SigmaProcess::MASSMARGIN    = 0.1;
  28
  29 // Parameters of momentum rescaling procedure: maximally allowed
  30 // relative energy error and number of iterations.
  31 const double SigmaProcess::COMPRELERR = 1e-10;
  32 const int    SigmaProcess::NCOMPSTEP  = 10;
  33
  34 //--------------------------------------------------------------------------
  35
  36 // Perform simple initialization and store pointers.
  37
  38 void SigmaProcess::init(Info* infoPtrIn, Settings* settingsPtrIn,
  39   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, BeamParticle* beamAPtrIn,
  40   BeamParticle* beamBPtrIn, Couplings* couplingsPtrIn,
  41   SigmaTotal* sigmaTotPtrIn, SusyLesHouches* slhaPtrIn) {
  42
  43   // Store pointers.
  44   infoPtr         = infoPtrIn;
  45   settingsPtr     = settingsPtrIn;
  46   particleDataPtr = particleDataPtrIn;
  47   rndmPtr         = rndmPtrIn;
  48   beamAPtr        = beamAPtrIn;
  49   beamBPtr        = beamBPtrIn;
  50   couplingsPtr    = couplingsPtrIn;
  51   sigmaTotPtr     = sigmaTotPtrIn;
  52   slhaPtr         = slhaPtrIn;
  53
  54   // Read out some properties of beams to allow shorthand.
  55   idA             = (beamAPtr != 0) ? beamAPtr->id() : 0;
  56   idB             = (beamBPtr != 0) ? beamBPtr->id() : 0;
  57   mA              = (beamAPtr != 0) ? beamAPtr->m() : 0.;
  58   mB              = (beamBPtr != 0) ? beamBPtr->m() : 0.;
  59   isLeptonA       = (beamAPtr != 0) ? beamAPtr->isLepton() : false;
  60   isLeptonB       = (beamBPtr != 0) ? beamBPtr->isLepton() : false;
  61   hasLeptonBeams  = isLeptonA || isLeptonB;
  62
  63   // K factor, multiplying resolved processes. (But not here for MPI.)
  64   Kfactor         = settingsPtr->parm("SigmaProcess:Kfactor");
  65
  66   // Maximum incoming quark flavour.
  67   nQuarkIn        = settingsPtr->mode("PDFinProcess:nQuarkIn");
  68
  69   // Medium heavy fermion masses set massless or not in ME expressions.
  70   mcME            = (settingsPtr->flag("SigmaProcess:cMassiveME"))
  71                   ? particleDataPtr->m0(4)  : 0.;
  72   mbME            = (settingsPtr->flag("SigmaProcess:bMassiveME"))
  73                   ? particleDataPtr->m0(5)  : 0.;
  74   mmuME           = (settingsPtr->flag("SigmaProcess:muMassiveME"))
  75                   ? particleDataPtr->m0(13) : 0.;
  76   mtauME          = (settingsPtr->flag("SigmaProcess:tauMassiveME"))
  77                   ? particleDataPtr->m0(15) : 0.;
  78
  79   // Renormalization scale choice.
  80   renormScale1    = settingsPtr->mode("SigmaProcess:renormScale1");
  81   renormScale2    = settingsPtr->mode("SigmaProcess:renormScale2");
  82   renormScale3    = settingsPtr->mode("SigmaProcess:renormScale3");
  83   renormScale3VV  = settingsPtr->mode("SigmaProcess:renormScale3VV");
  84   renormMultFac   = settingsPtr->parm("SigmaProcess:renormMultFac");
  85   renormFixScale  = settingsPtr->parm("SigmaProcess:renormFixScale");
  86
  87   // Factorization scale choice.
  88   factorScale1    = settingsPtr->mode("SigmaProcess:factorScale1");
  89   factorScale2    = settingsPtr->mode("SigmaProcess:factorScale2");
  90   factorScale3    = settingsPtr->mode("SigmaProcess:factorScale3");
  91   factorScale3VV  = settingsPtr->mode("SigmaProcess:factorScale3VV");
  92   factorMultFac   = settingsPtr->parm("SigmaProcess:factorMultFac");
  93   factorFixScale  = settingsPtr->parm("SigmaProcess:factorFixScale");
  94
  95   // CP violation parameters for the BSM Higgs sector.
  96   higgsH1parity   = settingsPtr->mode("HiggsH1:parity");
  97   higgsH1eta      = settingsPtr->parm("HiggsH1:etaParity");
  98   higgsH2parity   = settingsPtr->mode("HiggsH2:parity");
  99   higgsH2eta      = settingsPtr->parm("HiggsH2:etaParity");
 100   higgsA3parity   = settingsPtr->mode("HiggsA3:parity");
 101   higgsA3eta      = settingsPtr->parm("HiggsA3:etaParity");
 102
 103   // If BSM not switched on then H1 should have SM properties.
 104   if (!settingsPtr->flag("Higgs:useBSM")){
 105     higgsH1parity = 1;
 106     higgsH1eta    = 0.;
 107   }
 108
 109 }
 110
 111 //--------------------------------------------------------------------------
 112
 113 // Set up allowed flux of incoming partons.
 114 // addBeam: set up PDF's that need to be evaluated for the two beams.
 115 // addPair: set up pairs of incoming partons from the two beams.
 116
 117 bool SigmaProcess::initFlux() {
 118
 119   // Reset arrays (in case of several init's in same run).
 120   inBeamA.clear();
 121   inBeamB.clear();
 122   inPair.clear();
 123
 124   // Read in process-specific channel information.
 125   string fluxType = inFlux();
 126
 127   // Case with g g incoming state.
 128   if (fluxType == "gg") {
 129     addBeamA(21);
 130     addBeamB(21);
 131     addPair(21, 21);
 132   }
 133
 134   // Case with q g incoming state.
 135   else if (fluxType == "qg") {
 136     for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
 137       int idNow = (i == 0) ? 21 : i;
 138       addBeamA(idNow);
 139       addBeamB(idNow);
 140     }
 141     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 142     if (idNow != 0) {
 143       addPair(idNow, 21);
 144       addPair(21, idNow);
 145     }
 146   }
 147
 148   // Case with q q', q qbar' or qbar qbar' incoming state.
 149   else if (fluxType == "qq") {
 150     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 151     if (idNow != 0) {
 152       addBeamA(idNow);
 153       addBeamB(idNow);
 154     }
 155     for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 156     if (id1Now != 0)
 157     for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 158     if (id2Now != 0)
 159       addPair(id1Now, id2Now);
 160   }
 161
 162   // Case with q qbar incoming state.
 163   else if (fluxType == "qqbarSame") {
 164     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 165     if (idNow != 0) {
 166       addBeamA(idNow);
 167       addBeamB(idNow);
 168     }
 169     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 170     if (idNow != 0)
 171       addPair(idNow, -idNow);
 172   }
 173
 174   // Case with f f', f fbar', fbar fbar' incoming state.
 175   else if (fluxType == "ff") {
 176     // If beams are leptons then they are also the colliding partons.
 177     if ( isLeptonA && isLeptonB ) {
 178       addBeamA(idA);
 179       addBeamB(idB);
 180       addPair(idA, idB);
 181     // First beam is lepton and second is hadron.
 182     } else if ( isLeptonA ) {
 183       addBeamA(idA);
 184       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 185       if (idNow != 0) {
 186         addBeamB(idNow);
 187         addPair(idA, idNow);
 188       }
 189     // First beam is hadron and second is lepton.
 190     } else if ( isLeptonB ) {
 191       addBeamB(idB);
 192       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 193       if (idNow != 0) {
 194         addBeamA(idNow);
 195         addPair(idNow, idB);
 196       }
 197     // Hadron beams gives quarks.
 198     } else {
 199       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 200       if (idNow != 0) {
 201         addBeamA(idNow);
 202         addBeamB(idNow);
 203       }
 204       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 205       if (id1Now != 0)
 206       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 207       if (id2Now != 0)
 208         addPair(id1Now, id2Now);
 209     }
 210   }
 211
 212   // Case with f fbar incoming state.
 213   else if (fluxType == "ffbarSame") {
 214     // If beams are antiparticle pair and leptons then also colliding partons.
 215     if ( idA + idB == 0 && isLeptonA ) {
 216       addBeamA(idA);
 217       addBeamB(idB);
 218       addPair(idA, idB);
 219     // Else assume both to be hadrons, for better or worse.
 220     } else {
 221       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 222       if (idNow != 0) {
 223         addBeamA(idNow);
 224         addBeamB(idNow);
 225       }
 226       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 227       if (idNow != 0)
 228         addPair(idNow, -idNow);
 229     }
 230   }
 231
 232   // Case with f fbar' charged(+-1) incoming state.
 233   else if (fluxType == "ffbarChg") {
 234     // If beams are leptons then also colliding partons.
 235     if ( isLeptonA && isLeptonB && abs( particleDataPtr->chargeType(idA)
 236              + particleDataPtr->chargeType(idB) ) == 3 ) {
 237       addBeamA(idA);
 238       addBeamB(idB);
 239       addPair(idA, idB);
 240     // Hadron beams gives quarks.
 241     } else {
 242       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 243       if (idNow != 0) {
 244         addBeamA(idNow);
 245         addBeamB(idNow);
 246       }
 247       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 248       if (id1Now != 0)
 249       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 250       if (id2Now != 0 && id1Now * id2Now < 0
 251         && (abs(id1Now) + abs(id2Now))%2 == 1) addPair(id1Now, id2Now);
 252     }
 253   }
 254
 255   // Case with f fbar' generic incoming state.
 256   else if (fluxType == "ffbar") {
 257     // If beams are leptons then also colliding partons.
 258     if (isLeptonA && isLeptonB && idA * idB < 0) {
 259       addBeamA(idA);
 260       addBeamB(idB);
 261       addPair(idA, idB);
 262     // Hadron beams gives quarks.
 263     } else {
 264       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 265       if (idNow != 0) {
 266         addBeamA(idNow);
 267         addBeamB(idNow);
 268       }
 269       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 270       if (id1Now != 0)
 271       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 272       if (id2Now != 0 && id1Now * id2Now < 0)
 273         addPair(id1Now, id2Now);
 274     }
 275   }
 276
 277   // Case with f gamma incoming state.
 278   else if (fluxType == "fgm") {
 279     // Fermion from incoming side A.
 280     if ( isLeptonA ) {
 281       addBeamA(idA);
 282       addPair(idA, 22);
 283     } else {
 284       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 285       if (idNow != 0) {
 286         addBeamA(idNow);
 287         addPair(idNow, 22);
 288       }
 289     }
 290     // Fermion from incoming side B.
 291     if ( isLeptonB ) {
 292       addBeamB( idB);
 293       addPair(22, idB);
 294     } else {
 295       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 296       if (idNow != 0) {
 297         addBeamB(idNow);
 298         addPair(22, idNow);
 299       }
 300     }
 301     // Photons in the beams.
 302     addBeamA(22);
 303     addBeamB(22);
 304   }
 305
 306   // Case with gamma gamma incoming state.
 307   else if (fluxType == "ggm") {
 308     addBeamA(21);
 309     addBeamA(22);
 310     addBeamB(21);
 311     addBeamB(22);
 312     addPair(21, 22);
 313     addPair(22, 21);
 314   }
 315
 316   // Case with gamma gamma incoming state.
 317   else if (fluxType == "gmgm") {
 318     addBeamA(22);
 319     addBeamB(22);
 320     addPair(22, 22);
 321   }
 322
 323   // Unrecognized fluxType is bad sign. Else done.
 324   else {
 325     infoPtr->errorMsg("Error in SigmaProcess::initFlux: "
 326     "unrecognized inFlux type", fluxType);
 327     return false;
 328   }
 329   return true;
 330
 331 }
 332
 333 //--------------------------------------------------------------------------
 334
 335 // Convolute matrix-element expression(s) with parton flux and K factor.
 336
 337 double SigmaProcess::sigmaPDF() {
 338
 339   // Evaluate and store the required parton densities.
 340   for (int j = 0; j < sizeBeamA(); ++j)
 341     inBeamA[j].pdf = beamAPtr->xfHard( inBeamA[j].id, x1Save, Q2FacSave);
 342   for (int j = 0; j < sizeBeamB(); ++j)
 343     inBeamB[j].pdf = beamBPtr->xfHard( inBeamB[j].id, x2Save, Q2FacSave);
 344
 345   // Loop over allowed incoming channels.
 346   sigmaSumSave = 0.;
 347   for (int i = 0; i < sizePair(); ++i) {
 348
 349     // Evaluate hard-scattering cross section. Include K factor.
 350     inPair[i].pdfSigma = Kfactor
 351                        * sigmaHatWrap(inPair[i].idA, inPair[i].idB);
 352
 353     // Multiply by respective parton densities.
 354     for (int j = 0; j < sizeBeamA(); ++j)
 355     if (inPair[i].idA == inBeamA[j].id) {
 356       inPair[i].pdfA      = inBeamA[j].pdf;
 357       inPair[i].pdfSigma *= inBeamA[j].pdf;
 358       break;
 359     }
 360     for (int j = 0; j < sizeBeamB(); ++j)
 361     if (inPair[i].idB == inBeamB[j].id) {
 362       inPair[i].pdfB      = inBeamB[j].pdf;
 363       inPair[i].pdfSigma *= inBeamB[j].pdf;
 364       break;
 365     }
 366
 367     // Sum for all channels.
 368     sigmaSumSave += inPair[i].pdfSigma;
 369   }
 370
 371   // Done.
 372   return sigmaSumSave;
 373
 374 }
 375
 376 //--------------------------------------------------------------------------
 377
 378 // Select incoming parton channel and extract parton densities (resolved).
 379
 380 void SigmaProcess::pickInState(int id1in, int id2in) {
 381
 382   // Multiparton interactions: partons already selected.
 383   if (id1in != 0 && id2in != 0) {
 384     id1 = id1in;
 385     id2 = id2in;
 386   }
 387
 388   // Pick channel. Extract channel flavours and pdf's.
 389   double sigmaRand =  sigmaSumSave * rndmPtr->flat();
 390   for (int i = 0; i < sizePair(); ++i) {
 391     sigmaRand -= inPair[i].pdfSigma;
 392     if (sigmaRand <= 0.) {
 393       id1      = inPair[i].idA;
 394       id2      = inPair[i].idB;
 395       pdf1Save = inPair[i].pdfA;
 396       pdf2Save = inPair[i].pdfB;
 397       break;
 398     }
 399   }
 400
 401 }
 402
 403 //--------------------------------------------------------------------------
 404
 405 // Calculate incoming modified masses and four-vectors for matrix elements.
 406
 407 bool SigmaProcess::setupForMEin() {
 408
 409   // Initially assume it will work out to set up modified kinematics.
 410   bool allowME = true;
 411
 412   // Correct incoming c, b, mu and tau to be massive or not.
 413   mME[0] = 0.;
 414   int id1Tmp = abs(id1);
 415   if (id1Tmp ==  4) mME[0] = mcME;
 416   if (id1Tmp ==  5) mME[0] = mbME;
 417   if (id1Tmp == 13) mME[0] = mmuME;
 418   if (id1Tmp == 15) mME[0] = mtauME;
 419   mME[1] = 0.;
 420   int id2Tmp = abs(id2);
 421   if (id2Tmp ==  4) mME[1] = mcME;
 422   if (id2Tmp ==  5) mME[1] = mbME;
 423   if (id2Tmp == 13) mME[1] = mmuME;
 424   if (id2Tmp == 15) mME[1] = mtauME;
 425
 426   // If kinematically impossible return to massless case, but set error.
 427   if (mME[0] + mME[1] >= mH) {
 428     mME[0] = 0.;
 429     mME[1] = 0.;
 430     allowME = false;
 431   }
 432
 433   // Do incoming two-body kinematics for massless or massive cases.
 434   if (mME[0] == 0. && mME[1] == 0.) {
 435   pME[0] = 0.5 * mH * Vec4( 0., 0.,  1., 1.);
 436   pME[1] = 0.5 * mH * Vec4( 0., 0., -1., 1.);
 437   } else {
 438     double e0   = 0.5 * (mH * mH + mME[0] * mME[0] - mME[1] * mME[1]) / mH;
 439     double pz0  = sqrtpos(e0 * e0 - mME[0] * mME[0]);
 440     pME[0] = Vec4( 0., 0.,  pz0, e0);
 441     pME[1] = Vec4( 0., 0., -pz0, mH - e0);
 442   }
 443
 444   // Done.
 445   return allowME;
 446
 447 }
 448
 449 //--------------------------------------------------------------------------
 450
 451 // Evaluate weight for W decay distribution in t -> W b -> f fbar b.
 452
 453 double SigmaProcess::weightTopDecay( Event& process, int iResBeg,
 454   int iResEnd) {
 455
 456   // If not pair W d/s/b and mother t then return unit weight.
 457   if (iResEnd - iResBeg != 1) return 1.;
 458   int iW1  = iResBeg;
 459   int iB2  = iResBeg + 1;
 460   int idW1 = process[iW1].idAbs();
 461   int idB2 = process[iB2].idAbs();
 462   if (idW1 != 24) {
 463     swap(iW1, iB2);
 464     swap(idW1, idB2);
 465   }
 466   if (idW1 != 24 || (idB2 != 1 && idB2 != 3 && idB2 != 5)) return 1.;
 467   int iT   = process[iW1].mother1();
 468   if (iT <= 0 || process[iT].idAbs() != 6) return 1.;
 469
 470   // Find sign-matched order of W decay products.
 471   int iF    = process[iW1].daughter1();
 472   int iFbar = process[iW1].daughter2();
 473   if (iFbar - iF != 1) return 1.;
 474   if (process[iT].id() * process[iF].id() < 0) swap(iF, iFbar);
 475
 476   // Weight and maximum weight.
 477   double wt    = (process[iT].p() * process[iFbar].p())
 478                * (process[iF].p() * process[iB2].p());
 479   double wtMax = ( pow4(process[iT].m()) - pow4(process[iW1].m()) ) / 8.;
 480
 481   // Done.
 482   return wt / wtMax;
 483
 484 }
 485
 486 //--------------------------------------------------------------------------
 487
 488 // Evaluate weight for Z0/W+- decay distributions in H -> Z0/W+ Z0/W- -> 4f
 489 // and H -> gamma Z0 -> gamma f fbar.
 490
 491 double SigmaProcess::weightHiggsDecay( Event& process, int iResBeg,
 492   int iResEnd) {
 493
 494   // If not pair Z0 Z0, W+ W- or gamma Z0 then return unit weight.
 495   if (iResEnd - iResBeg != 1) return 1.;
 496   int iZW1  = iResBeg;
 497   int iZW2  = iResBeg + 1;
 498   int idZW1 = process[iZW1].id();
 499   int idZW2 = process[iZW2].id();
 500   if (idZW1 < 0 || idZW2 == 22) {
 501     swap(iZW1, iZW2);
 502     swap(idZW1, idZW2);
 503   }
 504   if ( (idZW1 != 23 || idZW2 != 23) && (idZW1 != 24 || idZW2 != -24)
 505     && (idZW1 != 22 || idZW2 != 23) ) return 1.;
 506
 507   // If mother is not Higgs then return unit weight.
 508   int iH  = process[iZW1].mother1();
 509   if (iH <= 0) return 1.;
 510   int idH = process[iH].id();
 511   if (idH != 25 && idH != 35 && idH !=36) return 1.;
 512
 513   // H -> gamma Z0 -> gamma f fbar is 1 + cos^2(theta) in Z rest frame.
 514   if (idZW1 == 22) {
 515     int i5 = process[iZW2].daughter1();
 516     int i6 = process[iZW2].daughter2();
 517     double pgmZ = process[iZW1].p() * process[iZW2].p();
 518     double pgm5 = process[iZW1].p() * process[i5].p();
 519     double pgm6 = process[iZW1].p() * process[i6].p();
 520     return (pow2(pgm5) + pow2(pgm6)) / pow2(pgmZ);
 521   }
 522
 523   // Parameters depend on Higgs type: H0(H_1), H^0(H_2) or A^0(H_3).
 524   int    higgsParity = higgsH1parity;
 525   double higgsEta    = higgsH1eta;
 526   if (idH == 35) {
 527     higgsParity      = higgsH2parity;
 528     higgsEta         = higgsH2eta;
 529   } else if (idH == 36) {
 530     higgsParity      = higgsA3parity;
 531     higgsEta         = higgsA3eta;
 532   }
 533
 534   // Option with isotropic decays.
 535   if (higgsParity == 0) return 1.;
 536
 537   // Maximum and initial weight.
 538   double wtMax = pow4(process[iH].m());
 539   double wt    = wtMax;
 540
 541   // Find sign-matched order of Z0/W+- decay products.
 542   int i3 = process[iZW1].daughter1();
 543   int i4 = process[iZW1].daughter2();
 544   if (process[i3].id() < 0) swap( i3, i4);
 545   int i5 = process[iZW2].daughter1();
 546   int i6 = process[iZW2].daughter2();
 547   if (process[i5].id() < 0) swap( i5, i6);
 548
 549   // Evaluate four-vector products and find masses..
 550   double p35  = 2. * process[i3].p() * process[i5].p();
 551   double p36  = 2. * process[i3].p() * process[i6].p();
 552   double p45  = 2. * process[i4].p() * process[i5].p();
 553   double p46  = 2. * process[i4].p() * process[i6].p();
 554   double p34  = 2. * process[i3].p() * process[i4].p();
 555   double p56  = 2. * process[i5].p() * process[i6].p();
 556   double mZW1 = process[iZW1].m();
 557   double mZW2 = process[iZW2].m();
 558
 559   // For mixed CP states need epsilon product and gauge boson masses.
 560   double epsilonProd = 0.;
 561   if (higgsParity == 3) {
 562     double p[4][4];
 563     for (int i = 0; i < 4; ++i) {
 564       int         ii = i3;
 565       if (i == 1) ii = i4;
 566       if (i == 2) ii = i5;
 567       if (i == 3) ii = i6;
 568       p[i][0] = process[ii].e();
 569       p[i][1] = process[ii].px();
 570       p[i][2] = process[ii].py();
 571       p[i][3] = process[ii].pz();
 572     }
 573     epsilonProd
 574       = p[0][0]*p[1][1]*p[2][2]*p[3][3] - p[0][0]*p[1][1]*p[2][3]*p[3][2]
 575       - p[0][0]*p[1][2]*p[2][1]*p[3][3] + p[0][0]*p[1][2]*p[2][3]*p[3][1]
 576       + p[0][0]*p[1][3]*p[2][1]*p[3][2] - p[0][0]*p[1][3]*p[2][2]*p[3][1]
 577       - p[0][1]*p[1][0]*p[2][2]*p[3][3] + p[0][1]*p[1][0]*p[2][3]*p[3][2]
 578       + p[0][1]*p[1][2]*p[2][0]*p[3][3] - p[0][1]*p[1][2]*p[2][3]*p[3][0]
 579       - p[0][1]*p[1][3]*p[2][0]*p[3][2] + p[0][1]*p[1][3]*p[2][2]*p[3][0]
 580       + p[0][2]*p[1][0]*p[2][1]*p[3][3] - p[0][2]*p[1][0]*p[2][3]*p[3][1]
 581       - p[0][2]*p[1][1]*p[2][0]*p[3][3] + p[0][2]*p[1][1]*p[2][3]*p[3][0]
 582       + p[0][2]*p[1][3]*p[2][0]*p[3][1] - p[0][2]*p[1][3]*p[2][1]*p[3][0]
 583       - p[0][3]*p[1][0]*p[2][1]*p[3][2] + p[0][3]*p[1][0]*p[2][2]*p[3][1]
 584       + p[0][3]*p[1][1]*p[2][0]*p[3][2] - p[0][3]*p[1][1]*p[2][2]*p[3][0]
 585       - p[0][3]*p[1][2]*p[2][0]*p[3][1] + p[0][3]*p[1][2]*p[2][1]*p[3][0];
 586   }
 587
 588   // Z0 Z0 decay: vector and axial couplings of two fermion pairs.
 589   if (idZW1 == 23) {
 590     double vf1 = couplingsPtr->vf(process[i3].idAbs());
 591     double af1 = couplingsPtr->af(process[i3].idAbs());
 592     double vf2 = couplingsPtr->vf(process[i5].idAbs());
 593     double af2 = couplingsPtr->af(process[i5].idAbs());
 594     double va12asym = 4. * vf1 * af1 * vf2 * af2
 595       / ( (vf1*vf1 + af1*af1) * (vf2*vf2 + af2*af2) );
 596     double etaMod = higgsEta / pow2( particleDataPtr->m0(23) );
 597
 598     // Normal CP-even decay.
 599     if (higgsParity == 1) wt = 8. * (1. + va12asym) * p35 * p46
 600       + 8. * (1. - va12asym) * p36 * p45;
 601
 602     // CP-odd decay (normal for A0(H_3)).
 603     else if (higgsParity == 2) wt = ( pow2(p35 + p46)
 604       + pow2(p36 + p45) - 2. * p34 * p56
 605       - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
 606       + va12asym * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) )
 607       / (1. +  va12asym);
 608
 609     // Mixed CP states.
 610     else wt = 32. * ( 0.25 * ( (1. + va12asym) * p35 * p46
 611       + (1. - va12asym) * p36 * p45 ) - 0.5 * etaMod * epsilonProd
 612       * ( (1. + va12asym) * (p35 + p46) - (1. - va12asym) * (p36 + p45) )
 613       + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
 614       - 2. * pow2(p35 * p46 - p36 * p45)
 615       + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
 616       + va12asym * p34 * p56 * (p35 + p36 - p45 - p46)
 617       * (p35 + p45 - p36 - p46) ) ) / ( 1. + 2. * etaMod * mZW1 * mZW2
 618       + 2. * pow2(etaMod * mZW1 * mZW2) * (1. + va12asym) );
 619
 620   // W+ W- decay.
 621   } else if (idZW1 == 24) {
 622     double etaMod = higgsEta / pow2( particleDataPtr->m0(24) );
 623
 624     // Normal CP-even decay.
 625     if (higgsParity == 1) wt = 16. * p35 * p46;
 626
 627     // CP-odd decay (normal for A0(H_3)).
 628     else if (higgsParity == 2) wt = 0.5 * ( pow2(p35 + p46)
 629       + pow2(p36 + p45) - 2. * p34 * p56
 630       - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
 631       + (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) );
 632
 633     // Mixed CP states.
 634     else wt = 32. * ( 0.25 * 2. * p35 * p46
 635       - 0.5 * etaMod * epsilonProd * 2. * (p35 + p46)
 636       + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
 637       - 2. * pow2(p35 * p46 - p36 * p45)
 638       + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
 639       + p34 * p56 * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) )
 640       / ( 1. * 2. * etaMod * mZW1 * mZW2 + 2. * pow2(etaMod * mZW1 * mZW2) );
 641   }
 642
 643   // Done.
 644   return wt / wtMax;
 645
 646 }
 647
 648 //==========================================================================
 649
 650 // The Sigma1Process class.
 651 // Base class for resolved 2 -> 1 cross sections; derived from SigmaProcess.
 652
 653 //--------------------------------------------------------------------------
 654
 655 // Wrapper to sigmaHat, to (a) store current incoming flavours,
 656 // (b) convert from GeV^-2 to mb where required, and
 657 // (c) convert from |M|^2 to d(sigmaHat)/d(tHat) where required.
 658
 659 double Sigma1Process::sigmaHatWrap(int id1in, int id2in) {
 660
 661   id1 = id1in;
 662   id2 = id2in;
 663   double sigmaTmp = sigmaHat();
 664   if (convertM2()) {
 665     sigmaTmp /= 2. * sH;
 666     // Convert 2 * pi * delta(p^2 - m^2) to Breit-Wigner with same area.
 667     int idTmp     = resonanceA();
 668     double mTmp   = particleDataPtr->m0(idTmp);
 669     double GamTmp = particleDataPtr->mWidth(idTmp);
 670     sigmaTmp *= 2. * mTmp * GamTmp / ( pow2(sH - mTmp * mTmp)
 671       + pow2(mTmp * GamTmp) );
 672   }
 673   if (convert2mb()) sigmaTmp *= CONVERT2MB;
 674   return sigmaTmp;
 675
 676 }
 677
 678 //--------------------------------------------------------------------------
 679
 680 // Input and complement kinematics for resolved 2 -> 1 process.
 681
 682 void Sigma1Process::store1Kin( double x1in, double x2in, double sHin) {
 683
 684   // Default value only sensible for these processes.
 685   swapTU = false;
 686
 687   // Incoming parton momentum fractions and sHat.
 688   x1Save = x1in;
 689   x2Save = x2in;
 690   sH     = sHin;
 691   mH     = sqrt(sH);
 692   sH2    = sH * sH;
 693
 694   // Different options for renormalization scale, but normally sHat.
 695   Q2RenSave                        = renormMultFac * sH;
 696   if (renormScale1 == 2) Q2RenSave = renormFixScale;
 697
 698   // Different options for factorization scale, but normally sHat.
 699   Q2FacSave                        = factorMultFac * sH;
 700   if (factorScale1 == 2) Q2FacSave = factorFixScale;
 701
 702   // Evaluate alpha_strong and alpha_EM.
 703   alpS   = couplingsPtr->alphaS(Q2RenSave);
 704   alpEM  = couplingsPtr->alphaEM(Q2RenSave);
 705
 706 }
 707
 708 //--------------------------------------------------------------------------
 709
 710 // Calculate modified masses and four-vectors for matrix elements.
 711
 712 bool Sigma1Process::setupForME() {
 713
 714   // Common initial-state handling.
 715   bool allowME = setupForMEin();
 716
 717   // Final state trivial here.
 718   mME[2] = mH;
 719   pME[2] = Vec4( 0., 0., 0., mH);
 720
 721   // Done.
 722   return allowME;
 723
 724 }
 725
 726 //==========================================================================
 727
 728 // The Sigma2Process class.
 729 // Base class for resolved 2 -> 2 cross sections; derived from SigmaProcess.
 730
 731 //--------------------------------------------------------------------------
 732
 733 // Input and complement kinematics for resolved 2 -> 2 process.
 734
 735 void Sigma2Process::store2Kin( double x1in, double x2in, double sHin,
 736   double tHin, double m3in, double m4in, double runBW3in, double runBW4in) {
 737
 738   // Default ordering of particles 3 and 4.
 739   swapTU   = false;
 740
 741   // Incoming parton momentum fractions.
 742   x1Save   = x1in;
 743   x2Save   = x2in;
 744
 745   // Incoming masses and their squares.
 746   bool masslessKin = (id3Mass() == 0) && (id4Mass() == 0);
 747   if (masslessKin) {
 748     m3     = 0.;
 749     m4     = 0.;
 750   } else {
 751     m3     = m3in;
 752     m4     = m4in;
 753   }
 754   mSave[3] = m3;
 755   mSave[4] = m4;
 756   s3       = m3 * m3;
 757   s4       = m4 * m4;
 758
 759   // Standard Mandelstam variables and their squares.
 760   sH       = sHin;
 761   tH       = tHin;
 762   uH       = (masslessKin) ? -(sH + tH) : s3 + s4 - (sH + tH);
 763   mH       = sqrt(sH);
 764   sH2      = sH * sH;
 765   tH2      = tH * tH;
 766   uH2      = uH * uH;
 767
 768   // The nominal Breit-Wigner factors with running width.
 769   runBW3   = runBW3in;
 770   runBW4   = runBW4in;
 771
 772   // Calculate squared transverse momentum.
 773   pT2 = (masslessKin) ?  tH * uH / sH : (tH * uH - s3 * s4) / sH;
 774
 775   // Special case: pick scale as if 2 -> 1 process in disguise.
 776   if (isSChannel()) {
 777
 778     // Different options for renormalization scale, but normally sHat.
 779     Q2RenSave                        = renormMultFac * sH;
 780     if (renormScale1 == 2) Q2RenSave = renormFixScale;
 781
 782     // Different options for factorization scale, but normally sHat.
 783     Q2FacSave                        = factorMultFac * sH;
 784     if (factorScale1 == 2) Q2FacSave = factorFixScale;
 785
 786   // Normal case with "true" 2 -> 2.
 787   } else {
 788
 789     // Different options for renormalization scale.
 790     if (masslessKin)            Q2RenSave = (renormScale2 < 4) ? pT2 : sH;
 791     else if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
 792     else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
 793     else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
 794     else                        Q2RenSave = sH;
 795     Q2RenSave                            *= renormMultFac;
 796     if      (renormScale2 == 5) Q2RenSave = renormFixScale;
 797
 798     // Different options for factorization scale.
 799     if (masslessKin)            Q2FacSave = (factorScale2 < 4) ? pT2 : sH;
 800     else if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
 801     else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
 802     else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
 803     else                        Q2FacSave = sH;
 804     Q2FacSave                            *= factorMultFac;
 805     if      (factorScale2 == 5) Q2FacSave = factorFixScale;
 806   }
 807
 808   // Evaluate alpha_strong and alpha_EM.
 809   alpS  = couplingsPtr->alphaS(Q2RenSave);
 810   alpEM = couplingsPtr->alphaEM(Q2RenSave);
 811
 812 }
 813
 814 //--------------------------------------------------------------------------
 815
 816 // As above, special kinematics for multiparton interactions.
 817
 818 void Sigma2Process::store2KinMPI( double x1in, double x2in,
 819   double sHin, double tHin, double uHin, double alpSin, double alpEMin,
 820   bool needMasses, double m3in, double m4in) {
 821
 822   // Default ordering of particles 3 and 4.
 823   swapTU    = false;
 824
 825   // Incoming x values.
 826   x1Save    = x1in;
 827   x2Save    = x2in;
 828
 829   // Standard Mandelstam variables and their squares.
 830   sH        = sHin;
 831   tH        = tHin;
 832   uH        = uHin;
 833   mH        = sqrt(sH);
 834   sH2       = sH * sH;
 835   tH2       = tH * tH;
 836   uH2       = uH * uH;
 837
 838   // Strong and electroweak couplings.
 839   alpS      = alpSin;
 840   alpEM     = alpEMin;
 841
 842   // Assume vanishing masses. (Will be modified in final kinematics.)
 843   m3        = 0.;
 844   s3        = 0.;
 845   m4        = 0.;
 846   s4        = 0.;
 847   sHBeta    = sH;
 848
 849   // Scattering angle.
 850   cosTheta  = (tH - uH) / sH;
 851   sinTheta  = 2. * sqrtpos( tH * uH ) / sH;
 852
 853   // In some cases must use masses and redefine meaning of tHat and uHat.
 854   if (needMasses) {
 855     m3      = m3in;
 856     s3      = m3 * m3;
 857     m4      = m4in;
 858     s4      = m4 * m4;
 859     sHMass  = sH - s3 - s4;
 860     sHBeta  = sqrtpos(sHMass*sHMass - 4. * s3 * s4);
 861     tH      = -0.5 * (sHMass - sHBeta * cosTheta);
 862     uH      = -0.5 * (sHMass + sHBeta * cosTheta);
 863     tH2     = tH * tH;
 864     uH2     = uH * uH;
 865   }
 866
 867   // pT2 with masses (at this stage) included.
 868   pT2Mass   = 0.25 * sHBeta * pow2(sinTheta);
 869
 870 }
 871
 872 //--------------------------------------------------------------------------
 873
 874 // Perform kinematics for a multiparton interaction, including a rescattering.
 875
 876 bool Sigma2Process::final2KinMPI( int i1Res, int i2Res, Vec4 p1Res, Vec4 p2Res,
 877   double m1Res, double m2Res) {
 878
 879   // Have to set flavours and colours.
 880   setIdColAcol();
 881
 882   // Check that masses of outgoing particles not too big.
 883   m3           = particleDataPtr->m0(idSave[3]);
 884   m4           = particleDataPtr->m0(idSave[4]);
 885   mH           = sqrt(sH);
 886   if (m3 + m4 + MASSMARGIN > mH) return false;
 887   s3           = m3 * m3;
 888   s4           = m4 * m4;
 889
 890   // Do kinematics of the production; without or with masses.
 891   double e1In  = 0.5 * mH;
 892   double e2In  = e1In;
 893   double pzIn  = e1In;
 894   if (i1Res > 0 || i2Res > 0) {
 895     double s1  = m1Res * m1Res;
 896     double s2  = m2Res * m2Res;
 897     e1In       = 0.5 * (sH + s1 - s2) / mH;
 898     e2In       = 0.5 * (sH + s2 - s1) / mH;
 899     pzIn       = sqrtpos( e1In*e1In - s1 );
 900   }
 901
 902   // Do kinematics of the decay.
 903   double e3    = 0.5 * (sH + s3 - s4) / mH;
 904   double e4    = 0.5 * (sH + s4 - s3) / mH;
 905   double pAbs  = sqrtpos( e3*e3 - s3 );
 906   phi          = 2. * M_PI * rndmPtr->flat();
 907   double pZ    = pAbs * cosTheta;
 908   pTFin        = pAbs * sinTheta;
 909   double pX    = pTFin * sin(phi);
 910   double pY    = pTFin * cos(phi);
 911   double scale = 0.5 * mH * sinTheta;
 912
 913   // Fill particle info.
 914   int status1  = (i1Res == 0) ? -31 : -34;
 915   int status2  = (i2Res == 0) ? -31 : -34;
 916   parton[1]    = Particle( idSave[1], status1, 0, 0, 3, 4,
 917     colSave[1], acolSave[1],  0.,  0.,  pzIn, e1In, m1Res, scale);
 918   parton[2]    = Particle( idSave[2], status2, 0, 0, 3, 4,
 919     colSave[2], acolSave[2],  0.,  0., -pzIn, e2In, m2Res, scale);
 920   parton[3]    = Particle( idSave[3],      33, 1, 2, 0, 0,
 921     colSave[3], acolSave[3],  pX,  pY,    pZ,   e3,    m3, scale);
 922   parton[4]    = Particle( idSave[4],      33, 1, 2, 0, 0,
 923     colSave[4], acolSave[4], -pX, -pY,   -pZ,   e4,    m4, scale);
 924
 925   // Boost particles from subprocess rest frame to event rest frame.
 926   // Normal multiparton interaction: only longitudinal boost.
 927   if (i1Res == 0 && i2Res == 0) {
 928     double betaZ = (x1Save - x2Save) / (x1Save + x2Save);
 929     for (int i = 1; i <= 4; ++i) parton[i].bst(0., 0., betaZ);
 930   // Rescattering: generic rotation and boost required.
 931   } else {
 932     RotBstMatrix M;
 933     M.fromCMframe( p1Res, p2Res);
 934     for (int i = 1; i <= 4; ++i) parton[i].rotbst(M);
 935   }
 936
 937   // Done.
 938   return true;
 939
 940 }
 941
 942 //--------------------------------------------------------------------------
 943
 944 // Calculate modified masses and four-vectors for matrix elements.
 945
 946 bool Sigma2Process::setupForME() {
 947
 948   // Common initial-state handling.
 949   bool allowME = setupForMEin();
 950
 951   // Correct outgoing c, b, mu and tau to be massive or not.
 952   mME[2] = m3;
 953   int id3Tmp = abs(id3Mass());
 954   if (id3Tmp ==  4) mME[2] = mcME;
 955   if (id3Tmp ==  5) mME[2] = mbME;
 956   if (id3Tmp == 13) mME[2] = mmuME;
 957   if (id3Tmp == 15) mME[2] = mtauME;
 958   mME[3] = m4;
 959   int id4Tmp = abs(id4Mass());
 960   if (id4Tmp ==  4) mME[3] = mcME;
 961   if (id4Tmp ==  5) mME[3] = mbME;
 962   if (id4Tmp == 13) mME[3] = mmuME;
 963   if (id4Tmp == 15) mME[3] = mtauME;
 964
 965   // If kinematically impossible turn to massless case, but set error.
 966   if (mME[2] + mME[3] >= mH) {
 967     mME[2] = 0.;
 968     mME[3] = 0.;
 969     allowME = false;
 970   }
 971
 972   // Calculate scattering angle in subsystem rest frame.
 973   double sH34 = sqrtpos( pow2(sH - s3 - s4) - 4. * s3 * s4);
 974   double cThe = (tH - uH) / sH34;
 975   double sThe = sqrtpos(1. - cThe * cThe);
 976
 977   // Setup massive kinematics with preserved scattering angle.
 978   double s3ME   = pow2(mME[2]);
 979   double s4ME   = pow2(mME[3]);
 980   double sH34ME = sqrtpos( pow2(sH - s3ME - s4ME) - 4. * s3ME * s4ME);
 981   double pAbsME = 0.5 * sH34ME / mH;
 982
 983   // Normally allowed with unequal (or vanishing) masses.
 984   if (id3Tmp == 0 || id3Tmp != id4Tmp) {
 985     pME[2] = Vec4(  pAbsME * sThe, 0.,  pAbsME * cThe,
 986              0.5 * (sH + s3ME - s4ME) / mH);
 987     pME[3] = Vec4( -pAbsME * sThe, 0., -pAbsME * cThe,
 988              0.5 * (sH + s4ME - s3ME) / mH);
 989
 990   // For equal (anti)particles (e.g. W+ W-) use averaged mass.
 991   } else {
 992     mME[2] = sqrtpos(0.5 * (s3ME + s4ME) - 0.25 * pow2(s3ME - s4ME) / sH);
 993     mME[3] = mME[2];
 994     pME[2] = Vec4(  pAbsME * sThe, 0.,  pAbsME * cThe, 0.5 * mH);
 995     pME[3] = Vec4( -pAbsME * sThe, 0., -pAbsME * cThe, 0.5 * mH);
 996   }
 997
 998   // Done.
 999   return allowME;
1000
1001 }
1002
1003 //==========================================================================
1004
1005 // The Sigma3Process class.
1006 // Base class for resolved 2 -> 3 cross sections; derived from SigmaProcess.
1007
1008 //--------------------------------------------------------------------------
1009
1010 // Input and complement kinematics for resolved 2 -> 3 process.
1011
1012 void Sigma3Process::store3Kin( double x1in, double x2in, double sHin,
1013   Vec4 p3cmIn, Vec4 p4cmIn, Vec4 p5cmIn, double m3in, double m4in,
1014   double m5in, double runBW3in, double runBW4in, double runBW5in) {
1015
1016   // Default ordering of particles 3 and 4 - not relevant here.
1017   swapTU   = false;
1018
1019   // Incoming parton momentum fractions.
1020   x1Save   = x1in;
1021   x2Save   = x2in;
1022
1023   // Incoming masses and their squares.
1024   if (id3Mass() == 0 && id4Mass() == 0 && id5Mass() == 0) {
1025     m3     = 0.;
1026     m4     = 0.;
1027     m5     = 0.;
1028   } else {
1029     m3     = m3in;
1030     m4     = m4in;
1031     m5     = m5in;
1032   }
1033   mSave[3] = m3;
1034   mSave[4] = m4;
1035   mSave[5] = m5;
1036   s3       = m3 * m3;
1037   s4       = m4 * m4;
1038   s5       = m5 * m5;
1039
1040   // Standard Mandelstam variables and four-momenta in rest frame.
1041   sH       = sHin;
1042   mH       = sqrt(sH);
1043   sH2      = sH * sH;
1044   p3cm     = p3cmIn;
1045   p4cm     = p4cmIn;
1046   p5cm     = p5cmIn;
1047
1048   // The nominal Breit-Wigner factors with running width.
1049   runBW3   = runBW3in;
1050   runBW4   = runBW4in;
1051   runBW5   = runBW5in;
1052
1053   // Special case: pick scale as if 2 -> 1 process in disguise.
1054   if (isSChannel()) {
1055
1056     // Different options for renormalization scale, but normally sHat.
1057     Q2RenSave = renormMultFac * sH;
1058     if (renormScale1 == 2) Q2RenSave = renormFixScale;
1059
1060     // Different options for factorization scale, but normally sHat.
1061     Q2FacSave = factorMultFac * sH;
1062     if (factorScale1 == 2) Q2RenSave = factorFixScale;
1063
1064   // "Normal" 2 -> 3 processes, i.e. not vector boson fusion.
1065   } else if ( idTchan1() != 23 && idTchan1() != 24 && idTchan2() != 23
1066     && idTchan2() != 24 ) {
1067     double mT3S = s3 + p3cm.pT2();
1068     double mT4S = s4 + p4cm.pT2();
1069     double mT5S = s5 + p5cm.pT2();
1070
1071     // Different options for renormalization scale.
1072     if      (renormScale3 == 1) Q2RenSave = min( mT3S, min(mT4S, mT5S) );
1073     else if (renormScale3 == 2) Q2RenSave = sqrt( mT3S * mT4S * mT5S
1074       / max( mT3S, max(mT4S, mT5S) ) );
1075     else if (renormScale3 == 3) Q2RenSave = pow( mT3S * mT4S * mT5S,
1076                                             1./3. );
1077     else if (renormScale3 == 4) Q2RenSave = (mT3S + mT4S + mT5S) / 3.;
1078     else                        Q2RenSave = sH;
1079     Q2RenSave                            *= renormMultFac;
1080     if      (renormScale3 == 6) Q2RenSave = renormFixScale;
1081
1082     // Different options for factorization scale.
1083     if      (factorScale3 == 1) Q2FacSave = min( mT3S, min(mT4S, mT5S) );
1084     else if (factorScale3 == 2) Q2FacSave = sqrt( mT3S * mT4S * mT5S
1085       / max( mT3S, max(mT4S, mT5S) ) );
1086     else if (factorScale3 == 3) Q2FacSave = pow( mT3S * mT4S * mT5S,
1087                                             1./3. );
1088     else if (factorScale3 == 4) Q2FacSave = (mT3S + mT4S + mT5S) / 3.;
1089     else                        Q2FacSave = sH;
1090     Q2FacSave                            *= factorMultFac;
1091     if      (factorScale3 == 6) Q2FacSave = factorFixScale;
1092
1093   // Vector boson fusion 2 -> 3 processes; recoils in positions 4 and 5.
1094   } else {
1095     double sV4   = pow2( particleDataPtr->m0(idTchan1()) );
1096     double sV5   = pow2( particleDataPtr->m0(idTchan2()) );
1097     double mT3S  = s3  + p3cm.pT2();
1098     double mTV4S = sV4 + p4cm.pT2();
1099     double mTV5S = sV5 + p5cm.pT2();
1100
1101     // Different options for renormalization scale.
1102     if      (renormScale3VV == 1) Q2RenSave = max( sV4, sV5);
1103     else if (renormScale3VV == 2) Q2RenSave = sqrt( mTV4S * mTV5S );
1104     else if (renormScale3VV == 3) Q2RenSave = pow( mT3S * mTV4S * mTV5S,
1105                                               1./3. );
1106     else if (renormScale3VV == 4) Q2RenSave = (mT3S * mTV4S * mTV5S) / 3.;
1107     else                          Q2RenSave = sH;
1108     Q2RenSave                              *= renormMultFac;
1109     if      (renormScale3VV == 6) Q2RenSave = renormFixScale;
1110
1111     // Different options for factorization scale.
1112     if      (factorScale3VV == 1) Q2FacSave = max( sV4, sV5);
1113     else if (factorScale3VV == 2) Q2FacSave = sqrt( mTV4S * mTV5S );
1114     else if (factorScale3VV == 3) Q2FacSave = pow( mT3S * mTV4S * mTV5S,
1115                                               1./3. );
1116     else if (factorScale3VV == 4) Q2FacSave = (mT3S * mTV4S * mTV5S) / 3.;
1117     else                          Q2FacSave = sH;
1118     Q2FacSave                              *= factorMultFac;
1119     if      (factorScale3VV == 6) Q2FacSave = factorFixScale;
1120   }
1121
1122   // Evaluate alpha_strong and alpha_EM.
1123   alpS  = couplingsPtr->alphaS(Q2RenSave);
1124   alpEM = couplingsPtr->alphaEM(Q2RenSave);
1125
1126 }
1127
1128 //--------------------------------------------------------------------------
1129
1130 // Calculate modified masses and four-vectors for matrix elements.
1131
1132 bool Sigma3Process::setupForME() {
1133
1134   // Common initial-state handling.
1135   bool allowME = setupForMEin();
1136
1137   // Correct outgoing c, b, mu and tau to be massive or not.
1138   mME[2] = m3;
1139   int id3Tmp = abs(id3Mass());
1140   if (id3Tmp ==  4) mME[2] = mcME;
1141   if (id3Tmp ==  5) mME[2] = mbME;
1142   if (id3Tmp == 13) mME[2] = mmuME;
1143   if (id3Tmp == 15) mME[2] = mtauME;
1144   mME[3] = m4;
1145   int id4Tmp = abs(id4Mass());
1146   if (id4Tmp ==  4) mME[3] = mcME;
1147   if (id4Tmp ==  5) mME[3] = mbME;
1148   if (id4Tmp == 13) mME[3] = mmuME;
1149   if (id4Tmp == 15) mME[3] = mtauME;
1150   mME[4] = m5;
1151   int id5Tmp = abs(id5Mass());
1152   if (id5Tmp ==  4) mME[4] = mcME;
1153   if (id5Tmp ==  5) mME[4] = mbME;
1154   if (id5Tmp == 13) mME[4] = mmuME;
1155   if (id5Tmp == 15) mME[4] = mtauME;
1156
1157   // If kinematically impossible turn to massless case, but set error.
1158   if (mME[2] + mME[3] + mME[4] >= mH) {
1159     mME[2] = 0.;
1160     mME[3] = 0.;
1161     mME[4] = 0.;
1162     allowME = false;
1163   }
1164
1165   // Form new average masses if identical particles.
1166   if (id3Tmp != 0 && id4Tmp == id3Tmp && id5Tmp == id3Tmp) {
1167     double mAvg = (mME[2] + mME[3] + mME[4]) / 3.;
1168     mME[2] = mAvg;
1169     mME[3] = mAvg;
1170     mME[4] = mAvg;
1171   } else if (id3Tmp != 0 && id4Tmp == id3Tmp) {
1172     mME[2] = sqrtpos(0.5 * (pow2(mME[2]) + pow2(mME[3]))
1173            - 0.25 * pow2(pow2(mME[2]) - pow2(mME[3])) / sH);
1174     mME[3] = mME[2];
1175   } else if (id3Tmp != 0 && id5Tmp == id3Tmp) {
1176     mME[2] = sqrtpos(0.5 * (pow2(mME[2]) + pow2(mME[4]))
1177            - 0.25 * pow2(pow2(mME[2]) - pow2(mME[4])) / sH);
1178     mME[4] = mME[2];
1179   } else if (id4Tmp != 0 && id5Tmp == id4Tmp) {
1180     mME[3] = sqrtpos(0.5 * (pow2(mME[3]) + pow2(mME[4]))
1181            - 0.25 * pow2(pow2(mME[3]) - pow2(mME[4])) / sH);
1182     mME[4] = mME[2];
1183   }
1184
1185   // Iterate rescaled three-momenta until convergence.
1186   double m2ME3 = pow2(mME[2]);
1187   double m2ME4 = pow2(mME[3]);
1188   double m2ME5 = pow2(mME[4]);
1189   double p2ME3 = p3cm.pAbs2();
1190   double p2ME4 = p4cm.pAbs2();
1191   double p2ME5 = p5cm.pAbs2();
1192   double p2sum = p2ME3 + p2ME4 + p2ME5;
1193   double eME3  = sqrt(m2ME3 + p2ME3);
1194   double eME4  = sqrt(m2ME4 + p2ME4);
1195   double eME5  = sqrt(m2ME5 + p2ME5);
1196   double esum  = eME3 + eME4 + eME5;
1197   double p2rat = p2ME3 / eME3 + p2ME4 / eME4 + p2ME5 / eME5;
1198   int iStep = 0;
1199   while ( abs(esum - mH) > COMPRELERR * mH && iStep < NCOMPSTEP ) {
1200     ++iStep;
1201     double compFac = 1. + 2. * (mH - esum) / p2rat;
1202     p2ME3 *= compFac;
1203     p2ME4 *= compFac;
1204     p2ME5 *= compFac;
1205     eME3   = sqrt(m2ME3 + p2ME3);
1206     eME4   = sqrt(m2ME4 + p2ME4);
1207     eME5   = sqrt(m2ME5 + p2ME5);
1208     esum   = eME3 + eME4 + eME5;
1209     p2rat  = p2ME3 / eME3 + p2ME4 / eME4 + p2ME5 / eME5;
1210   }
1211
1212   // If failed convergence set error flag.
1213   if (abs(esum - mH) > COMPRELERR * mH) allowME = false;
1214
1215   // Set up accepted kinematics.
1216   double totFac = sqrt( (p2ME3 + p2ME4 + p2ME5) / p2sum);
1217   pME[2] = totFac * p3cm;
1218   pME[2].e( eME3);
1219   pME[3] = totFac * p4cm;
1220   pME[3].e( eME4);
1221   pME[4] = totFac * p5cm;
1222   pME[4].e( eME5);
1223
1224   // Done.
1225   return allowME;
1226
1227 }
1228
1229 //==========================================================================
1230
1231 // The SigmaLHAProcess class.
1232 // Wrapper for Les Houches Accord external input; derived from SigmaProcess.
1233 // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
1234
1235 //--------------------------------------------------------------------------
1236
1237 // Evaluate weight for decay angles.
1238
1239 double SigmaLHAProcess::weightDecay( Event& process, int iResBeg,
1240   int iResEnd) {
1241
1242   // Do nothing if decays present already at input.
1243   if (iResBeg < process.savedSizeValue()) return 1.;
1244
1245   // Identity of mother of decaying reseonance(s).
1246   int idMother = process[process[iResBeg].mother1()].idAbs();
1247
1248   // For Higgs decay hand over to standard routine.
1249   if (idMother == 25 || idMother == 35 || idMother == 36)
1250     return weightHiggsDecay( process, iResBeg, iResEnd);
1251
1252   // For top decay hand over to standard routine.
1253   if (idMother == 6)
1254     return weightTopDecay( process, iResBeg, iResEnd);
1255
1256   // Else done.
1257   return 1.;
1258
1259 }
1260
1261 //--------------------------------------------------------------------------
1262
1263 // Set scale, alpha_strong and alpha_EM when not set.
1264
1265 void SigmaLHAProcess::setScale() {
1266
1267   // If scale has not been set, then to set.
1268   double scaleLHA = lhaUpPtr->scale();
1269   if (scaleLHA < 0.) {
1270
1271     // Final-state partons and their invariant mass.
1272     vector<int> iFin;
1273     Vec4 pFinSum;
1274     for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
1275     if (lhaUpPtr->mother1(i) == 1) {
1276       iFin.push_back(i);
1277       pFinSum += Vec4( lhaUpPtr->px(i), lhaUpPtr->py(i),
1278         lhaUpPtr->pz(i), lhaUpPtr->e(i) );
1279     }
1280     int nFin = iFin.size();
1281     sH       = pFinSum * pFinSum;
1282     mH       = sqrt(sH);
1283     sH2      = sH * sH;
1284
1285     // If 1 final-state particle then use Sigma1Process logic.
1286     if (nFin == 1) {
1287       Q2RenSave                             = renormMultFac * sH;
1288       if (renormScale1 == 2) Q2RenSave      = renormFixScale;
1289       Q2FacSave                             = factorMultFac * sH;
1290       if (factorScale1 == 2) Q2FacSave      = factorFixScale;
1291
1292     // If 2 final-state particles then use Sigma2Process logic.
1293     } else if (nFin == 2) {
1294       double s3  = pow2(lhaUpPtr->m(iFin[0]));
1295       double s4  = pow2(lhaUpPtr->m(iFin[1]));
1296       double pT2 = pow2(lhaUpPtr->px(iFin[0])) + pow2(lhaUpPtr->py(iFin[0]));
1297       if      (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
1298       else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
1299       else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
1300       else                        Q2RenSave = sH;
1301       Q2RenSave                            *= renormMultFac;
1302       if      (renormScale2 == 5) Q2RenSave = renormFixScale;
1303       if      (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
1304       else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
1305       else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
1306       else                        Q2FacSave = sH;
1307       Q2FacSave                            *= factorMultFac;
1308       if      (factorScale2 == 5) Q2FacSave = factorFixScale;
1309
1310     // If 3 or more final-state particles then use Sigma3Process logic.
1311     } else {
1312       double mTSlow  = sH;
1313       double mTSmed  = sH;
1314       double mTSprod = 1.;
1315       double mTSsum  = 0.;
1316       for (int i = 0; i < nFin; ++i) {
1317         double mTSnow = pow2(lhaUpPtr->m(iFin[i]))
1318           + pow2(lhaUpPtr->px(iFin[i])) + pow2(lhaUpPtr->py(iFin[i]));
1319         if      (mTSnow < mTSlow) {mTSmed = mTSlow; mTSlow = mTSnow;}
1320         else if (mTSnow < mTSmed) mTSmed = mTSnow;
1321         mTSprod *= mTSnow;
1322         mTSsum  += mTSnow;
1323       }
1324       if      (renormScale3 == 1) Q2RenSave = mTSlow;
1325       else if (renormScale3 == 2) Q2RenSave = sqrt(mTSlow * mTSmed);
1326       else if (renormScale3 == 3) Q2RenSave = pow(mTSprod, 1. / nFin);
1327       else if (renormScale3 == 4) Q2RenSave = mTSsum / nFin;
1328       else                        Q2RenSave = sH;
1329       Q2RenSave                            *= renormMultFac;
1330       if      (renormScale3 == 6) Q2RenSave = renormFixScale;
1331       if      (factorScale3 == 1) Q2FacSave = mTSlow;
1332       else if (factorScale3 == 2) Q2FacSave = sqrt(mTSlow * mTSmed);
1333       else if (factorScale3 == 3) Q2FacSave = pow(mTSprod, 1. / nFin);
1334       else if (factorScale3 == 4) Q2FacSave = mTSsum / nFin;
1335       else                        Q2FacSave = sH;
1336       Q2FacSave                            *= factorMultFac;
1337       if      (factorScale3 == 6) Q2FacSave = factorFixScale;
1338     }
1339   }
1340
1341   // If alpha_strong and alpha_EM have not been set, then set them.
1342   if (lhaUpPtr->alphaQCD() < 0.001) {
1343     double Q2RenNow = (scaleLHA < 0.) ? Q2RenSave : pow2(scaleLHA);
1344     alpS = couplingsPtr->alphaS(Q2RenNow);
1345   }
1346   if (lhaUpPtr->alphaQED() < 0.001) {
1347     double Q2RenNow = (scaleLHA < 0.) ? Q2RenSave : pow2(scaleLHA);
1348     alpEM = couplingsPtr->alphaEM(Q2RenNow);
1349   }
1350
1351 }
1352
1353 //--------------------------------------------------------------------------
1354
1355 // Obtain number of final-state partons from LHA object.
1356
1357 int SigmaLHAProcess::nFinal() const {
1358
1359   // At initialization size unknown, so return 0.
1360   if (lhaUpPtr->sizePart() <= 0) return 0;
1361
1362   // Sum up all particles that has first mother = 1.
1363   int nFin = 0;
1364   for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
1365     if (lhaUpPtr->mother1(i) == 1) ++nFin;
1366   return nFin;
1367
1368 }
1369
1370 //==========================================================================
1371
1372 } // end namespace Pythia8