PYTHIA8/pythia8130/src/SigmaProcess.cxx

   1 // SigmaProcess.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   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 //*********
  30
  31 // Perform simple initialization and store pointers.
  32
  33 void SigmaProcess::init(Info* infoPtrIn, BeamParticle* beamAPtrIn,
  34   BeamParticle* beamBPtrIn, AlphaStrong* alphaSPtrIn,
  35   AlphaEM* alphaEMPtrIn,SigmaTotal* sigmaTotPtrIn,
  36   SusyLesHouches* slhaPtrIn) {
  37
  38   // Store pointers.
  39   infoPtr     = infoPtrIn;
  40   beamAPtr    = beamAPtrIn;
  41   beamBPtr    = beamBPtrIn;
  42   alphaSPtr   = alphaSPtrIn;
  43   alphaEMPtr  = alphaEMPtrIn;
  44   sigmaTotPtr = sigmaTotPtrIn;
  45   slhaPtr     = slhaPtrIn;
  46
  47   // Read out some properties of beams to allow shorthand.
  48   idA         = beamAPtr->id();
  49   idB         = beamBPtr->id();
  50   mA          = beamAPtr->m();
  51   mB          = beamBPtr->m();
  52   isLeptonA   = beamAPtr->isLepton();
  53   isLeptonB   = beamBPtr->isLepton();
  54   hasLeptonBeams = isLeptonA || isLeptonB;
  55
  56   // K factor, multiplying resolved processes. (But not here for MI.)
  57   Kfactor         = Settings::parm("SigmaProcess:Kfactor");
  58
  59   // Maximum incoming quark flavour.
  60   nQuarkIn        = Settings::mode("PDFinProcess:nQuarkIn");
  61
  62   // Renormalization scale choice.
  63   renormScale1    = Settings::mode("SigmaProcess:renormScale1");
  64   renormScale2    = Settings::mode("SigmaProcess:renormScale2");
  65   renormScale3    = Settings::mode("SigmaProcess:renormScale3");
  66   renormScale3VV  = Settings::mode("SigmaProcess:renormScale3VV");
  67   renormMultFac   = Settings::parm("SigmaProcess:renormMultFac");
  68   renormFixScale  = Settings::parm("SigmaProcess:renormFixScale");
  69
  70   // Factorization scale choice.
  71   factorScale1    = Settings::mode("SigmaProcess:factorScale1");
  72   factorScale2    = Settings::mode("SigmaProcess:factorScale2");
  73   factorScale3    = Settings::mode("SigmaProcess:factorScale3");
  74   factorScale3VV  = Settings::mode("SigmaProcess:factorScale3VV");
  75   factorMultFac   = Settings::parm("SigmaProcess:factorMultFac");
  76   factorFixScale  = Settings::parm("SigmaProcess:factorFixScale");
  77
  78   // CP violation parameters for the BSM Higgs sector.
  79   higgsH1parity   = Settings::mode("HiggsH1:parity");
  80   higgsH1eta      = Settings::parm("HiggsH1:etaParity");
  81   higgsH2parity   = Settings::mode("HiggsH2:parity");
  82   higgsH2eta      = Settings::parm("HiggsH2:etaParity");
  83   higgsA3parity   = Settings::mode("HiggsA3:parity");
  84   higgsA3eta      = Settings::parm("HiggsA3:etaParity");
  85
  86   // If BSM not switched on then H1 should have SM properties.
  87   if (!Settings::flag("Higgs:useBSM")){
  88     higgsH1parity = 1;
  89     higgsH1eta    = 0.;
  90   }
  91
  92 }
  93
  94 //*********
  95
  96 // Set up allowed flux of incoming partons.
  97 // addBeam: set up PDF's that need to be evaluated for the two beams.
  98 // addPair: set up pairs of incoming partons from the two beams.
  99
 100 bool SigmaProcess::initFlux() {
 101
 102   // Read in process-specific channel information.
 103   string fluxType = inFlux();
 104
 105   // Case with g g incoming state.
 106   if (fluxType == "gg") {
 107     addBeamA(21);
 108     addBeamB(21);
 109     addPair(21, 21);
 110   }
 111
 112   // Case with q g incoming state.
 113   else if (fluxType == "qg") {
 114     for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
 115       int idNow = (i == 0) ? 21 : i;
 116       addBeamA(idNow);
 117       addBeamB(idNow);
 118     }
 119     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 120     if (idNow != 0) {
 121       addPair(idNow, 21);
 122       addPair(21, idNow);
 123     }
 124   }
 125
 126   // Case with q q', q qbar' or qbar qbar' incoming state.
 127   else if (fluxType == "qq") {
 128     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 129     if (idNow != 0) {
 130       addBeamA(idNow);
 131       addBeamB(idNow);
 132     }
 133     for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 134     if (id1Now != 0)
 135     for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 136     if (id2Now != 0)
 137       addPair(id1Now, id2Now);
 138   }
 139
 140   // Case with q qbar incoming state.
 141   else if (fluxType == "qqbarSame") {
 142     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 143     if (idNow != 0) {
 144       addBeamA(idNow);
 145       addBeamB(idNow);
 146     }
 147     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 148     if (idNow != 0)
 149       addPair(idNow, -idNow);
 150   }
 151
 152   // Case with f f', f fbar', fbar fbar' incoming state.
 153   else if (fluxType == "ff") {
 154     // If beams are leptons then they are also the colliding partons.
 155     if ( isLeptonA && isLeptonB ) {
 156       addBeamA(idA);
 157       addBeamB(idB);
 158       addPair(idA, idB);
 159     // First beam is lepton and second is hadron.
 160     } else if ( isLeptonA ) {
 161       addBeamA(idA);
 162       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 163       if (idNow != 0) {
 164         addBeamB(idNow);
 165         addPair(idA, idNow);
 166       }
 167     // First beam is hadron and second is lepton.
 168     } else if ( isLeptonB ) {
 169       addBeamB(idB);
 170       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 171       if (idNow != 0) {
 172         addBeamA(idNow);
 173         addPair(idNow, idB);
 174       }
 175     // Hadron beams gives quarks.
 176     } else {
 177       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 178       if (idNow != 0) {
 179         addBeamA(idNow);
 180         addBeamB(idNow);
 181       }
 182       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 183       if (id1Now != 0)
 184       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 185       if (id2Now != 0)
 186         addPair(id1Now, id2Now);
 187     }
 188   }
 189
 190   // Case with f fbar incoming state.
 191   else if (fluxType == "ffbarSame") {
 192     // If beams are antiparticle pair and leptons then also colliding partons.
 193     if ( idA + idB == 0 && isLeptonA ) {
 194       addBeamA(idA);
 195       addBeamB(idB);
 196       addPair(idA, idB);
 197     // Else assume both to be hadrons, for better or worse.
 198     } else {
 199       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 200       if (idNow != 0) {
 201         addBeamA(idNow);
 202         addBeamB(idNow);
 203       }
 204       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 205       if (idNow != 0)
 206         addPair(idNow, -idNow);
 207     }
 208   }
 209
 210   // Case with f fbar' charged(+-1) incoming state.
 211   else if (fluxType == "ffbarChg") {
 212     // If beams are leptons then also colliding partons.
 213     if ( isLeptonA && isLeptonB && abs( ParticleDataTable::chargeType(idA)
 214              + ParticleDataTable::chargeType(idB) ) == 3 ) {
 215       addBeamA(idA);
 216       addBeamB(idB);
 217       addPair(idA, idB);
 218     // Hadron beams gives quarks.
 219     } else {
 220       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 221       if (idNow != 0) {
 222         addBeamA(idNow);
 223         addBeamB(idNow);
 224       }
 225       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
 226       if (id1Now != 0)
 227       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
 228       if (id2Now != 0 && id1Now * id2Now < 0
 229         && (abs(id1Now) + abs(id2Now))%2 == 1) addPair(id1Now, id2Now);
 230     }
 231   }
 232
 233   // Case with f fbar' generic incoming state.
 234   else if (fluxType == "ffbar") {
 235     // If beams are leptons then also colliding partons.
 236     if (isLeptonA && isLeptonB && idA * idB < 0) {
 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         addPair(id1Now, id2Now);
 252     }
 253   }
 254
 255   // Case with f gamma incoming state.
 256   else if (fluxType == "fgm") {
 257     // Fermion from incoming side A.
 258     if ( isLeptonA ) {
 259       addBeamA(idA);
 260       addPair(idA, 22);
 261     } else {
 262       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 263       if (idNow != 0) {
 264         addBeamA(idNow);
 265         addPair(idNow, 22);
 266       }
 267     }
 268     // Fermion from incoming side B.
 269     if ( isLeptonB ) {
 270       addBeamB( idB);
 271       addPair(22, idB);
 272     } else {
 273       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
 274       if (idNow != 0) {
 275         addBeamB(idNow);
 276         addPair(22, idNow);
 277       }
 278     }
 279     // Photons in the beams.
 280     addBeamA(22);
 281     addBeamB(22);
 282   }
 283
 284   // Case with gamma gamma incoming state.
 285   else if (fluxType == "ggm") {
 286     addBeamA(21);
 287     addBeamA(22);
 288     addBeamB(21);
 289     addBeamB(22);
 290     addPair(21, 22);
 291     addPair(22, 21);
 292   }
 293
 294   // Case with gamma gamma incoming state.
 295   else if (fluxType == "gmgm") {
 296     addBeamA(22);
 297     addBeamB(22);
 298     addPair(22, 22);
 299   }
 300
 301   // Unrecognized fluxType is bad sign. Else done.
 302   else {
 303     infoPtr->errorMsg("Error in SigmaProcess::initFlux: "
 304     "unrecognized inFlux type", fluxType);
 305     return false;
 306   }
 307   return true;
 308
 309 }
 310
 311 //*********
 312
 313 // Convolute matrix-element expression(s) with parton flux and K factor.
 314
 315 double SigmaProcess::sigmaPDF() {
 316
 317   // Evaluate and store the required parton densities.
 318   for (int j = 0; j < sizeBeamA(); ++j)
 319     inBeamA[j].pdf = beamAPtr->xfHard( inBeamA[j].id, x1Save, Q2FacSave);
 320   for (int j = 0; j < sizeBeamB(); ++j)
 321     inBeamB[j].pdf = beamBPtr->xfHard( inBeamB[j].id, x2Save, Q2FacSave);
 322
 323   // Loop over allowed incoming channels.
 324   sigmaSumSave = 0.;
 325   for (int i = 0; i < sizePair(); ++i) {
 326
 327     // Evaluate hard-scattering cross section. Include K factor.
 328     inPair[i].pdfSigma = Kfactor
 329                        * sigmaHatWrap(inPair[i].idA, inPair[i].idB);
 330
 331     // Multiply by respective parton densities.
 332     for (int j = 0; j < sizeBeamA(); ++j)
 333     if (inPair[i].idA == inBeamA[j].id) {
 334       inPair[i].pdfA      = inBeamA[j].pdf;
 335       inPair[i].pdfSigma *= inBeamA[j].pdf;
 336       break;
 337     }
 338     for (int j = 0; j < sizeBeamB(); ++j)
 339     if (inPair[i].idB == inBeamB[j].id) {
 340       inPair[i].pdfB      = inBeamB[j].pdf;
 341       inPair[i].pdfSigma *= inBeamB[j].pdf;
 342       break;
 343     }
 344
 345     // Sum for all channels.
 346     sigmaSumSave += inPair[i].pdfSigma;
 347   }
 348
 349   // Done.
 350   return sigmaSumSave;
 351
 352 }
 353
 354 //*********
 355
 356 // Select incoming parton channel and extract parton densities (resolved).
 357
 358 void SigmaProcess::pickInState(int id1in, int id2in) {
 359
 360   // Multiple interactions: partons already selected.
 361   if (id1in != 0 && id2in != 0) {
 362     id1 = id1in;
 363     id2 = id2in;
 364   }
 365
 366   // Pick channel. Extract channel flavours and pdf's.
 367   double sigmaRand =  sigmaSumSave * Rndm::flat();
 368   for (int i = 0; i < sizePair(); ++i) {
 369     sigmaRand -= inPair[i].pdfSigma;
 370     if (sigmaRand <= 0.) {
 371       id1      = inPair[i].idA;
 372       id2      = inPair[i].idB;
 373       pdf1Save = inPair[i].pdfA;
 374       pdf2Save = inPair[i].pdfB;
 375       break;
 376     }
 377   }
 378
 379 }
 380
 381 //*********
 382
 383 // Evaluate weight for W decay distribution in t -> W b -> f fbar b.
 384
 385 double SigmaProcess::weightTopDecay( Event& process, int iResBeg,
 386   int iResEnd) {
 387
 388   // If not pair W d/s/b and mother t then return unit weight.
 389   if (iResEnd - iResBeg != 1) return 1.;
 390   int iW1  = iResBeg;
 391   int iB2  = iResBeg + 1;
 392   int idW1 = process[iW1].idAbs();
 393   int idB2 = process[iB2].idAbs();
 394   if (idW1 != 24) {
 395     swap(iW1, iB2);
 396     swap(idW1, idB2);
 397   }
 398   if (idW1 != 24 || (idB2 != 1 && idB2 != 3 && idB2 != 5)) return 1.;
 399   int iT   = process[iW1].mother1();
 400   if (iT <= 0 || process[iT].idAbs() != 6) return 1.;
 401
 402   // Find sign-matched order of W decay products.
 403   int iF    = process[iW1].daughter1();
 404   int iFbar = process[iW1].daughter2();
 405   if (iFbar - iF != 1) return 1.;
 406   if (process[iT].id() * process[iF].id() < 0) swap(iF, iFbar);
 407
 408   // Weight and maximum weight.
 409   double wt    = (process[iT].p() * process[iFbar].p())
 410                * (process[iF].p() * process[iB2].p());
 411   double wtMax = ( pow4(process[iT].m()) - pow4(process[iW1].m()) ) / 8.;
 412
 413   // Done.
 414   return wt / wtMax;
 415
 416 }
 417
 418
 419 //*********
 420
 421 // Evaluate weight for Z0/W+- decay distributions in H -> Z0/W+ Z0/W- -> 4f.
 422
 423 double SigmaProcess::weightHiggsDecay( Event& process, int iResBeg,
 424   int iResEnd) {
 425
 426   // If not pair Z0 Z0 or W+ W- then return unit weight.
 427   if (iResEnd - iResBeg != 1) return 1.;
 428   int iZW1  = iResBeg;
 429   int iZW2  = iResBeg + 1;
 430   int idZW1 = process[iZW1].id();
 431   int idZW2 = process[iZW2].id();
 432   if (idZW1 < 0) {
 433     swap(iZW1, iZW2);
 434     swap(idZW1, idZW2);
 435   }
 436   if ( (idZW1 != 23 || idZW2 != 23) && (idZW1 != 24 || idZW2 != -24) )
 437     return 1.;
 438
 439   // If mother is not Higgs then return unit weight.
 440   int iH  = process[iZW1].mother1();
 441   if (iH <= 0) return 1.;
 442   int idH = process[iH].id();
 443   if (idH != 25 && idH != 35 && idH !=36) return 1.;
 444
 445   // Parameters depend on Higgs type: H0(H_1), H^0(H_2) or A^0(H_3).
 446   int    higgsParity = higgsH1parity;
 447   double higgsEta    = higgsH1eta;
 448   if (idH == 35) {
 449     higgsParity      = higgsH2parity;
 450     higgsEta         = higgsH2eta;
 451   } else if (idH == 36) {
 452     higgsParity      = higgsA3parity;
 453     higgsEta         = higgsA3eta;
 454   }
 455
 456   // Option with isotropic decays.
 457   if (higgsParity == 0) return 1.;
 458
 459   // Maximum and initial weight.
 460   double wtMax = pow4(process[iH].m());
 461   double wt    = wtMax;
 462
 463   // Find sign-matched order of Z0/W+- decay products.
 464   int i3 = process[iZW1].daughter1();
 465   int i4 = process[iZW1].daughter2();
 466   if (process[i3].id() < 0) swap( i3, i4);
 467   int i5 = process[iZW2].daughter1();
 468   int i6 = process[iZW2].daughter2();
 469   if (process[i5].id() < 0) swap( i5, i6);
 470
 471   // Evaluate four-vector products and find masses..
 472   double p35  = 2. * process[i3].p() * process[i5].p();
 473   double p36  = 2. * process[i3].p() * process[i6].p();
 474   double p45  = 2. * process[i4].p() * process[i5].p();
 475   double p46  = 2. * process[i4].p() * process[i6].p();
 476   double p34  = 2. * process[i3].p() * process[i4].p();
 477   double p56  = 2. * process[i5].p() * process[i6].p();
 478   double mZW1 = process[iZW1].m();
 479   double mZW2 = process[iZW2].m();
 480
 481   // For mixed CP states need epsilon product and gauge boson masses.
 482   double epsilonProd = 0.;
 483   if (higgsParity == 3) {
 484     double p[4][4];
 485     for (int i = 0; i < 4; ++i) {
 486       int         ii = i3;
 487       if (i == 1) ii = i4;
 488       if (i == 2) ii = i5;
 489       if (i == 3) ii = i6;
 490       p[i][0] = process[ii].e();
 491       p[i][1] = process[ii].px();
 492       p[i][2] = process[ii].py();
 493       p[i][3] = process[ii].pz();
 494     }
 495     epsilonProd
 496       = p[0][0]*p[1][1]*p[2][2]*p[3][3] - p[0][0]*p[1][1]*p[2][3]*p[3][2]
 497       - p[0][0]*p[1][2]*p[2][1]*p[3][3] + p[0][0]*p[1][2]*p[2][3]*p[3][1]
 498       + p[0][0]*p[1][3]*p[2][1]*p[3][2] - p[0][0]*p[1][3]*p[2][2]*p[3][1]
 499       - p[0][1]*p[1][0]*p[2][2]*p[3][3] + p[0][1]*p[1][0]*p[2][3]*p[3][2]
 500       + p[0][1]*p[1][2]*p[2][0]*p[3][3] - p[0][1]*p[1][2]*p[2][3]*p[3][0]
 501       - p[0][1]*p[1][3]*p[2][0]*p[3][2] + p[0][1]*p[1][3]*p[2][2]*p[3][0]
 502       + p[0][2]*p[1][0]*p[2][1]*p[3][3] - p[0][2]*p[1][0]*p[2][3]*p[3][1]
 503       - p[0][2]*p[1][1]*p[2][0]*p[3][3] + p[0][2]*p[1][1]*p[2][3]*p[3][0]
 504       + p[0][2]*p[1][3]*p[2][0]*p[3][1] - p[0][2]*p[1][3]*p[2][1]*p[3][0]
 505       - p[0][3]*p[1][0]*p[2][1]*p[3][2] + p[0][3]*p[1][0]*p[2][2]*p[3][1]
 506       + p[0][3]*p[1][1]*p[2][0]*p[3][2] - p[0][3]*p[1][1]*p[2][2]*p[3][0]
 507       - p[0][3]*p[1][2]*p[2][0]*p[3][1] + p[0][3]*p[1][2]*p[2][1]*p[3][0];
 508   }
 509
 510   // Z0 Z0 decay: vector and axial couplings of two fermion pairs.
 511   if (idZW1 == 23) {
 512     double vf1 = CoupEW::vf(process[i3].idAbs());
 513     double af1 = CoupEW::af(process[i3].idAbs());
 514     double vf2 = CoupEW::vf(process[i5].idAbs());
 515     double af2 = CoupEW::af(process[i5].idAbs());
 516     double va12asym = 4. * vf1 * af1 * vf2 * af2
 517       / ( (vf1*vf1 + af1*af1) * (vf2*vf2 + af2*af2) );
 518     double etaMod = higgsEta / pow2( ParticleDataTable::m0(23) );
 519
 520     // Normal CP-even decay.
 521     if (higgsParity == 1) wt = 8. * (1. + va12asym) * p35 * p46
 522       + 8. * (1. - va12asym) * p36 * p45;
 523
 524     // CP-odd decay (normal for A0(H_3)).
 525     else if (higgsParity == 2) wt = ( pow2(p35 + p46)
 526       + pow2(p36 + p45) - 2. * p34 * p56
 527       - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
 528       + va12asym * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) )
 529       / (1. +  va12asym);
 530
 531     // Mixed CP states.
 532     else wt = 32. * ( 0.25 * ( (1. + va12asym) * p35 * p46
 533       + (1. - va12asym) * p36 * p45 ) - 0.5 * etaMod * epsilonProd
 534       * ( (1. + va12asym) * (p35 + p46) - (1. - va12asym) * (p36 + p45) )
 535       + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
 536       - 2. * pow2(p35 * p46 - p36 * p45)
 537       + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
 538       + va12asym * p34 * p56 * (p35 + p36 - p45 - p46)
 539       * (p35 + p45 - p36 - p46) ) ) / ( 1. + 2. * etaMod * mZW1 * mZW2
 540       + 2. * pow2(etaMod * mZW1 * mZW2) * (1. + va12asym) );
 541
 542   // W+ W- decay.
 543   } else if (idZW1 == 24) {
 544     double etaMod = higgsEta / pow2( ParticleDataTable::m0(24) );
 545
 546     // Normal CP-even decay.
 547     if (higgsParity == 1) wt = 16. * p35 * p46;
 548
 549     // CP-odd decay (normal for A0(H_3)).
 550     else if (higgsParity == 2) wt = 0.5 * ( pow2(p35 + p46)
 551       + pow2(p36 + p45) - 2. * p34 * p56
 552       - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
 553       + (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) );
 554
 555     // Mixed CP states.
 556     else wt = 32. * ( 0.25 * 2. * p35 * p46
 557       - 0.5 * etaMod * epsilonProd * 2. * (p35 + p46)
 558       + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
 559       - 2. * pow2(p35 * p46 - p36 * p45)
 560       + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
 561       + p34 * p56 * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) )
 562       / ( 1. * 2. * etaMod * mZW1 * mZW2 + 2. * pow2(etaMod * mZW1 * mZW2) );
 563   }
 564
 565   // Done.
 566   return wt / wtMax;
 567
 568 }
 569
 570 //**************************************************************************
 571
 572 // The Sigma1Process class.
 573 // Base class for resolved 2 -> 1 cross sections; derived from SigmaProcess.
 574
 575 //*********
 576
 577 // Input and complement kinematics for resolved 2 -> 1 process.
 578
 579 void Sigma1Process::store1Kin( double x1in, double x2in, double sHin) {
 580
 581   // Default value only sensible for these processes.
 582   swapTU = false;
 583
 584   // Incoming parton momentum fractions and sHat.
 585   x1Save = x1in;
 586   x2Save = x2in;
 587   sH     = sHin;
 588   mH     = sqrt(sH);
 589   sH2    = sH * sH;
 590
 591   // Different options for renormalization scale, but normally sHat.
 592   Q2RenSave = renormMultFac * sH;
 593   if (renormScale1 == 2) Q2RenSave = renormFixScale;
 594
 595   // Different options for factorization scale, but normally sHat.
 596   Q2FacSave = factorMultFac * sH;
 597   if (factorScale1 == 2) Q2RenSave = factorFixScale;
 598
 599   // Evaluate alpha_strong and alpha_EM.
 600   alpS   = alphaSPtr->alphaS(Q2RenSave);
 601   alpEM  = alphaEMPtr->alphaEM(Q2RenSave);
 602
 603 }
 604
 605 //**************************************************************************
 606
 607 // The Sigma2Process class.
 608 // Base class for resolved 2 -> 2 cross sections; derived from SigmaProcess.
 609
 610 //*********
 611
 612 // Input and complement kinematics for resolved 2 -> 2 process.
 613
 614 void Sigma2Process::store2Kin( double x1in, double x2in, double sHin,
 615   double tHin, double m3in, double m4in, double runBW3in, double runBW4in) {
 616
 617   // Default ordering of particles 3 and 4.
 618   swapTU   = false;
 619
 620   // Incoming parton momentum fractions.
 621   x1Save   = x1in;
 622   x2Save   = x2in;
 623
 624   // Incoming masses and their squares.
 625   bool masslessKin = (id3Mass() == 0) && (id4Mass() == 0);
 626   if (masslessKin) {
 627     m3     = 0.;
 628     m4     = 0.;
 629   } else {
 630     m3     = m3in;
 631     m4     = m4in;
 632   }
 633   mSave[3] = m3;
 634   mSave[4] = m4;
 635   s3       = m3 * m3;
 636   s4       = m4 * m4;
 637
 638   // Standard Mandelstam variables and their squares.
 639   sH       = sHin;
 640   tH       = tHin;
 641   uH       = (masslessKin) ? -(sH + tH) : s3 + s4 - (sH + tH);
 642   mH       = sqrt(sH);
 643   sH2      = sH * sH;
 644   tH2      = tH * tH;
 645   uH2      = uH * uH;
 646
 647   // The nominal Breit-Wigner factors with running width.
 648   runBW3   = runBW3in;
 649   runBW4   = runBW4in;
 650
 651   // Calculate squared transverse momentum.
 652   pT2 = (masslessKin) ?  tH * uH / sH : (tH * uH - s3 * s4) / sH;
 653
 654   // Special case: pick scale as if 2 -> 1 process in disguise.
 655   if (isSChannel()) {
 656
 657     // Different options for renormalization scale, but normally sHat.
 658     Q2RenSave = renormMultFac * sH;
 659     if (renormScale1 == 2) Q2RenSave = renormFixScale;
 660
 661     // Different options for factorization scale, but normally sHat.
 662     Q2FacSave = factorMultFac * sH;
 663     if (factorScale1 == 2) Q2RenSave = factorFixScale;
 664
 665   // Normal case with "true" 2 -> 2.
 666   } else {
 667
 668     // Different options for renormalization scale.
 669     if (masslessKin)            Q2RenSave = (renormScale2 < 4) ? pT2 : sH;
 670     else if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
 671     else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
 672     else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
 673     else                        Q2RenSave = sH;
 674     Q2RenSave *= renormMultFac;
 675     if      (renormScale2 == 5) Q2RenSave = renormFixScale;
 676
 677     // Different options for factorization scale.
 678     if (masslessKin)            Q2FacSave = (factorScale2 < 4) ? pT2 : sH;
 679     else if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
 680     else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
 681     else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
 682     else                        Q2FacSave = sH;
 683     Q2FacSave *= factorMultFac;
 684     if      (factorScale2 == 5) Q2FacSave = factorFixScale;
 685   }
 686
 687   // Evaluate alpha_strong and alpha_EM.
 688   alpS = alphaSPtr->alphaS(Q2RenSave);
 689   alpEM = alphaEMPtr->alphaEM(Q2RenSave);
 690
 691 }
 692
 693 //*********
 694
 695 // As above, special kinematics for multiple interactions.
 696
 697 void Sigma2Process::store2KinMI( double x1in, double x2in,
 698   double sHin, double tHin, double uHin, double alpSin, double alpEMin,
 699   bool needMasses, double m3in, double m4in) {
 700
 701   // Default ordering of particles 3 and 4.
 702   swapTU    = false;
 703
 704   // Incoming x values.
 705   x1Save    = x1in;
 706   x2Save    = x2in;
 707
 708   // Standard Mandelstam variables and their squares.
 709   sH        = sHin;
 710   tH        = tHin;
 711   uH        = uHin;
 712   mH        = sqrt(sH);
 713   sH2       = sH * sH;
 714   tH2       = tH * tH;
 715   uH2       = uH * uH;
 716
 717   // Strong and electroweak couplings.
 718   alpS      = alpSin;
 719   alpEM     = alpEMin;
 720
 721   // Assume vanishing masses. (Will be modified in final kinematics.)
 722   m3        = 0.;
 723   s3        = 0.;
 724   m4        = 0.;
 725   s4        = 0.;
 726   sHBeta    = sH;
 727
 728   // Scattering angle.
 729   cosTheta  = (tH - uH) / sH;
 730   sinTheta  = 2. * sqrtpos( tH * uH ) / sH;
 731
 732   // In some cases must use masses and redefine meaning of tHat and uHat.
 733   if (needMasses) {
 734     m3      = m3in;
 735     s3      = m3 * m3;
 736     m4      = m4in;
 737     s4      = m4 * m4;
 738     sHMass  = sH - s3 - s4;
 739     sHBeta  = sqrtpos(sHMass*sHMass - 4. * s3 * s4);
 740     tH      = -0.5 * (sHMass - sHBeta * cosTheta);
 741     uH      = -0.5 * (sHMass + sHBeta * cosTheta);
 742     tH2     = tH * tH;
 743     uH2     = uH * uH;
 744   }
 745
 746   // pT2 with masses (at this stage) included.
 747   pT2Mass   = 0.25 * sHBeta * pow2(sinTheta);
 748
 749 }
 750
 751 //*********
 752
 753 // Perform kinematics for a Multiple Interaction.
 754
 755 bool Sigma2Process::final2KinMI() {
 756
 757   // Have to set flavours and colours.
 758   setIdColAcol();
 759
 760   // Check that masses of outgoing particles not too big.
 761   m3           = ParticleDataTable::m0(idSave[3]);
 762   m4           = ParticleDataTable::m0(idSave[4]);
 763   mH           = sqrt(sH);
 764   if (m3 + m4 + MASSMARGIN > mH) return false;
 765   s3           = m3 * m3;
 766   s4           = m4 * m4;
 767
 768   // Do kinematics of the decay.
 769   double eIn   = 0.5 * mH;
 770   double e3    = 0.5 * (sH + s3 - s4) / mH;
 771   double e4    = 0.5 * (sH + s4 - s3) / mH;
 772   double pAbs  = sqrtpos( e3*e3 - s3 );
 773   phi          = 2. * M_PI * Rndm::flat();
 774   double pZ    = pAbs * cosTheta;
 775   double pX    = pAbs * sinTheta * sin(phi);
 776   double pY    = pAbs * sinTheta * cos(phi);
 777   double scale = eIn * sinTheta;
 778
 779   // Fill particle info.
 780   parton[1] = Particle( idSave[1], -31, 0, 0, 3, 4, colSave[1], acolSave[1],
 781     0., 0., eIn, eIn, 0., scale);
 782   parton[2] = Particle( idSave[2], -31, 0, 0, 3, 4, colSave[2], acolSave[2],
 783     0., 0., -eIn, eIn, 0., scale);
 784   parton[3] = Particle( idSave[3],  33, 1, 2, 0, 0, colSave[3], acolSave[3],
 785     pX, pY, pZ, e3, m3, scale);
 786   parton[4] = Particle( idSave[4],  33, 1, 2, 0, 0, colSave[4], acolSave[4],
 787     -pX, -pY, -pZ, e4, m4, scale);
 788
 789   // Boost particles from subprocess rest frame to event rest frame.
 790   double betaZ = (x1Save - x2Save) / (x1Save + x2Save);
 791   for (int i = 1; i <= 4; ++i) parton[i].bst(0., 0., betaZ);
 792
 793   // Done.
 794   return true;
 795
 796 }
 797
 798 //**************************************************************************
 799
 800 // The Sigma3Process class.
 801 // Base class for resolved 2 -> 3 cross sections; derived from SigmaProcess.
 802
 803 //*********
 804
 805 // Input and complement kinematics for resolved 2 -> 3 process.
 806
 807 void Sigma3Process::store3Kin( double x1in, double x2in, double sHin,
 808   Vec4 p3cmIn, Vec4 p4cmIn, Vec4 p5cmIn, double m3in, double m4in,
 809   double m5in, double runBW3in, double runBW4in, double runBW5in) {
 810
 811   // Default ordering of particles 3 and 4 - not relevant here.
 812   swapTU   = false;
 813
 814   // Incoming parton momentum fractions.
 815   x1Save   = x1in;
 816   x2Save   = x2in;
 817
 818   // Incoming masses and their squares.
 819   if (id3Mass() == 0 && id4Mass() == 0 && id5Mass() == 0) {
 820     m3     = 0.;
 821     m4     = 0.;
 822     m5     = 0.;
 823   } else {
 824     m3     = m3in;
 825     m4     = m4in;
 826     m5     = m5in;
 827   }
 828   mSave[3] = m3;
 829   mSave[4] = m4;
 830   mSave[5] = m5;
 831   s3       = m3 * m3;
 832   s4       = m4 * m4;
 833   s5       = m5 * m5;
 834
 835   // Standard Mandelstam variables and four-momenta in rest frame.
 836   sH       = sHin;
 837   mH       = sqrt(sH);
 838   p3cm     = p3cmIn;
 839   p4cm     = p4cmIn;
 840   p5cm     = p5cmIn;
 841
 842   // The nominal Breit-Wigner factors with running width.
 843   runBW3   = runBW3in;
 844   runBW4   = runBW4in;
 845   runBW5   = runBW5in;
 846
 847   // Special case: pick scale as if 2 -> 1 process in disguise.
 848   if (isSChannel()) {
 849
 850     // Different options for renormalization scale, but normally sHat.
 851     Q2RenSave = renormMultFac * sH;
 852     if (renormScale1 == 2) Q2RenSave = renormFixScale;
 853
 854     // Different options for factorization scale, but normally sHat.
 855     Q2FacSave = factorMultFac * sH;
 856     if (factorScale1 == 2) Q2RenSave = factorFixScale;
 857
 858   // "Normal" 2 -> 3 processes, i.e. not vector boson fusion.
 859   } else if ( idTchan1() != 23 && idTchan1() != 24 && idTchan2() != 23
 860     && idTchan1() != 24 ) {
 861     double mT3S = s3 + p3cm.pT2();
 862     double mT4S = s4 + p4cm.pT2();
 863     double mT5S = s5 + p5cm.pT2();
 864
 865     // Different options for renormalization scale.
 866     if      (renormScale3 == 1) Q2RenSave = min( mT3S, min(mT4S, mT5S) );
 867     else if (renormScale3 == 2) Q2RenSave = sqrt( mT3S * mT4S * mT5S
 868       / max( mT3S, max(mT4S, mT5S) ) );
 869     else if (renormScale3 == 3) Q2RenSave = pow( mT3S * mT4S * mT5S,
 870                                             1./3. );
 871     else if (renormScale3 == 4) Q2RenSave = (mT3S * mT4S * mT5S) / 3.;
 872     else                        Q2RenSave = sH;
 873     Q2RenSave *= renormMultFac;
 874     if      (renormScale3 == 6) Q2RenSave = renormFixScale;
 875
 876     // Different options for factorization scale.
 877     if      (factorScale3 == 1) Q2FacSave = min( mT3S, min(mT4S, mT5S) );
 878     else if (factorScale3 == 2) Q2FacSave = sqrt( mT3S * mT4S * mT5S
 879       / max( mT3S, max(mT4S, mT5S) ) );
 880     else if (factorScale3 == 3) Q2FacSave = pow( mT3S * mT4S * mT5S,
 881                                             1./3. );
 882     else if (factorScale3 == 4) Q2FacSave = (mT3S * mT4S * mT5S) / 3.;
 883     else                        Q2FacSave = sH;
 884     Q2RenSave *= factorMultFac;
 885     if      (factorScale3 == 6) Q2FacSave = factorFixScale;
 886
 887   // Vector boson fusion 2 -> 3 processes; recoils in positions 4 and 5.
 888   } else {
 889     double sV4   = pow2( ParticleDataTable::m0(idTchan1()) );
 890     double sV5   = pow2( ParticleDataTable::m0(idTchan2()) );
 891     double mT3S  = s3  + p3cm.pT2();
 892     double mTV4S = sV4 + p4cm.pT2();
 893     double mTV5S = sV5 + p5cm.pT2();
 894
 895     // Different options for renormalization scale.
 896     if      (renormScale3VV == 1) Q2RenSave = max( sV4, sV5);
 897     else if (renormScale3VV == 2) Q2RenSave = sqrt( mTV4S * mTV5S );
 898     else if (renormScale3VV == 3) Q2RenSave = pow( mT3S * mTV4S * mTV5S,
 899                                               1./3. );
 900     else if (renormScale3VV == 4) Q2RenSave = (mT3S * mTV4S * mTV5S) / 3.;
 901     else                          Q2RenSave = sH;
 902     Q2RenSave *= renormMultFac;
 903     if      (renormScale3VV == 6) Q2RenSave = renormFixScale;
 904
 905     // Different options for factorization scale.
 906     if      (factorScale3VV == 1) Q2FacSave = max( sV4, sV5);
 907     else if (factorScale3VV == 2) Q2FacSave = sqrt( mTV4S * mTV5S );
 908     else if (factorScale3VV == 3) Q2FacSave = pow( mT3S * mTV4S * mTV5S,
 909                                               1./3. );
 910     else if (factorScale3VV == 4) Q2FacSave = (mT3S * mTV4S * mTV5S) / 3.;
 911     else                          Q2FacSave = sH;
 912     Q2RenSave *= factorMultFac;
 913     if      (factorScale3VV == 6) Q2FacSave = factorFixScale;
 914   }
 915
 916   // Evaluate alpha_strong and alpha_EM.
 917   alpS = alphaSPtr->alphaS(Q2RenSave);
 918   alpEM = alphaEMPtr->alphaEM(Q2RenSave);
 919
 920 }
 921
 922 //**************************************************************************
 923
 924 // The SigmaLHAProcess class.
 925 // Wrapper for Les Houches Accord external input; derived from SigmaProcess.
 926 // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
 927
 928 //*********
 929
 930 // Obtain number of final-state partons from LHA object.
 931
 932 int SigmaLHAProcess::nFinal() const {
 933
 934   // At initialization size unknown, so return 0.
 935   if (lhaUpPtr->sizePart() <= 0) return 0;
 936
 937   // Sum up all particles that has first mother = 1.
 938   int nFin = 0;
 939   for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
 940     if (lhaUpPtr->mother1(i) == 1) ++nFin;
 941   return nFin;
 942
 943 }
 944
 945 //**************************************************************************
 946
 947 } // end namespace Pythia8