+++ /dev/null
-// SigmaProcess.cc is a part of the PYTHIA event generator.
-// Copyright (C) 2012 Torbjorn Sjostrand.
-// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
-// Please respect the MCnet Guidelines, see GUIDELINES for details.
-
-// Function definitions (not found in the header) for the
-// SigmaProcess class, and classes derived from it.
-
-#include "SigmaProcess.h"
-
-namespace Pythia8 {
-
-//==========================================================================
-
-// The SigmaProcess class.
-// Base class for cross sections.
-
-//--------------------------------------------------------------------------
-
-// Constants: could be changed here if desired, but normally should not.
-// These are of technical nature, as described for each.
-
-// Conversion of GeV^{-2} to mb for cross section.
-const double SigmaProcess::CONVERT2MB = 0.389380;
-
-// The sum of outgoing masses must not be too close to the cm energy.
-const double SigmaProcess::MASSMARGIN = 0.1;
-
-// Parameters of momentum rescaling procedure: maximally allowed
-// relative energy error and number of iterations.
-const double SigmaProcess::COMPRELERR = 1e-10;
-const int SigmaProcess::NCOMPSTEP = 10;
-
-//--------------------------------------------------------------------------
-
-// Perform simple initialization and store pointers.
-
-void SigmaProcess::init(Info* infoPtrIn, Settings* settingsPtrIn,
- ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, BeamParticle* beamAPtrIn,
- BeamParticle* beamBPtrIn, Couplings* couplingsPtrIn,
- SigmaTotal* sigmaTotPtrIn, SusyLesHouches* slhaPtrIn) {
-
- // Store pointers.
- infoPtr = infoPtrIn;
- settingsPtr = settingsPtrIn;
- particleDataPtr = particleDataPtrIn;
- rndmPtr = rndmPtrIn;
- beamAPtr = beamAPtrIn;
- beamBPtr = beamBPtrIn;
- couplingsPtr = couplingsPtrIn;
- sigmaTotPtr = sigmaTotPtrIn;
- slhaPtr = slhaPtrIn;
-
- // Read out some properties of beams to allow shorthand.
- idA = (beamAPtr != 0) ? beamAPtr->id() : 0;
- idB = (beamBPtr != 0) ? beamBPtr->id() : 0;
- mA = (beamAPtr != 0) ? beamAPtr->m() : 0.;
- mB = (beamBPtr != 0) ? beamBPtr->m() : 0.;
- isLeptonA = (beamAPtr != 0) ? beamAPtr->isLepton() : false;
- isLeptonB = (beamBPtr != 0) ? beamBPtr->isLepton() : false;
- hasLeptonBeams = isLeptonA || isLeptonB;
-
- // K factor, multiplying resolved processes. (But not here for MPI.)
- Kfactor = settingsPtr->parm("SigmaProcess:Kfactor");
-
- // Maximum incoming quark flavour.
- nQuarkIn = settingsPtr->mode("PDFinProcess:nQuarkIn");
-
- // Medium heavy fermion masses set massless or not in ME expressions.
- mcME = (settingsPtr->flag("SigmaProcess:cMassiveME"))
- ? particleDataPtr->m0(4) : 0.;
- mbME = (settingsPtr->flag("SigmaProcess:bMassiveME"))
- ? particleDataPtr->m0(5) : 0.;
- mmuME = (settingsPtr->flag("SigmaProcess:muMassiveME"))
- ? particleDataPtr->m0(13) : 0.;
- mtauME = (settingsPtr->flag("SigmaProcess:tauMassiveME"))
- ? particleDataPtr->m0(15) : 0.;
-
- // Renormalization scale choice.
- renormScale1 = settingsPtr->mode("SigmaProcess:renormScale1");
- renormScale2 = settingsPtr->mode("SigmaProcess:renormScale2");
- renormScale3 = settingsPtr->mode("SigmaProcess:renormScale3");
- renormScale3VV = settingsPtr->mode("SigmaProcess:renormScale3VV");
- renormMultFac = settingsPtr->parm("SigmaProcess:renormMultFac");
- renormFixScale = settingsPtr->parm("SigmaProcess:renormFixScale");
-
- // Factorization scale choice.
- factorScale1 = settingsPtr->mode("SigmaProcess:factorScale1");
- factorScale2 = settingsPtr->mode("SigmaProcess:factorScale2");
- factorScale3 = settingsPtr->mode("SigmaProcess:factorScale3");
- factorScale3VV = settingsPtr->mode("SigmaProcess:factorScale3VV");
- factorMultFac = settingsPtr->parm("SigmaProcess:factorMultFac");
- factorFixScale = settingsPtr->parm("SigmaProcess:factorFixScale");
-
- // CP violation parameters for the BSM Higgs sector.
- higgsH1parity = settingsPtr->mode("HiggsH1:parity");
- higgsH1eta = settingsPtr->parm("HiggsH1:etaParity");
- higgsH2parity = settingsPtr->mode("HiggsH2:parity");
- higgsH2eta = settingsPtr->parm("HiggsH2:etaParity");
- higgsA3parity = settingsPtr->mode("HiggsA3:parity");
- higgsA3eta = settingsPtr->parm("HiggsA3:etaParity");
-
- // If BSM not switched on then H1 should have SM properties.
- if (!settingsPtr->flag("Higgs:useBSM")){
- higgsH1parity = 1;
- higgsH1eta = 0.;
- }
-
-}
-
-//--------------------------------------------------------------------------
-
-// Set up allowed flux of incoming partons.
-// addBeam: set up PDF's that need to be evaluated for the two beams.
-// addPair: set up pairs of incoming partons from the two beams.
-
-bool SigmaProcess::initFlux() {
-
- // Reset arrays (in case of several init's in same run).
- inBeamA.clear();
- inBeamB.clear();
- inPair.clear();
-
- // Read in process-specific channel information.
- string fluxType = inFlux();
-
- // Case with g g incoming state.
- if (fluxType == "gg") {
- addBeamA(21);
- addBeamB(21);
- addPair(21, 21);
- }
-
- // Case with q g incoming state.
- else if (fluxType == "qg") {
- for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
- int idNow = (i == 0) ? 21 : i;
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addPair(idNow, 21);
- addPair(21, idNow);
- }
- }
-
- // Case with q q', q qbar' or qbar qbar' incoming state.
- else if (fluxType == "qq") {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
- if (id1Now != 0)
- for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
- if (id2Now != 0)
- addPair(id1Now, id2Now);
- }
-
- // Case with q qbar incoming state.
- else if (fluxType == "qqbarSame") {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0)
- addPair(idNow, -idNow);
- }
-
- // Case with f f', f fbar', fbar fbar' incoming state.
- else if (fluxType == "ff") {
- // If beams are leptons then they are also the colliding partons.
- if ( isLeptonA && isLeptonB ) {
- addBeamA(idA);
- addBeamB(idB);
- addPair(idA, idB);
- // First beam is lepton and second is hadron.
- } else if ( isLeptonA ) {
- addBeamA(idA);
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamB(idNow);
- addPair(idA, idNow);
- }
- // First beam is hadron and second is lepton.
- } else if ( isLeptonB ) {
- addBeamB(idB);
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addPair(idNow, idB);
- }
- // Hadron beams gives quarks.
- } else {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
- if (id1Now != 0)
- for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
- if (id2Now != 0)
- addPair(id1Now, id2Now);
- }
- }
-
- // Case with f fbar incoming state.
- else if (fluxType == "ffbarSame") {
- // If beams are antiparticle pair and leptons then also colliding partons.
- if ( idA + idB == 0 && isLeptonA ) {
- addBeamA(idA);
- addBeamB(idB);
- addPair(idA, idB);
- // Else assume both to be hadrons, for better or worse.
- } else {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0)
- addPair(idNow, -idNow);
- }
- }
-
- // Case with f fbar' charged(+-1) incoming state.
- else if (fluxType == "ffbarChg") {
- // If beams are leptons then also colliding partons.
- if ( isLeptonA && isLeptonB && abs( particleDataPtr->chargeType(idA)
- + particleDataPtr->chargeType(idB) ) == 3 ) {
- addBeamA(idA);
- addBeamB(idB);
- addPair(idA, idB);
- // Hadron beams gives quarks.
- } else {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
- if (id1Now != 0)
- for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
- if (id2Now != 0 && id1Now * id2Now < 0
- && (abs(id1Now) + abs(id2Now))%2 == 1) addPair(id1Now, id2Now);
- }
- }
-
- // Case with f fbar' generic incoming state.
- else if (fluxType == "ffbar") {
- // If beams are leptons then also colliding partons.
- if (isLeptonA && isLeptonB && idA * idB < 0) {
- addBeamA(idA);
- addBeamB(idB);
- addPair(idA, idB);
- // Hadron beams gives quarks.
- } else {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addBeamB(idNow);
- }
- for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
- if (id1Now != 0)
- for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
- if (id2Now != 0 && id1Now * id2Now < 0)
- addPair(id1Now, id2Now);
- }
- }
-
- // Case with f gamma incoming state.
- else if (fluxType == "fgm") {
- // Fermion from incoming side A.
- if ( isLeptonA ) {
- addBeamA(idA);
- addPair(idA, 22);
- } else {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamA(idNow);
- addPair(idNow, 22);
- }
- }
- // Fermion from incoming side B.
- if ( isLeptonB ) {
- addBeamB( idB);
- addPair(22, idB);
- } else {
- for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
- if (idNow != 0) {
- addBeamB(idNow);
- addPair(22, idNow);
- }
- }
- // Photons in the beams.
- addBeamA(22);
- addBeamB(22);
- }
-
- // Case with gamma gamma incoming state.
- else if (fluxType == "ggm") {
- addBeamA(21);
- addBeamA(22);
- addBeamB(21);
- addBeamB(22);
- addPair(21, 22);
- addPair(22, 21);
- }
-
- // Case with gamma gamma incoming state.
- else if (fluxType == "gmgm") {
- addBeamA(22);
- addBeamB(22);
- addPair(22, 22);
- }
-
- // Unrecognized fluxType is bad sign. Else done.
- else {
- infoPtr->errorMsg("Error in SigmaProcess::initFlux: "
- "unrecognized inFlux type", fluxType);
- return false;
- }
- return true;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Convolute matrix-element expression(s) with parton flux and K factor.
-
-double SigmaProcess::sigmaPDF() {
-
- // Evaluate and store the required parton densities.
- for (int j = 0; j < sizeBeamA(); ++j)
- inBeamA[j].pdf = beamAPtr->xfHard( inBeamA[j].id, x1Save, Q2FacSave);
- for (int j = 0; j < sizeBeamB(); ++j)
- inBeamB[j].pdf = beamBPtr->xfHard( inBeamB[j].id, x2Save, Q2FacSave);
-
- // Loop over allowed incoming channels.
- sigmaSumSave = 0.;
- for (int i = 0; i < sizePair(); ++i) {
-
- // Evaluate hard-scattering cross section. Include K factor.
- inPair[i].pdfSigma = Kfactor
- * sigmaHatWrap(inPair[i].idA, inPair[i].idB);
-
- // Multiply by respective parton densities.
- for (int j = 0; j < sizeBeamA(); ++j)
- if (inPair[i].idA == inBeamA[j].id) {
- inPair[i].pdfA = inBeamA[j].pdf;
- inPair[i].pdfSigma *= inBeamA[j].pdf;
- break;
- }
- for (int j = 0; j < sizeBeamB(); ++j)
- if (inPair[i].idB == inBeamB[j].id) {
- inPair[i].pdfB = inBeamB[j].pdf;
- inPair[i].pdfSigma *= inBeamB[j].pdf;
- break;
- }
-
- // Sum for all channels.
- sigmaSumSave += inPair[i].pdfSigma;
- }
-
- // Done.
- return sigmaSumSave;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Select incoming parton channel and extract parton densities (resolved).
-
-void SigmaProcess::pickInState(int id1in, int id2in) {
-
- // Multiparton interactions: partons already selected.
- if (id1in != 0 && id2in != 0) {
- id1 = id1in;
- id2 = id2in;
- }
-
- // Pick channel. Extract channel flavours and pdf's.
- double sigmaRand = sigmaSumSave * rndmPtr->flat();
- for (int i = 0; i < sizePair(); ++i) {
- sigmaRand -= inPair[i].pdfSigma;
- if (sigmaRand <= 0.) {
- id1 = inPair[i].idA;
- id2 = inPair[i].idB;
- pdf1Save = inPair[i].pdfA;
- pdf2Save = inPair[i].pdfB;
- break;
- }
- }
-
-}
-
-//--------------------------------------------------------------------------
-
-// Calculate incoming modified masses and four-vectors for matrix elements.
-
-bool SigmaProcess::setupForMEin() {
-
- // Initially assume it will work out to set up modified kinematics.
- bool allowME = true;
-
- // Correct incoming c, b, mu and tau to be massive or not.
- mME[0] = 0.;
- int id1Tmp = abs(id1);
- if (id1Tmp == 4) mME[0] = mcME;
- if (id1Tmp == 5) mME[0] = mbME;
- if (id1Tmp == 13) mME[0] = mmuME;
- if (id1Tmp == 15) mME[0] = mtauME;
- mME[1] = 0.;
- int id2Tmp = abs(id2);
- if (id2Tmp == 4) mME[1] = mcME;
- if (id2Tmp == 5) mME[1] = mbME;
- if (id2Tmp == 13) mME[1] = mmuME;
- if (id2Tmp == 15) mME[1] = mtauME;
-
- // If kinematically impossible return to massless case, but set error.
- if (mME[0] + mME[1] >= mH) {
- mME[0] = 0.;
- mME[1] = 0.;
- allowME = false;
- }
-
- // Do incoming two-body kinematics for massless or massive cases.
- if (mME[0] == 0. && mME[1] == 0.) {
- pME[0] = 0.5 * mH * Vec4( 0., 0., 1., 1.);
- pME[1] = 0.5 * mH * Vec4( 0., 0., -1., 1.);
- } else {
- double e0 = 0.5 * (mH * mH + mME[0] * mME[0] - mME[1] * mME[1]) / mH;
- double pz0 = sqrtpos(e0 * e0 - mME[0] * mME[0]);
- pME[0] = Vec4( 0., 0., pz0, e0);
- pME[1] = Vec4( 0., 0., -pz0, mH - e0);
- }
-
- // Done.
- return allowME;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Evaluate weight for W decay distribution in t -> W b -> f fbar b.
-
-double SigmaProcess::weightTopDecay( Event& process, int iResBeg,
- int iResEnd) {
-
- // If not pair W d/s/b and mother t then return unit weight.
- if (iResEnd - iResBeg != 1) return 1.;
- int iW1 = iResBeg;
- int iB2 = iResBeg + 1;
- int idW1 = process[iW1].idAbs();
- int idB2 = process[iB2].idAbs();
- if (idW1 != 24) {
- swap(iW1, iB2);
- swap(idW1, idB2);
- }
- if (idW1 != 24 || (idB2 != 1 && idB2 != 3 && idB2 != 5)) return 1.;
- int iT = process[iW1].mother1();
- if (iT <= 0 || process[iT].idAbs() != 6) return 1.;
-
- // Find sign-matched order of W decay products.
- int iF = process[iW1].daughter1();
- int iFbar = process[iW1].daughter2();
- if (iFbar - iF != 1) return 1.;
- if (process[iT].id() * process[iF].id() < 0) swap(iF, iFbar);
-
- // Weight and maximum weight.
- double wt = (process[iT].p() * process[iFbar].p())
- * (process[iF].p() * process[iB2].p());
- double wtMax = ( pow4(process[iT].m()) - pow4(process[iW1].m()) ) / 8.;
-
- // Done.
- return wt / wtMax;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Evaluate weight for Z0/W+- decay distributions in H -> Z0/W+ Z0/W- -> 4f
-// and H -> gamma Z0 -> gamma f fbar.
-
-double SigmaProcess::weightHiggsDecay( Event& process, int iResBeg,
- int iResEnd) {
-
- // If not pair Z0 Z0, W+ W- or gamma Z0 then return unit weight.
- if (iResEnd - iResBeg != 1) return 1.;
- int iZW1 = iResBeg;
- int iZW2 = iResBeg + 1;
- int idZW1 = process[iZW1].id();
- int idZW2 = process[iZW2].id();
- if (idZW1 < 0 || idZW2 == 22) {
- swap(iZW1, iZW2);
- swap(idZW1, idZW2);
- }
- if ( (idZW1 != 23 || idZW2 != 23) && (idZW1 != 24 || idZW2 != -24)
- && (idZW1 != 22 || idZW2 != 23) ) return 1.;
-
- // If mother is not Higgs then return unit weight.
- int iH = process[iZW1].mother1();
- if (iH <= 0) return 1.;
- int idH = process[iH].id();
- if (idH != 25 && idH != 35 && idH !=36) return 1.;
-
- // H -> gamma Z0 -> gamma f fbar is 1 + cos^2(theta) in Z rest frame.
- if (idZW1 == 22) {
- int i5 = process[iZW2].daughter1();
- int i6 = process[iZW2].daughter2();
- double pgmZ = process[iZW1].p() * process[iZW2].p();
- double pgm5 = process[iZW1].p() * process[i5].p();
- double pgm6 = process[iZW1].p() * process[i6].p();
- return (pow2(pgm5) + pow2(pgm6)) / pow2(pgmZ);
- }
-
- // Parameters depend on Higgs type: H0(H_1), H^0(H_2) or A^0(H_3).
- int higgsParity = higgsH1parity;
- double higgsEta = higgsH1eta;
- if (idH == 35) {
- higgsParity = higgsH2parity;
- higgsEta = higgsH2eta;
- } else if (idH == 36) {
- higgsParity = higgsA3parity;
- higgsEta = higgsA3eta;
- }
-
- // Option with isotropic decays.
- if (higgsParity == 0) return 1.;
-
- // Maximum and initial weight.
- double wtMax = pow4(process[iH].m());
- double wt = wtMax;
-
- // Find sign-matched order of Z0/W+- decay products.
- int i3 = process[iZW1].daughter1();
- int i4 = process[iZW1].daughter2();
- if (process[i3].id() < 0) swap( i3, i4);
- int i5 = process[iZW2].daughter1();
- int i6 = process[iZW2].daughter2();
- if (process[i5].id() < 0) swap( i5, i6);
-
- // Evaluate four-vector products and find masses..
- double p35 = 2. * process[i3].p() * process[i5].p();
- double p36 = 2. * process[i3].p() * process[i6].p();
- double p45 = 2. * process[i4].p() * process[i5].p();
- double p46 = 2. * process[i4].p() * process[i6].p();
- double p34 = 2. * process[i3].p() * process[i4].p();
- double p56 = 2. * process[i5].p() * process[i6].p();
- double mZW1 = process[iZW1].m();
- double mZW2 = process[iZW2].m();
-
- // For mixed CP states need epsilon product and gauge boson masses.
- double epsilonProd = 0.;
- if (higgsParity == 3) {
- double p[4][4];
- for (int i = 0; i < 4; ++i) {
- int ii = i3;
- if (i == 1) ii = i4;
- if (i == 2) ii = i5;
- if (i == 3) ii = i6;
- p[i][0] = process[ii].e();
- p[i][1] = process[ii].px();
- p[i][2] = process[ii].py();
- p[i][3] = process[ii].pz();
- }
- epsilonProd
- = p[0][0]*p[1][1]*p[2][2]*p[3][3] - p[0][0]*p[1][1]*p[2][3]*p[3][2]
- - p[0][0]*p[1][2]*p[2][1]*p[3][3] + p[0][0]*p[1][2]*p[2][3]*p[3][1]
- + p[0][0]*p[1][3]*p[2][1]*p[3][2] - p[0][0]*p[1][3]*p[2][2]*p[3][1]
- - p[0][1]*p[1][0]*p[2][2]*p[3][3] + p[0][1]*p[1][0]*p[2][3]*p[3][2]
- + p[0][1]*p[1][2]*p[2][0]*p[3][3] - p[0][1]*p[1][2]*p[2][3]*p[3][0]
- - p[0][1]*p[1][3]*p[2][0]*p[3][2] + p[0][1]*p[1][3]*p[2][2]*p[3][0]
- + p[0][2]*p[1][0]*p[2][1]*p[3][3] - p[0][2]*p[1][0]*p[2][3]*p[3][1]
- - p[0][2]*p[1][1]*p[2][0]*p[3][3] + p[0][2]*p[1][1]*p[2][3]*p[3][0]
- + p[0][2]*p[1][3]*p[2][0]*p[3][1] - p[0][2]*p[1][3]*p[2][1]*p[3][0]
- - p[0][3]*p[1][0]*p[2][1]*p[3][2] + p[0][3]*p[1][0]*p[2][2]*p[3][1]
- + p[0][3]*p[1][1]*p[2][0]*p[3][2] - p[0][3]*p[1][1]*p[2][2]*p[3][0]
- - p[0][3]*p[1][2]*p[2][0]*p[3][1] + p[0][3]*p[1][2]*p[2][1]*p[3][0];
- }
-
- // Z0 Z0 decay: vector and axial couplings of two fermion pairs.
- if (idZW1 == 23) {
- double vf1 = couplingsPtr->vf(process[i3].idAbs());
- double af1 = couplingsPtr->af(process[i3].idAbs());
- double vf2 = couplingsPtr->vf(process[i5].idAbs());
- double af2 = couplingsPtr->af(process[i5].idAbs());
- double va12asym = 4. * vf1 * af1 * vf2 * af2
- / ( (vf1*vf1 + af1*af1) * (vf2*vf2 + af2*af2) );
- double etaMod = higgsEta / pow2( particleDataPtr->m0(23) );
-
- // Normal CP-even decay.
- if (higgsParity == 1) wt = 8. * (1. + va12asym) * p35 * p46
- + 8. * (1. - va12asym) * p36 * p45;
-
- // CP-odd decay (normal for A0(H_3)).
- else if (higgsParity == 2) wt = ( pow2(p35 + p46)
- + pow2(p36 + p45) - 2. * p34 * p56
- - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
- + va12asym * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) )
- / (1. + va12asym);
-
- // Mixed CP states.
- else wt = 32. * ( 0.25 * ( (1. + va12asym) * p35 * p46
- + (1. - va12asym) * p36 * p45 ) - 0.5 * etaMod * epsilonProd
- * ( (1. + va12asym) * (p35 + p46) - (1. - va12asym) * (p36 + p45) )
- + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
- - 2. * pow2(p35 * p46 - p36 * p45)
- + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
- + va12asym * p34 * p56 * (p35 + p36 - p45 - p46)
- * (p35 + p45 - p36 - p46) ) ) / ( 1. + 2. * etaMod * mZW1 * mZW2
- + 2. * pow2(etaMod * mZW1 * mZW2) * (1. + va12asym) );
-
- // W+ W- decay.
- } else if (idZW1 == 24) {
- double etaMod = higgsEta / pow2( particleDataPtr->m0(24) );
-
- // Normal CP-even decay.
- if (higgsParity == 1) wt = 16. * p35 * p46;
-
- // CP-odd decay (normal for A0(H_3)).
- else if (higgsParity == 2) wt = 0.5 * ( pow2(p35 + p46)
- + pow2(p36 + p45) - 2. * p34 * p56
- - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
- + (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) );
-
- // Mixed CP states.
- else wt = 32. * ( 0.25 * 2. * p35 * p46
- - 0.5 * etaMod * epsilonProd * 2. * (p35 + p46)
- + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
- - 2. * pow2(p35 * p46 - p36 * p45)
- + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
- + p34 * p56 * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) )
- / ( 1. * 2. * etaMod * mZW1 * mZW2 + 2. * pow2(etaMod * mZW1 * mZW2) );
- }
-
- // Done.
- return wt / wtMax;
-
-}
-
-//==========================================================================
-
-// The Sigma1Process class.
-// Base class for resolved 2 -> 1 cross sections; derived from SigmaProcess.
-
-//--------------------------------------------------------------------------
-
-// Wrapper to sigmaHat, to (a) store current incoming flavours,
-// (b) convert from GeV^-2 to mb where required, and
-// (c) convert from |M|^2 to d(sigmaHat)/d(tHat) where required.
-
-double Sigma1Process::sigmaHatWrap(int id1in, int id2in) {
-
- id1 = id1in;
- id2 = id2in;
- double sigmaTmp = sigmaHat();
- if (convertM2()) {
- sigmaTmp /= 2. * sH;
- // Convert 2 * pi * delta(p^2 - m^2) to Breit-Wigner with same area.
- int idTmp = resonanceA();
- double mTmp = particleDataPtr->m0(idTmp);
- double GamTmp = particleDataPtr->mWidth(idTmp);
- sigmaTmp *= 2. * mTmp * GamTmp / ( pow2(sH - mTmp * mTmp)
- + pow2(mTmp * GamTmp) );
- }
- if (convert2mb()) sigmaTmp *= CONVERT2MB;
- return sigmaTmp;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Input and complement kinematics for resolved 2 -> 1 process.
-
-void Sigma1Process::store1Kin( double x1in, double x2in, double sHin) {
-
- // Default value only sensible for these processes.
- swapTU = false;
-
- // Incoming parton momentum fractions and sHat.
- x1Save = x1in;
- x2Save = x2in;
- sH = sHin;
- mH = sqrt(sH);
- sH2 = sH * sH;
-
- // Different options for renormalization scale, but normally sHat.
- Q2RenSave = renormMultFac * sH;
- if (renormScale1 == 2) Q2RenSave = renormFixScale;
-
- // Different options for factorization scale, but normally sHat.
- Q2FacSave = factorMultFac * sH;
- if (factorScale1 == 2) Q2FacSave = factorFixScale;
-
- // Evaluate alpha_strong and alpha_EM.
- alpS = couplingsPtr->alphaS(Q2RenSave);
- alpEM = couplingsPtr->alphaEM(Q2RenSave);
-
-}
-
-//--------------------------------------------------------------------------
-
-// Calculate modified masses and four-vectors for matrix elements.
-
-bool Sigma1Process::setupForME() {
-
- // Common initial-state handling.
- bool allowME = setupForMEin();
-
- // Final state trivial here.
- mME[2] = mH;
- pME[2] = Vec4( 0., 0., 0., mH);
-
- // Done.
- return allowME;
-
-}
-
-//==========================================================================
-
-// The Sigma2Process class.
-// Base class for resolved 2 -> 2 cross sections; derived from SigmaProcess.
-
-//--------------------------------------------------------------------------
-
-// Input and complement kinematics for resolved 2 -> 2 process.
-
-void Sigma2Process::store2Kin( double x1in, double x2in, double sHin,
- double tHin, double m3in, double m4in, double runBW3in, double runBW4in) {
-
- // Default ordering of particles 3 and 4.
- swapTU = false;
-
- // Incoming parton momentum fractions.
- x1Save = x1in;
- x2Save = x2in;
-
- // Incoming masses and their squares.
- bool masslessKin = (id3Mass() == 0) && (id4Mass() == 0);
- if (masslessKin) {
- m3 = 0.;
- m4 = 0.;
- } else {
- m3 = m3in;
- m4 = m4in;
- }
- mSave[3] = m3;
- mSave[4] = m4;
- s3 = m3 * m3;
- s4 = m4 * m4;
-
- // Standard Mandelstam variables and their squares.
- sH = sHin;
- tH = tHin;
- uH = (masslessKin) ? -(sH + tH) : s3 + s4 - (sH + tH);
- mH = sqrt(sH);
- sH2 = sH * sH;
- tH2 = tH * tH;
- uH2 = uH * uH;
-
- // The nominal Breit-Wigner factors with running width.
- runBW3 = runBW3in;
- runBW4 = runBW4in;
-
- // Calculate squared transverse momentum.
- pT2 = (masslessKin) ? tH * uH / sH : (tH * uH - s3 * s4) / sH;
-
- // Special case: pick scale as if 2 -> 1 process in disguise.
- if (isSChannel()) {
-
- // Different options for renormalization scale, but normally sHat.
- Q2RenSave = renormMultFac * sH;
- if (renormScale1 == 2) Q2RenSave = renormFixScale;
-
- // Different options for factorization scale, but normally sHat.
- Q2FacSave = factorMultFac * sH;
- if (factorScale1 == 2) Q2FacSave = factorFixScale;
-
- // Normal case with "true" 2 -> 2.
- } else {
-
- // Different options for renormalization scale.
- if (masslessKin) Q2RenSave = (renormScale2 < 4) ? pT2 : sH;
- else if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
- else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
- else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
- else Q2RenSave = sH;
- Q2RenSave *= renormMultFac;
- if (renormScale2 == 5) Q2RenSave = renormFixScale;
-
- // Different options for factorization scale.
- if (masslessKin) Q2FacSave = (factorScale2 < 4) ? pT2 : sH;
- else if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
- else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
- else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
- else Q2FacSave = sH;
- Q2FacSave *= factorMultFac;
- if (factorScale2 == 5) Q2FacSave = factorFixScale;
- }
-
- // Evaluate alpha_strong and alpha_EM.
- alpS = couplingsPtr->alphaS(Q2RenSave);
- alpEM = couplingsPtr->alphaEM(Q2RenSave);
-
-}
-
-//--------------------------------------------------------------------------
-
-// As above, special kinematics for multiparton interactions.
-
-void Sigma2Process::store2KinMPI( double x1in, double x2in,
- double sHin, double tHin, double uHin, double alpSin, double alpEMin,
- bool needMasses, double m3in, double m4in) {
-
- // Default ordering of particles 3 and 4.
- swapTU = false;
-
- // Incoming x values.
- x1Save = x1in;
- x2Save = x2in;
-
- // Standard Mandelstam variables and their squares.
- sH = sHin;
- tH = tHin;
- uH = uHin;
- mH = sqrt(sH);
- sH2 = sH * sH;
- tH2 = tH * tH;
- uH2 = uH * uH;
-
- // Strong and electroweak couplings.
- alpS = alpSin;
- alpEM = alpEMin;
-
- // Assume vanishing masses. (Will be modified in final kinematics.)
- m3 = 0.;
- s3 = 0.;
- m4 = 0.;
- s4 = 0.;
- sHBeta = sH;
-
- // Scattering angle.
- cosTheta = (tH - uH) / sH;
- sinTheta = 2. * sqrtpos( tH * uH ) / sH;
-
- // In some cases must use masses and redefine meaning of tHat and uHat.
- if (needMasses) {
- m3 = m3in;
- s3 = m3 * m3;
- m4 = m4in;
- s4 = m4 * m4;
- sHMass = sH - s3 - s4;
- sHBeta = sqrtpos(sHMass*sHMass - 4. * s3 * s4);
- tH = -0.5 * (sHMass - sHBeta * cosTheta);
- uH = -0.5 * (sHMass + sHBeta * cosTheta);
- tH2 = tH * tH;
- uH2 = uH * uH;
- }
-
- // pT2 with masses (at this stage) included.
- pT2Mass = 0.25 * sHBeta * pow2(sinTheta);
-
-}
-
-//--------------------------------------------------------------------------
-
-// Perform kinematics for a multiparton interaction, including a rescattering.
-
-bool Sigma2Process::final2KinMPI( int i1Res, int i2Res, Vec4 p1Res, Vec4 p2Res,
- double m1Res, double m2Res) {
-
- // Have to set flavours and colours.
- setIdColAcol();
-
- // Check that masses of outgoing particles not too big.
- m3 = particleDataPtr->m0(idSave[3]);
- m4 = particleDataPtr->m0(idSave[4]);
- mH = sqrt(sH);
- if (m3 + m4 + MASSMARGIN > mH) return false;
- s3 = m3 * m3;
- s4 = m4 * m4;
-
- // Do kinematics of the production; without or with masses.
- double e1In = 0.5 * mH;
- double e2In = e1In;
- double pzIn = e1In;
- if (i1Res > 0 || i2Res > 0) {
- double s1 = m1Res * m1Res;
- double s2 = m2Res * m2Res;
- e1In = 0.5 * (sH + s1 - s2) / mH;
- e2In = 0.5 * (sH + s2 - s1) / mH;
- pzIn = sqrtpos( e1In*e1In - s1 );
- }
-
- // Do kinematics of the decay.
- double e3 = 0.5 * (sH + s3 - s4) / mH;
- double e4 = 0.5 * (sH + s4 - s3) / mH;
- double pAbs = sqrtpos( e3*e3 - s3 );
- phi = 2. * M_PI * rndmPtr->flat();
- double pZ = pAbs * cosTheta;
- pTFin = pAbs * sinTheta;
- double pX = pTFin * sin(phi);
- double pY = pTFin * cos(phi);
- double scale = 0.5 * mH * sinTheta;
-
- // Fill particle info.
- int status1 = (i1Res == 0) ? -31 : -34;
- int status2 = (i2Res == 0) ? -31 : -34;
- parton[1] = Particle( idSave[1], status1, 0, 0, 3, 4,
- colSave[1], acolSave[1], 0., 0., pzIn, e1In, m1Res, scale);
- parton[2] = Particle( idSave[2], status2, 0, 0, 3, 4,
- colSave[2], acolSave[2], 0., 0., -pzIn, e2In, m2Res, scale);
- parton[3] = Particle( idSave[3], 33, 1, 2, 0, 0,
- colSave[3], acolSave[3], pX, pY, pZ, e3, m3, scale);
- parton[4] = Particle( idSave[4], 33, 1, 2, 0, 0,
- colSave[4], acolSave[4], -pX, -pY, -pZ, e4, m4, scale);
-
- // Boost particles from subprocess rest frame to event rest frame.
- // Normal multiparton interaction: only longitudinal boost.
- if (i1Res == 0 && i2Res == 0) {
- double betaZ = (x1Save - x2Save) / (x1Save + x2Save);
- for (int i = 1; i <= 4; ++i) parton[i].bst(0., 0., betaZ);
- // Rescattering: generic rotation and boost required.
- } else {
- RotBstMatrix M;
- M.fromCMframe( p1Res, p2Res);
- for (int i = 1; i <= 4; ++i) parton[i].rotbst(M);
- }
-
- // Done.
- return true;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Calculate modified masses and four-vectors for matrix elements.
-
-bool Sigma2Process::setupForME() {
-
- // Common initial-state handling.
- bool allowME = setupForMEin();
-
- // Correct outgoing c, b, mu and tau to be massive or not.
- mME[2] = m3;
- int id3Tmp = abs(id3Mass());
- if (id3Tmp == 4) mME[2] = mcME;
- if (id3Tmp == 5) mME[2] = mbME;
- if (id3Tmp == 13) mME[2] = mmuME;
- if (id3Tmp == 15) mME[2] = mtauME;
- mME[3] = m4;
- int id4Tmp = abs(id4Mass());
- if (id4Tmp == 4) mME[3] = mcME;
- if (id4Tmp == 5) mME[3] = mbME;
- if (id4Tmp == 13) mME[3] = mmuME;
- if (id4Tmp == 15) mME[3] = mtauME;
-
- // If kinematically impossible turn to massless case, but set error.
- if (mME[2] + mME[3] >= mH) {
- mME[2] = 0.;
- mME[3] = 0.;
- allowME = false;
- }
-
- // Calculate scattering angle in subsystem rest frame.
- double sH34 = sqrtpos( pow2(sH - s3 - s4) - 4. * s3 * s4);
- double cThe = (tH - uH) / sH34;
- double sThe = sqrtpos(1. - cThe * cThe);
-
- // Setup massive kinematics with preserved scattering angle.
- double s3ME = pow2(mME[2]);
- double s4ME = pow2(mME[3]);
- double sH34ME = sqrtpos( pow2(sH - s3ME - s4ME) - 4. * s3ME * s4ME);
- double pAbsME = 0.5 * sH34ME / mH;
-
- // Normally allowed with unequal (or vanishing) masses.
- if (id3Tmp == 0 || id3Tmp != id4Tmp) {
- pME[2] = Vec4( pAbsME * sThe, 0., pAbsME * cThe,
- 0.5 * (sH + s3ME - s4ME) / mH);
- pME[3] = Vec4( -pAbsME * sThe, 0., -pAbsME * cThe,
- 0.5 * (sH + s4ME - s3ME) / mH);
-
- // For equal (anti)particles (e.g. W+ W-) use averaged mass.
- } else {
- mME[2] = sqrtpos(0.5 * (s3ME + s4ME) - 0.25 * pow2(s3ME - s4ME) / sH);
- mME[3] = mME[2];
- pME[2] = Vec4( pAbsME * sThe, 0., pAbsME * cThe, 0.5 * mH);
- pME[3] = Vec4( -pAbsME * sThe, 0., -pAbsME * cThe, 0.5 * mH);
- }
-
- // Done.
- return allowME;
-
-}
-
-//==========================================================================
-
-// The Sigma3Process class.
-// Base class for resolved 2 -> 3 cross sections; derived from SigmaProcess.
-
-//--------------------------------------------------------------------------
-
-// Input and complement kinematics for resolved 2 -> 3 process.
-
-void Sigma3Process::store3Kin( double x1in, double x2in, double sHin,
- Vec4 p3cmIn, Vec4 p4cmIn, Vec4 p5cmIn, double m3in, double m4in,
- double m5in, double runBW3in, double runBW4in, double runBW5in) {
-
- // Default ordering of particles 3 and 4 - not relevant here.
- swapTU = false;
-
- // Incoming parton momentum fractions.
- x1Save = x1in;
- x2Save = x2in;
-
- // Incoming masses and their squares.
- if (id3Mass() == 0 && id4Mass() == 0 && id5Mass() == 0) {
- m3 = 0.;
- m4 = 0.;
- m5 = 0.;
- } else {
- m3 = m3in;
- m4 = m4in;
- m5 = m5in;
- }
- mSave[3] = m3;
- mSave[4] = m4;
- mSave[5] = m5;
- s3 = m3 * m3;
- s4 = m4 * m4;
- s5 = m5 * m5;
-
- // Standard Mandelstam variables and four-momenta in rest frame.
- sH = sHin;
- mH = sqrt(sH);
- sH2 = sH * sH;
- p3cm = p3cmIn;
- p4cm = p4cmIn;
- p5cm = p5cmIn;
-
- // The nominal Breit-Wigner factors with running width.
- runBW3 = runBW3in;
- runBW4 = runBW4in;
- runBW5 = runBW5in;
-
- // Special case: pick scale as if 2 -> 1 process in disguise.
- if (isSChannel()) {
-
- // Different options for renormalization scale, but normally sHat.
- Q2RenSave = renormMultFac * sH;
- if (renormScale1 == 2) Q2RenSave = renormFixScale;
-
- // Different options for factorization scale, but normally sHat.
- Q2FacSave = factorMultFac * sH;
- if (factorScale1 == 2) Q2RenSave = factorFixScale;
-
- // "Normal" 2 -> 3 processes, i.e. not vector boson fusion.
- } else if ( idTchan1() != 23 && idTchan1() != 24 && idTchan2() != 23
- && idTchan2() != 24 ) {
- double mT3S = s3 + p3cm.pT2();
- double mT4S = s4 + p4cm.pT2();
- double mT5S = s5 + p5cm.pT2();
-
- // Different options for renormalization scale.
- if (renormScale3 == 1) Q2RenSave = min( mT3S, min(mT4S, mT5S) );
- else if (renormScale3 == 2) Q2RenSave = sqrt( mT3S * mT4S * mT5S
- / max( mT3S, max(mT4S, mT5S) ) );
- else if (renormScale3 == 3) Q2RenSave = pow( mT3S * mT4S * mT5S,
- 1./3. );
- else if (renormScale3 == 4) Q2RenSave = (mT3S + mT4S + mT5S) / 3.;
- else Q2RenSave = sH;
- Q2RenSave *= renormMultFac;
- if (renormScale3 == 6) Q2RenSave = renormFixScale;
-
- // Different options for factorization scale.
- if (factorScale3 == 1) Q2FacSave = min( mT3S, min(mT4S, mT5S) );
- else if (factorScale3 == 2) Q2FacSave = sqrt( mT3S * mT4S * mT5S
- / max( mT3S, max(mT4S, mT5S) ) );
- else if (factorScale3 == 3) Q2FacSave = pow( mT3S * mT4S * mT5S,
- 1./3. );
- else if (factorScale3 == 4) Q2FacSave = (mT3S + mT4S + mT5S) / 3.;
- else Q2FacSave = sH;
- Q2FacSave *= factorMultFac;
- if (factorScale3 == 6) Q2FacSave = factorFixScale;
-
- // Vector boson fusion 2 -> 3 processes; recoils in positions 4 and 5.
- } else {
- double sV4 = pow2( particleDataPtr->m0(idTchan1()) );
- double sV5 = pow2( particleDataPtr->m0(idTchan2()) );
- double mT3S = s3 + p3cm.pT2();
- double mTV4S = sV4 + p4cm.pT2();
- double mTV5S = sV5 + p5cm.pT2();
-
- // Different options for renormalization scale.
- if (renormScale3VV == 1) Q2RenSave = max( sV4, sV5);
- else if (renormScale3VV == 2) Q2RenSave = sqrt( mTV4S * mTV5S );
- else if (renormScale3VV == 3) Q2RenSave = pow( mT3S * mTV4S * mTV5S,
- 1./3. );
- else if (renormScale3VV == 4) Q2RenSave = (mT3S * mTV4S * mTV5S) / 3.;
- else Q2RenSave = sH;
- Q2RenSave *= renormMultFac;
- if (renormScale3VV == 6) Q2RenSave = renormFixScale;
-
- // Different options for factorization scale.
- if (factorScale3VV == 1) Q2FacSave = max( sV4, sV5);
- else if (factorScale3VV == 2) Q2FacSave = sqrt( mTV4S * mTV5S );
- else if (factorScale3VV == 3) Q2FacSave = pow( mT3S * mTV4S * mTV5S,
- 1./3. );
- else if (factorScale3VV == 4) Q2FacSave = (mT3S * mTV4S * mTV5S) / 3.;
- else Q2FacSave = sH;
- Q2FacSave *= factorMultFac;
- if (factorScale3VV == 6) Q2FacSave = factorFixScale;
- }
-
- // Evaluate alpha_strong and alpha_EM.
- alpS = couplingsPtr->alphaS(Q2RenSave);
- alpEM = couplingsPtr->alphaEM(Q2RenSave);
-
-}
-
-//--------------------------------------------------------------------------
-
-// Calculate modified masses and four-vectors for matrix elements.
-
-bool Sigma3Process::setupForME() {
-
- // Common initial-state handling.
- bool allowME = setupForMEin();
-
- // Correct outgoing c, b, mu and tau to be massive or not.
- mME[2] = m3;
- int id3Tmp = abs(id3Mass());
- if (id3Tmp == 4) mME[2] = mcME;
- if (id3Tmp == 5) mME[2] = mbME;
- if (id3Tmp == 13) mME[2] = mmuME;
- if (id3Tmp == 15) mME[2] = mtauME;
- mME[3] = m4;
- int id4Tmp = abs(id4Mass());
- if (id4Tmp == 4) mME[3] = mcME;
- if (id4Tmp == 5) mME[3] = mbME;
- if (id4Tmp == 13) mME[3] = mmuME;
- if (id4Tmp == 15) mME[3] = mtauME;
- mME[4] = m5;
- int id5Tmp = abs(id5Mass());
- if (id5Tmp == 4) mME[4] = mcME;
- if (id5Tmp == 5) mME[4] = mbME;
- if (id5Tmp == 13) mME[4] = mmuME;
- if (id5Tmp == 15) mME[4] = mtauME;
-
- // If kinematically impossible turn to massless case, but set error.
- if (mME[2] + mME[3] + mME[4] >= mH) {
- mME[2] = 0.;
- mME[3] = 0.;
- mME[4] = 0.;
- allowME = false;
- }
-
- // Form new average masses if identical particles.
- if (id3Tmp != 0 && id4Tmp == id3Tmp && id5Tmp == id3Tmp) {
- double mAvg = (mME[2] + mME[3] + mME[4]) / 3.;
- mME[2] = mAvg;
- mME[3] = mAvg;
- mME[4] = mAvg;
- } else if (id3Tmp != 0 && id4Tmp == id3Tmp) {
- mME[2] = sqrtpos(0.5 * (pow2(mME[2]) + pow2(mME[3]))
- - 0.25 * pow2(pow2(mME[2]) - pow2(mME[3])) / sH);
- mME[3] = mME[2];
- } else if (id3Tmp != 0 && id5Tmp == id3Tmp) {
- mME[2] = sqrtpos(0.5 * (pow2(mME[2]) + pow2(mME[4]))
- - 0.25 * pow2(pow2(mME[2]) - pow2(mME[4])) / sH);
- mME[4] = mME[2];
- } else if (id4Tmp != 0 && id5Tmp == id4Tmp) {
- mME[3] = sqrtpos(0.5 * (pow2(mME[3]) + pow2(mME[4]))
- - 0.25 * pow2(pow2(mME[3]) - pow2(mME[4])) / sH);
- mME[4] = mME[2];
- }
-
- // Iterate rescaled three-momenta until convergence.
- double m2ME3 = pow2(mME[2]);
- double m2ME4 = pow2(mME[3]);
- double m2ME5 = pow2(mME[4]);
- double p2ME3 = p3cm.pAbs2();
- double p2ME4 = p4cm.pAbs2();
- double p2ME5 = p5cm.pAbs2();
- double p2sum = p2ME3 + p2ME4 + p2ME5;
- double eME3 = sqrt(m2ME3 + p2ME3);
- double eME4 = sqrt(m2ME4 + p2ME4);
- double eME5 = sqrt(m2ME5 + p2ME5);
- double esum = eME3 + eME4 + eME5;
- double p2rat = p2ME3 / eME3 + p2ME4 / eME4 + p2ME5 / eME5;
- int iStep = 0;
- while ( abs(esum - mH) > COMPRELERR * mH && iStep < NCOMPSTEP ) {
- ++iStep;
- double compFac = 1. + 2. * (mH - esum) / p2rat;
- p2ME3 *= compFac;
- p2ME4 *= compFac;
- p2ME5 *= compFac;
- eME3 = sqrt(m2ME3 + p2ME3);
- eME4 = sqrt(m2ME4 + p2ME4);
- eME5 = sqrt(m2ME5 + p2ME5);
- esum = eME3 + eME4 + eME5;
- p2rat = p2ME3 / eME3 + p2ME4 / eME4 + p2ME5 / eME5;
- }
-
- // If failed convergence set error flag.
- if (abs(esum - mH) > COMPRELERR * mH) allowME = false;
-
- // Set up accepted kinematics.
- double totFac = sqrt( (p2ME3 + p2ME4 + p2ME5) / p2sum);
- pME[2] = totFac * p3cm;
- pME[2].e( eME3);
- pME[3] = totFac * p4cm;
- pME[3].e( eME4);
- pME[4] = totFac * p5cm;
- pME[4].e( eME5);
-
- // Done.
- return allowME;
-
-}
-
-//==========================================================================
-
-// The SigmaLHAProcess class.
-// Wrapper for Les Houches Accord external input; derived from SigmaProcess.
-// Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
-
-//--------------------------------------------------------------------------
-
-// Evaluate weight for decay angles.
-
-double SigmaLHAProcess::weightDecay( Event& process, int iResBeg,
- int iResEnd) {
-
- // Do nothing if decays present already at input.
- if (iResBeg < process.savedSizeValue()) return 1.;
-
- // Identity of mother of decaying reseonance(s).
- int idMother = process[process[iResBeg].mother1()].idAbs();
-
- // For Higgs decay hand over to standard routine.
- if (idMother == 25 || idMother == 35 || idMother == 36)
- return weightHiggsDecay( process, iResBeg, iResEnd);
-
- // For top decay hand over to standard routine.
- if (idMother == 6)
- return weightTopDecay( process, iResBeg, iResEnd);
-
- // Else done.
- return 1.;
-
-}
-
-//--------------------------------------------------------------------------
-
-// Set scale, alpha_strong and alpha_EM when not set.
-
-void SigmaLHAProcess::setScale() {
-
- // If scale has not been set, then to set.
- double scaleLHA = lhaUpPtr->scale();
- if (scaleLHA < 0.) {
-
- // Final-state partons and their invariant mass.
- vector<int> iFin;
- Vec4 pFinSum;
- for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
- if (lhaUpPtr->mother1(i) == 1) {
- iFin.push_back(i);
- pFinSum += Vec4( lhaUpPtr->px(i), lhaUpPtr->py(i),
- lhaUpPtr->pz(i), lhaUpPtr->e(i) );
- }
- int nFin = iFin.size();
- sH = pFinSum * pFinSum;
- mH = sqrt(sH);
- sH2 = sH * sH;
-
- // If 1 final-state particle then use Sigma1Process logic.
- if (nFin == 1) {
- Q2RenSave = renormMultFac * sH;
- if (renormScale1 == 2) Q2RenSave = renormFixScale;
- Q2FacSave = factorMultFac * sH;
- if (factorScale1 == 2) Q2FacSave = factorFixScale;
-
- // If 2 final-state particles then use Sigma2Process logic.
- } else if (nFin == 2) {
- double s3 = pow2(lhaUpPtr->m(iFin[0]));
- double s4 = pow2(lhaUpPtr->m(iFin[1]));
- double pT2 = pow2(lhaUpPtr->px(iFin[0])) + pow2(lhaUpPtr->py(iFin[0]));
- if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
- else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
- else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
- else Q2RenSave = sH;
- Q2RenSave *= renormMultFac;
- if (renormScale2 == 5) Q2RenSave = renormFixScale;
- if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
- else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
- else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
- else Q2FacSave = sH;
- Q2FacSave *= factorMultFac;
- if (factorScale2 == 5) Q2FacSave = factorFixScale;
-
- // If 3 or more final-state particles then use Sigma3Process logic.
- } else {
- double mTSlow = sH;
- double mTSmed = sH;
- double mTSprod = 1.;
- double mTSsum = 0.;
- for (int i = 0; i < nFin; ++i) {
- double mTSnow = pow2(lhaUpPtr->m(iFin[i]))
- + pow2(lhaUpPtr->px(iFin[i])) + pow2(lhaUpPtr->py(iFin[i]));
- if (mTSnow < mTSlow) {mTSmed = mTSlow; mTSlow = mTSnow;}
- else if (mTSnow < mTSmed) mTSmed = mTSnow;
- mTSprod *= mTSnow;
- mTSsum += mTSnow;
- }
- if (renormScale3 == 1) Q2RenSave = mTSlow;
- else if (renormScale3 == 2) Q2RenSave = sqrt(mTSlow * mTSmed);
- else if (renormScale3 == 3) Q2RenSave = pow(mTSprod, 1. / nFin);
- else if (renormScale3 == 4) Q2RenSave = mTSsum / nFin;
- else Q2RenSave = sH;
- Q2RenSave *= renormMultFac;
- if (renormScale3 == 6) Q2RenSave = renormFixScale;
- if (factorScale3 == 1) Q2FacSave = mTSlow;
- else if (factorScale3 == 2) Q2FacSave = sqrt(mTSlow * mTSmed);
- else if (factorScale3 == 3) Q2FacSave = pow(mTSprod, 1. / nFin);
- else if (factorScale3 == 4) Q2FacSave = mTSsum / nFin;
- else Q2FacSave = sH;
- Q2FacSave *= factorMultFac;
- if (factorScale3 == 6) Q2FacSave = factorFixScale;
- }
- }
-
- // If alpha_strong and alpha_EM have not been set, then set them.
- if (lhaUpPtr->alphaQCD() < 0.001) {
- double Q2RenNow = (scaleLHA < 0.) ? Q2RenSave : pow2(scaleLHA);
- alpS = couplingsPtr->alphaS(Q2RenNow);
- }
- if (lhaUpPtr->alphaQED() < 0.001) {
- double Q2RenNow = (scaleLHA < 0.) ? Q2RenSave : pow2(scaleLHA);
- alpEM = couplingsPtr->alphaEM(Q2RenNow);
- }
-
-}
-
-//--------------------------------------------------------------------------
-
-// Obtain number of final-state partons from LHA object.
-
-int SigmaLHAProcess::nFinal() const {
-
- // At initialization size unknown, so return 0.
- if (lhaUpPtr->sizePart() <= 0) return 0;
-
- // Sum up all particles that has first mother = 1.
- int nFin = 0;
- for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
- if (lhaUpPtr->mother1(i) == 1) ++nFin;
- return nFin;
-
-}
-
-//==========================================================================
-
-} // end namespace Pythia8
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff