From b603a507960b8d0faa987dae28bdac948aaac9f6 Mon Sep 17 00:00:00 2001 From: Saptaparna Date: Mon, 24 Oct 2022 10:11:29 -0500 Subject: [PATCH 1/3] Pythia Update to fix tau decays --- src/HelicityMatrixElements.cc | 140 ++++----------------------- src/TauDecays.cc | 174 +++++++++++++++++++--------------- 2 files changed, 114 insertions(+), 200 deletions(-) diff --git a/src/HelicityMatrixElements.cc b/src/HelicityMatrixElements.cc index 076481091ee..d3dade812e4 100644 --- a/src/HelicityMatrixElements.cc +++ b/src/HelicityMatrixElements.cc @@ -1,5 +1,5 @@ // HelicityMatrixElements.cc is a part of the PYTHIA event generator. -// Copyright (C) 2021 Philip Ilten, Torbjorn Sjostrand. +// Copyright (C) 2019 Philip Ilten, Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. @@ -19,10 +19,10 @@ namespace Pythia8 { // Initialize the helicity matrix element. void HelicityMatrixElement::initPointers(ParticleData* particleDataPtrIn, - CoupSM* coupSMPtrIn, Settings* settingsPtrIn) { + Couplings* couplingsPtrIn, Settings* settingsPtrIn) { particleDataPtr = particleDataPtrIn; - coupSMPtr = coupSMPtrIn; + couplingsPtr = couplingsPtrIn; settingsPtr = settingsPtrIn; for(int i = 0; i <= 5; i++) gamma.push_back(GammaMatrix(i)); @@ -405,18 +405,18 @@ complex HMETwoFermions2W2TwoFermions::calculateME(vector h) { void HMETwoFermions2GammaZ2TwoFermions::initConstants() { // Set the Weinberg angle. - sin2W = coupSMPtr->sin2thetaW(); - cos2W = coupSMPtr->cos2thetaW(); + sin2W = couplingsPtr->sin2thetaW(); + cos2W = couplingsPtr->cos2thetaW(); // Set the on-shell Z/Z' mass and width. zG = particleDataPtr->mWidth(23); zM = particleDataPtr->m0(23); zpG = particleDataPtr->mWidth(32); zpM = particleDataPtr->m0(32); // Set the Z vector and axial couplings to the fermions. - p0CAZ = coupSMPtr->af(abs(pID[0])); - p2CAZ = coupSMPtr->af(abs(pID[2])); - p0CVZ = coupSMPtr->vf(abs(pID[0])); - p2CVZ = coupSMPtr->vf(abs(pID[2])); + p0CAZ = couplingsPtr->af(abs(pID[0])); + p2CAZ = couplingsPtr->af(abs(pID[2])); + p0CVZ = couplingsPtr->vf(abs(pID[0])); + p2CVZ = couplingsPtr->vf(abs(pID[2])); // Turn off the gamma/Z/Z' channels. includeGamma = false; includeZ = false; @@ -455,9 +455,6 @@ void HMETwoFermions2GammaZ2TwoFermions::initConstants() { else if (abs(pID[4]) == 23) includeZ = true; else if (abs(pID[4]) == 32) includeZp = true; } - // Massless approximation for ffbar -> Z -> ffbar. - sMin = settingsPtr->parm("TauDecays:mMinForZ"); - if (sMin > 0) sMin = pow2(sMin); } @@ -480,8 +477,8 @@ void HMETwoFermions2GammaZ2TwoFermions::initWaves(vector& p) // Center of mass energy. s = max( 1., pow2(p[4].m())); // Check if incoming fermions are oriented along z-axis. - zaxis = (p[0].pAbs() == abs(p[0].pz())) && - (p[1].pAbs() == abs(p[1].pz())); + zaxis = (p[0].pAbs() == fabs(p[0].pz())) && + (p[1].pAbs() == fabs(p[1].pz())); } @@ -525,13 +522,9 @@ complex HMETwoFermions2GammaZ2TwoFermions::calculateZME( vector h, double m, double g, double p0CA, double p2CA, double p0CV, double p2CV) { - // Return zero if correct helicity conditions. - if (h[0] == h[1] && zaxis) return complex(0,0); - - // Return massless approximation, otherwise full calculation. - if (sMin >= 0 && s > sMin) return calculateZMEMasslessFermions( - h, m, g, p0CA, p2CA, p0CV, p2CV); complex answer(0,0); + // Return zero if correct helicity conditions. + if (h[0] == h[1] && zaxis) return answer; for (int mu = 0; mu <= 3; mu++) { for (int nu = 0; nu <= 3; nu++) { answer += @@ -549,26 +542,6 @@ complex HMETwoFermions2GammaZ2TwoFermions::calculateZME( //-------------------------------------------------------------------------- -// Return Z/Z' element for helicity matrix element. - -complex HMETwoFermions2GammaZ2TwoFermions::calculateZMEMasslessFermions( - vector h, double m, double g, double p0CA, double p2CA, double p0CV, - double p2CV) { - - complex answer(0,0); - for (int mu = 0; mu <= 3; mu++) { - answer += - (u[1][h[pMap[1]]] * gamma[mu] * (p0CV - p0CA * gamma[5]) * - u[0][h[pMap[0]]]) * gamma[4](mu,mu) * - (u[3][h[pMap[3]]] * gamma[mu] * (p2CV - p2CA * gamma[5]) * - u[2][h[pMap[2]]]); - } - return answer / (16 * pow2(sin2W * cos2W) * (s - m*m + complex(0, s*g/m))); - -} - -//-------------------------------------------------------------------------- - // Return the Z' vector or axial coupling for a fermion. double HMETwoFermions2GammaZ2TwoFermions::zpCoupling(int id, string type) { @@ -597,85 +570,6 @@ double HMETwoFermions2GammaZ2TwoFermions::zpCoupling(int id, string type) { //========================================================================== -// Helicity matrix element for two photons -> two fermions. - -//-------------------------------------------------------------------------- - -// Initialize wave functions for the helicity matrix element. - -void HMETwoGammas2TwoFermions::initWaves(vector& p) { - - u.clear(); - pMap.resize(4); - // Initialize fermion wave functions. - pMap[0] = 0; pMap[1] = 1; pMap[2] = 2; pMap[3] = 3; - vector u0, u1; - for (int h = 0; h < p[0].spinStates(); h++) u0.push_back(p[0].wave(h)); - for (int h = 0; h < p[1].spinStates(); h++) u1.push_back(p[1].wave(h)); - u.push_back(u0); - u.push_back(u1); - setFermionLine(2, p[2], p[3]); - q0 = p[pID[2] > 0 ? 2 : 3].p() - p[0].p(); - q1 = p[pID[2] > 0 ? 2 : 3].p() - p[1].p(); - m = pM[2]; - tm = q0.m2Calc() - pow2(m); - um = q1.m2Calc() - pow2(m); - -} - -//-------------------------------------------------------------------------- - -// Return element for the helicity matrix element. - -complex HMETwoGammas2TwoFermions::calculateME(vector h) { - - Wave4 &u0(u[0][h[0]]), &u1(u[1][h[1]]), &u2(u[2][h[2]]), &u3(u[3][h[3]]); - const complex &u01(u0(1)), &u03(u0(3)), &u11(u1(1)), &u13(u1(3)); - const complex &u20(u2(0)), &u21(u2(1)), &u22(u2(2)), &u23(u2(3)); - const complex &u30(u3(0)), &u31(u3(1)), &u32(u3(2)), &u33(u3(3)); - const double &q00(q0[0]), &q01(q0[1]), &q03(q0[3]); - const double &q10(q1[0]), &q11(q1[1]), &q13(q1[3]); - complex iu02(-u0(2).imag(), u0(2).real()), iu12(-u1(2).imag(), u1(2).real()); - complex iq02(0, q0[2]), iq12(0, q1[2]); - return (u13*(-(m*(u03*u30 + (u01 + iu02)*u31)) + (q00 - q03)*(u03*u32 - + (u01 + iu02)*u33) - (q01 + iq02)*(u01*u32 - iu02*u32 - u03*u33))*u20 - + (u11 + iu12)*(m*(-(u01*u30) + iu02*u30 + u03*u31) - (q01 - iq02) - *(u03*u32 + (u01 + iu02)*u33) + (q00 + q03)*(u01*u32 - - iu02*u32 - u03*u33))*u20 + (u11 - iu12)*(-(m*(u03*u30 + (u01 + iu02) - *u31)) + (q00 - q03)*(u03*u32 + (u01 + iu02)*u33) - (q01 + iq02)*(u01*u32 - - iu02*u32 - u03*u33))*u21 - u13*(m*(-(u01*u30) + iu02*u30 + u03*u31) - - (q01 - iq02)*(u03*u32 + (u01 + iu02)*u33) + (q00 + q03)*(u01*u32 - - iu02*u32 - u03*u33))*u21 - u13*(-((q00 + q03)*(u03*u30 + (u01 + iu02) - *u31)) + (q01 + iq02)*(-(u01*u30) + iu02*u30 + u03*u31) + m*(u03*u32 - + (u01 + iu02)*u33))*u22 + (-u11 - iu12)*(-((q01 - iq02)*(u03*u30 - + (u01 + iu02)*u31)) + (q00 - q03)*(-(u01*u30) + iu02*u30 + u03*u31) - + m*(u01*u32 - iu02*u32 - u03*u33))*u22 + (-u11 + iu12)*(-((q00 + q03) - *(u03*u30 + (u01 + iu02)*u31)) + (q01 + iq02)*(-(u01*u30) + iu02*u30 - + u03*u31) + m*(u03*u32 + (u01 + iu02)*u33))*u23 + u13*(-((q01 - iq02) - *(u03*u30 + (u01 + iu02)*u31)) + (q00 - q03)*(-(u01*u30) + iu02*u30 - + u03*u31) + m*(u01*u32 - iu02*u32 - u03*u33))*u23)/tm + (u03 - *(-(m*(u13*u30 + (u11 + iu12)*u31)) + (q10 - q13)*(u13*u32 + (u11 + iu12) - *u33) - (q11 + iq12)*(u11*u32 - iu12*u32 - u13*u33))*u20 + (u01 + iu02) - *(m*(-(u11*u30) + iu12*u30 + u13*u31) - (q11 - iq12)*(u13*u32 + (u11 - + iu12)*u33) + (q10 + q13)*(u11*u32 - iu12*u32 - u13*u33))*u20 + (u01 - - iu02)*(-(m*(u13*u30 + (u11 + iu12)*u31)) + (q10 - q13)*(u13*u32 + (u11 - + iu12)*u33) - (q11 + iq12)*(u11*u32 - iu12*u32 - u13*u33))*u21 - u03 - *(m*(-(u11*u30) + iu12*u30 + u13*u31) - (q11 - iq12)*(u13*u32 + (u11 - + iu12)*u33) + (q10 + q13)*(u11*u32 - iu12*u32 - u13*u33))*u21 - u03 - *(-((q10 + q13)*(u13*u30 + (u11 + iu12)*u31)) + (q11 + iq12)*(-(u11*u30) - + iu12*u30 + u13*u31) + m*(u13*u32 + (u11 + iu12)*u33))*u22 + (-u01 - - iu02)*(-((q11 - iq12)*(u13*u30 + (u11 + iu12)*u31)) + (q10 - q13) - *(-(u11*u30) + iu12*u30 + u13*u31) + m*(u11*u32 - iu12*u32 - u13*u33)) - *u22 + (-u01 + iu02)*(-((q10 + q13)*(u13*u30 + (u11 + iu12)*u31)) + (q11 - + iq12)*(-(u11*u30) + iu12*u30 + u13*u31) + m*(u13*u32 + (u11 + iu12) - *u33))*u23 + u03*(-((q11 - iq12)*(u13*u30 + (u11 + iu12)*u31)) + (q10 - - q13)*(-(u11*u30) + iu12*u30 + u13*u31) + m*(u11*u32 - iu12*u32 - u13 - *u33))*u23)/um; - -} - -//========================================================================== - // Helicity matrix element for X -> two fermions. // Base class for the W, photon, and Z -> two fermions helicity matrix @@ -689,7 +583,7 @@ void HMEX2TwoFermions::initWaves(vector& p) { u.clear(); pMap.resize(4); - // Initialize boson wave function. + // Initialize W wave function. vector< Wave4 > u1; pMap[1] = 1; for (int h = 0; h < p[pMap[1]].spinStates(); h++) @@ -784,8 +678,8 @@ complex HMEGamma2TwoFermions::calculateME(vector h) { void HMEZ2TwoFermions::initConstants() { // Set the vector and axial couplings to the fermions. - p2CA = coupSMPtr->af(abs(pID[2])); - p2CV = coupSMPtr->vf(abs(pID[2])); + p2CA = couplingsPtr->af(abs(pID[2])); + p2CV = couplingsPtr->vf(abs(pID[2])); if (settingsPtr && abs(pID[0]) == 32) { p2CA = zpCoupling(abs(pID[2]), "a"); p2CV = zpCoupling(abs(pID[2]), "v"); @@ -983,7 +877,7 @@ double HMETauDecay::decayWeightMax(vector& p) { double on = real(p[0].rho[0][0]) > real(p[0].rho[1][1]) ? real(p[0].rho[0][0]) : real(p[0].rho[1][1]); // Determine the maximum off-diagonal element of rho. - double off = abs(real(p[0].rho[0][1])) + abs(imag(p[0].rho[0][1])); + double off = fabs(real(p[0].rho[0][1])) + fabs(imag(p[0].rho[0][1])); return DECAYWEIGHTMAX * (on + off); } diff --git a/src/TauDecays.cc b/src/TauDecays.cc index 53f27a52efb..dc971902e6c 100644 --- a/src/TauDecays.cc +++ b/src/TauDecays.cc @@ -1,5 +1,5 @@ // TauDecays.cc is a part of the PYTHIA event generator. -// Copyright (C) 2021 Philip Ilten, Torbjorn Sjostrand. +// Copyright (C) 2019 Philip Ilten, Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. @@ -36,53 +36,60 @@ const double TauDecays::WTCORRECTION[11] = { 1., 1., 1., // Additionally, the necessary matrix elements are initialized with the // Standard Model couplings, and particle data pointers. -void TauDecays::init() { +void TauDecays::init(Info* infoPtrIn, Settings* settingsPtrIn, + ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, + Couplings* couplingsPtrIn) { + + // Set the pointers. + infoPtr = infoPtrIn; + settingsPtr = settingsPtrIn; + particleDataPtr = particleDataPtrIn; + rndmPtr = rndmPtrIn; + couplingsPtr = couplingsPtrIn; // Initialize the hard matrix elements. hmeTwoFermions2W2TwoFermions - .initPointers(particleDataPtr, coupSMPtr, settingsPtr); + .initPointers(particleDataPtr, couplingsPtr, settingsPtr); hmeTwoFermions2GammaZ2TwoFermions - .initPointers(particleDataPtr, coupSMPtr, settingsPtr); - hmeTwoGammas2TwoFermions - .initPointers(particleDataPtr, coupSMPtr, settingsPtr); + .initPointers(particleDataPtr, couplingsPtr, settingsPtr); hmeW2TwoFermions - .initPointers(particleDataPtr, coupSMPtr, settingsPtr); + .initPointers(particleDataPtr, couplingsPtr, settingsPtr); hmeZ2TwoFermions - .initPointers(particleDataPtr, coupSMPtr, settingsPtr); + .initPointers(particleDataPtr, couplingsPtr, settingsPtr); hmeGamma2TwoFermions - .initPointers(particleDataPtr, coupSMPtr); + .initPointers(particleDataPtr, couplingsPtr); hmeHiggs2TwoFermions - .initPointers(particleDataPtr, coupSMPtr, settingsPtr); + .initPointers(particleDataPtr, couplingsPtr, settingsPtr); // Initialize the tau decay matrix elements. - hmeTau2Meson .initPointers(particleDataPtr, coupSMPtr); - hmeTau2TwoLeptons .initPointers(particleDataPtr, coupSMPtr); - hmeTau2TwoMesonsViaVector .initPointers(particleDataPtr, coupSMPtr); - hmeTau2TwoMesonsViaVectorScalar.initPointers(particleDataPtr, coupSMPtr); - hmeTau2ThreePions .initPointers(particleDataPtr, coupSMPtr); - hmeTau2ThreeMesonsWithKaons .initPointers(particleDataPtr, coupSMPtr); - hmeTau2ThreeMesonsGeneric .initPointers(particleDataPtr, coupSMPtr); - hmeTau2TwoPionsGamma .initPointers(particleDataPtr, coupSMPtr); - hmeTau2FourPions .initPointers(particleDataPtr, coupSMPtr); - hmeTau2FivePions .initPointers(particleDataPtr, coupSMPtr); - hmeTau2PhaseSpace .initPointers(particleDataPtr, coupSMPtr); + hmeTau2Meson .initPointers(particleDataPtr, couplingsPtr); + hmeTau2TwoLeptons .initPointers(particleDataPtr, couplingsPtr); + hmeTau2TwoMesonsViaVector .initPointers(particleDataPtr, couplingsPtr); + hmeTau2TwoMesonsViaVectorScalar.initPointers(particleDataPtr, couplingsPtr); + hmeTau2ThreePions .initPointers(particleDataPtr, couplingsPtr); + hmeTau2ThreeMesonsWithKaons .initPointers(particleDataPtr, couplingsPtr); + hmeTau2ThreeMesonsGeneric .initPointers(particleDataPtr, couplingsPtr); + hmeTau2TwoPionsGamma .initPointers(particleDataPtr, couplingsPtr); + hmeTau2FourPions .initPointers(particleDataPtr, couplingsPtr); + hmeTau2FivePions .initPointers(particleDataPtr, couplingsPtr); + hmeTau2PhaseSpace .initPointers(particleDataPtr, couplingsPtr); // User selected tau settings. - tauExt = mode("TauDecays:externalMode"); - tauMode = mode("TauDecays:mode"); - tauMother = mode("TauDecays:tauMother"); - tauPol = parm("TauDecays:tauPolarization"); + tauExt = settingsPtr->mode("TauDecays:externalMode"); + tauMode = settingsPtr->mode("TauDecays:mode"); + tauMother = settingsPtr->mode("TauDecays:tauMother"); + tauPol = settingsPtr->parm("TauDecays:tauPolarization"); // Parameters to determine if correlated partner should decay. - limitTau0 = flag("ParticleDecays:limitTau0"); - tau0Max = parm("ParticleDecays:tau0Max"); - limitTau = flag("ParticleDecays:limitTau"); - tauMax = parm("ParticleDecays:tauMax"); - limitRadius = flag("ParticleDecays:limitRadius"); - rMax = parm("ParticleDecays:rMax"); - limitCylinder = flag("ParticleDecays:limitCylinder"); - xyMax = parm("ParticleDecays:xyMax"); - zMax = parm("ParticleDecays:zMax"); + limitTau0 = settingsPtr->flag("ParticleDecays:limitTau0"); + tau0Max = settingsPtr->parm("ParticleDecays:tau0Max"); + limitTau = settingsPtr->flag("ParticleDecays:limitTau"); + tauMax = settingsPtr->parm("ParticleDecays:tauMax"); + limitRadius = settingsPtr->flag("ParticleDecays:limitRadius"); + rMax = settingsPtr->parm("ParticleDecays:rMax"); + limitCylinder = settingsPtr->flag("ParticleDecays:limitCylinder"); + xyMax = settingsPtr->parm("ParticleDecays:xyMax"); + zMax = settingsPtr->parm("ParticleDecays:zMax"); limitDecay = limitTau0 || limitTau || limitRadius || limitCylinder; } @@ -123,15 +130,10 @@ bool TauDecays::decay(int idxOut1, Event& event) { int idxOut2 = event[idxOut2Top].iBotCopyId(); out2 = HelicityParticle(event[idxOut2]); - // Set the mediator of the hard process (also handle no mediator). + // Set the mediator of the hard process. int idxMediator = event[idxOut1Top].mother1(); - if (event[idxOut1Top].mother2() > event[idxOut1Top].mother1() && - event[idxOut1Top].mother1() == event[idxOut2Top].mother1() && - event[idxOut1Top].mother2() == event[idxOut2Top].mother2()) - idxMediator = idxOut1Top; mediator = HelicityParticle(event[idxMediator]); mediator.direction = -1; - if (idxMediator == idxOut1Top) mediator.id(23); if (mediator.m() < out1.m() + out2.m()) { Vec4 p = out1.p() + out2.p(); mediator.p(p); @@ -170,12 +172,6 @@ bool TauDecays::decay(int idxOut1, Event& event) { // Check partner can decay. else if (!out2.canDecay()) correlated = false; else if (!out2.mayDecay()) correlated = false; - // Check partner not EW showered, set decay matrix otherwise. - else if (!out2.isFinal() && out2.statusAbs() > 40 - && out2.statusAbs() < 60) { - correlated = false; - out2.D = out2.rho; - } } // Set the production mechanism. @@ -196,7 +192,7 @@ bool TauDecays::decay(int idxOut1, Event& event) { } } - // Catch unknown production mechanisms. + // Catch unknown production mechanims. if (!known) { particles[1] = mediator; if (abs(mediator.id()) == 22) @@ -262,7 +258,7 @@ bool TauDecays::decay(int idxOut1, Event& event) { // Switch the taus. tau = &particles[idx]; // Calculate second tau's density matrix. - if (hardME) hardME->calculateRho(idx, particles); + hardME->calculateRho(idx, particles); // Decay the second tau. children.clear(); @@ -300,30 +296,29 @@ bool TauDecays::decay(int idxOut1, Event& event) { bool TauDecays::internalMechanism(Event&) { + // Flag if process is known. + bool known = true; + // Produced from a photon, Z, or Z'. if (abs(mediator.id()) == 22 || abs(mediator.id()) == 23 || abs(mediator.id()) == 32) { - // Produced from photons: t-channel. - if (in1.id() == 22 && in2.id() == 22) { - hardME = hmeTwoGammas2TwoFermions.initChannel(particles); // Produced from fermions: s-channel. - } else if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18 && - in1.daughter1() == in2.daughter1() && - in1.daughter2() == in2.daughter2()) { + if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18 && in1.daughter2() == 0 && + in2.daughter2() == 0 && in1.daughter1() == in2.daughter1()) { particles.push_back(mediator); hardME = hmeTwoFermions2GammaZ2TwoFermions.initChannel(particles); // Unknown photon production. - } else return false; + } else known = false; // Produced from a W or W'. } else if (abs(mediator.id()) == 24 || abs(mediator.id()) == 34) { // Produced from fermions: s-channel. if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18 && in1.daughter2() == 0 && - in2.daughter2() == 0 && in1.daughter1() == in2.daughter1()) { + in2.daughter2() == 0 && in1.daughter1() == in2.daughter1()) { particles.push_back(mediator); hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); - // Unknown W production. - } else return false; + // Unknown W production. + } else known = false; // Produced from a Higgs. } else if (abs(mediator.id()) == 25 || abs(mediator.id()) == 35 || @@ -362,29 +357,41 @@ bool TauDecays::internalMechanism(Event&) { hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); // Unknown production. - } else return false; - return true; + } else known = false; + return known; } //-------------------------------------------------------------------------- // Determine the tau polarization and tau decay correlation using the provided -// SPINUP digits interpreted as helicity states. +// SPINUP digits. bool TauDecays::externalMechanism(Event &event) { - // Uncorrelated, take directly from tau SPINUP if valid. + // Flag if process is known. + bool known = true; + + // Uncorrelated, take directly from SPINUP if valid. if (tauExt == 0) correlated = false; if (!correlated) { - if (particles[2].pol() == 9) - particles[2].pol(event[particles[2].iTopCopyId()].pol()); - if (particles[2].pol() == 9) return false; + double spinup = particles[2].pol(); + if (abs(spinup) > 1.001) spinup = event[particles[2].iTopCopyId()].pol(); + if (abs(spinup) > 1.001) known = false; + else { + particles[2].rho[0][0] = (1 - spinup) / 2; + particles[2].rho[1][1] = (1 + spinup) / 2; + } - // Correlated, take from mother SPINUP if valid. + // Correlated, try mother. } else if (tauExt == 1) { - if (mediator.pol() == 9) mediator.pol(event[mediator.iTopCopyId()].pol()); - if (mediator.pol() == 9) return false; + double spinup = mediator.pol(); + if (abs(spinup) > 1.001) spinup = event[mediator.iTopCopyId()].pol(); + if (abs(spinup) > 1.001) spinup = 0; + if (mediator.rho.size() > 1) { + mediator.rho[0][0] = (1 - spinup) / mediator.spinStates(); + mediator.rho[1][1] = (1 + spinup) / mediator.spinStates(); + } particles[1] = mediator; if (abs(mediator.id()) == 22) hardME = hmeGamma2TwoFermions.initChannel(particles); @@ -395,11 +402,11 @@ bool TauDecays::externalMechanism(Event &event) { else if (abs(mediator.id()) == 25 || abs(mediator.id()) == 35 || abs(mediator.id()) == 36 || abs(mediator.id()) == 37) hardME = hmeHiggs2TwoFermions.initChannel(particles); - else return false; + else known = false; - // Unknown mechanism. - } else return false; - return true; + // Unknown mechanism. + } else known = false; + return known; } @@ -608,11 +615,24 @@ void TauDecays::isotropicDecay(vector& children) { // Perform two-particle decays in the respective rest frame. vector pInv(decayMult + 1); for (int i = 1; i < decayMult; ++i) { - // Fill four-momenta - pair ps = - rndmPtr->phaseSpace2(mInv[i], mInv[i+1], children[i].m()); - pInv[i+1].p(ps.first); - children[i].p(ps.second); + double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - children[i].m()) + * (mInv[i] + mInv[i+1] + children[i].m()) + * (mInv[i] + mInv[i+1] - children[i].m()) + * (mInv[i] - mInv[i+1] + children[i].m()) ) / mInv[i]; + + // Isotropic angles give three-momentum. + double cosTheta = 2. * rndmPtr->flat() - 1.; + double sinTheta = sqrt(1. - cosTheta*cosTheta); + double phi = 2. * M_PI * rndmPtr->flat(); + double pX = pAbs * sinTheta * cos(phi); + double pY = pAbs * sinTheta * sin(phi); + double pZ = pAbs * cosTheta; + + // Calculate energies, fill four-momenta. + double eHad = sqrt( children[i].m()*children[i].m() + pAbs*pAbs); + double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs); + children[i].p( pX, pY, pZ, eHad); + pInv[i+1].p( -pX, -pY, -pZ, eInv); } // Boost decay products to the mother rest frame. From b4e3a52ec57256f0bec41950075ce1a5a96e0b5d Mon Sep 17 00:00:00 2001 From: Saptaparna Date: Mon, 24 Oct 2022 14:03:51 -0500 Subject: [PATCH 2/3] Update pythia bug fix -- correct timestamp --- src/HelicityMatrixElements.cc | 140 +++++++++++++++++++++++---- src/TauDecays.cc | 174 +++++++++++++++------------------- 2 files changed, 200 insertions(+), 114 deletions(-) diff --git a/src/HelicityMatrixElements.cc b/src/HelicityMatrixElements.cc index d3dade812e4..7bfadf66b44 100644 --- a/src/HelicityMatrixElements.cc +++ b/src/HelicityMatrixElements.cc @@ -1,5 +1,5 @@ // HelicityMatrixElements.cc is a part of the PYTHIA event generator. -// Copyright (C) 2019 Philip Ilten, Torbjorn Sjostrand. +// Copyright (C) 2022 Philip Ilten, Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. @@ -19,10 +19,10 @@ namespace Pythia8 { // Initialize the helicity matrix element. void HelicityMatrixElement::initPointers(ParticleData* particleDataPtrIn, - Couplings* couplingsPtrIn, Settings* settingsPtrIn) { + CoupSM* coupSMPtrIn, Settings* settingsPtrIn) { particleDataPtr = particleDataPtrIn; - couplingsPtr = couplingsPtrIn; + coupSMPtr = coupSMPtrIn; settingsPtr = settingsPtrIn; for(int i = 0; i <= 5; i++) gamma.push_back(GammaMatrix(i)); @@ -405,18 +405,18 @@ complex HMETwoFermions2W2TwoFermions::calculateME(vector h) { void HMETwoFermions2GammaZ2TwoFermions::initConstants() { // Set the Weinberg angle. - sin2W = couplingsPtr->sin2thetaW(); - cos2W = couplingsPtr->cos2thetaW(); + sin2W = coupSMPtr->sin2thetaW(); + cos2W = coupSMPtr->cos2thetaW(); // Set the on-shell Z/Z' mass and width. zG = particleDataPtr->mWidth(23); zM = particleDataPtr->m0(23); zpG = particleDataPtr->mWidth(32); zpM = particleDataPtr->m0(32); // Set the Z vector and axial couplings to the fermions. - p0CAZ = couplingsPtr->af(abs(pID[0])); - p2CAZ = couplingsPtr->af(abs(pID[2])); - p0CVZ = couplingsPtr->vf(abs(pID[0])); - p2CVZ = couplingsPtr->vf(abs(pID[2])); + p0CAZ = coupSMPtr->af(abs(pID[0])); + p2CAZ = coupSMPtr->af(abs(pID[2])); + p0CVZ = coupSMPtr->vf(abs(pID[0])); + p2CVZ = coupSMPtr->vf(abs(pID[2])); // Turn off the gamma/Z/Z' channels. includeGamma = false; includeZ = false; @@ -455,6 +455,9 @@ void HMETwoFermions2GammaZ2TwoFermions::initConstants() { else if (abs(pID[4]) == 23) includeZ = true; else if (abs(pID[4]) == 32) includeZp = true; } + // Massless approximation for ffbar -> Z -> ffbar. + sMin = settingsPtr->parm("TauDecays:mMinForZ"); + if (sMin > 0) sMin = pow2(sMin); } @@ -477,8 +480,8 @@ void HMETwoFermions2GammaZ2TwoFermions::initWaves(vector& p) // Center of mass energy. s = max( 1., pow2(p[4].m())); // Check if incoming fermions are oriented along z-axis. - zaxis = (p[0].pAbs() == fabs(p[0].pz())) && - (p[1].pAbs() == fabs(p[1].pz())); + zaxis = (p[0].pAbs() == abs(p[0].pz())) && + (p[1].pAbs() == abs(p[1].pz())); } @@ -522,9 +525,13 @@ complex HMETwoFermions2GammaZ2TwoFermions::calculateZME( vector h, double m, double g, double p0CA, double p2CA, double p0CV, double p2CV) { - complex answer(0,0); // Return zero if correct helicity conditions. - if (h[0] == h[1] && zaxis) return answer; + if (h[0] == h[1] && zaxis) return complex(0,0); + + // Return massless approximation, otherwise full calculation. + if (sMin >= 0 && s > sMin) return calculateZMEMasslessFermions( + h, m, g, p0CA, p2CA, p0CV, p2CV); + complex answer(0,0); for (int mu = 0; mu <= 3; mu++) { for (int nu = 0; nu <= 3; nu++) { answer += @@ -542,6 +549,26 @@ complex HMETwoFermions2GammaZ2TwoFermions::calculateZME( //-------------------------------------------------------------------------- +// Return Z/Z' element for helicity matrix element. + +complex HMETwoFermions2GammaZ2TwoFermions::calculateZMEMasslessFermions( + vector h, double m, double g, double p0CA, double p2CA, double p0CV, + double p2CV) { + + complex answer(0,0); + for (int mu = 0; mu <= 3; mu++) { + answer += + (u[1][h[pMap[1]]] * gamma[mu] * (p0CV - p0CA * gamma[5]) * + u[0][h[pMap[0]]]) * gamma[4](mu,mu) * + (u[3][h[pMap[3]]] * gamma[mu] * (p2CV - p2CA * gamma[5]) * + u[2][h[pMap[2]]]); + } + return answer / (16 * pow2(sin2W * cos2W) * (s - m*m + complex(0, s*g/m))); + +} + +//-------------------------------------------------------------------------- + // Return the Z' vector or axial coupling for a fermion. double HMETwoFermions2GammaZ2TwoFermions::zpCoupling(int id, string type) { @@ -570,6 +597,85 @@ double HMETwoFermions2GammaZ2TwoFermions::zpCoupling(int id, string type) { //========================================================================== +// Helicity matrix element for two photons -> two fermions. + +//-------------------------------------------------------------------------- + +// Initialize wave functions for the helicity matrix element. + +void HMETwoGammas2TwoFermions::initWaves(vector& p) { + + u.clear(); + pMap.resize(4); + // Initialize fermion wave functions. + pMap[0] = 0; pMap[1] = 1; pMap[2] = 2; pMap[3] = 3; + vector u0, u1; + for (int h = 0; h < p[0].spinStates(); h++) u0.push_back(p[0].wave(h)); + for (int h = 0; h < p[1].spinStates(); h++) u1.push_back(p[1].wave(h)); + u.push_back(u0); + u.push_back(u1); + setFermionLine(2, p[2], p[3]); + q0 = p[pID[2] > 0 ? 2 : 3].p() - p[0].p(); + q1 = p[pID[2] > 0 ? 2 : 3].p() - p[1].p(); + m = pM[2]; + tm = q0.m2Calc() - pow2(m); + um = q1.m2Calc() - pow2(m); + +} + +//-------------------------------------------------------------------------- + +// Return element for the helicity matrix element. + +complex HMETwoGammas2TwoFermions::calculateME(vector h) { + + Wave4 &u0(u[0][h[0]]), &u1(u[1][h[1]]), &u2(u[2][h[2]]), &u3(u[3][h[3]]); + const complex &u01(u0(1)), &u03(u0(3)), &u11(u1(1)), &u13(u1(3)); + const complex &u20(u2(0)), &u21(u2(1)), &u22(u2(2)), &u23(u2(3)); + const complex &u30(u3(0)), &u31(u3(1)), &u32(u3(2)), &u33(u3(3)); + const double &q00(q0[0]), &q01(q0[1]), &q03(q0[3]); + const double &q10(q1[0]), &q11(q1[1]), &q13(q1[3]); + complex iu02(-u0(2).imag(), u0(2).real()), iu12(-u1(2).imag(), u1(2).real()); + complex iq02(0, q0[2]), iq12(0, q1[2]); + return (u13*(-(m*(u03*u30 + (u01 + iu02)*u31)) + (q00 - q03)*(u03*u32 + + (u01 + iu02)*u33) - (q01 + iq02)*(u01*u32 - iu02*u32 - u03*u33))*u20 + + (u11 + iu12)*(m*(-(u01*u30) + iu02*u30 + u03*u31) - (q01 - iq02) + *(u03*u32 + (u01 + iu02)*u33) + (q00 + q03)*(u01*u32 + - iu02*u32 - u03*u33))*u20 + (u11 - iu12)*(-(m*(u03*u30 + (u01 + iu02) + *u31)) + (q00 - q03)*(u03*u32 + (u01 + iu02)*u33) - (q01 + iq02)*(u01*u32 + - iu02*u32 - u03*u33))*u21 - u13*(m*(-(u01*u30) + iu02*u30 + u03*u31) + - (q01 - iq02)*(u03*u32 + (u01 + iu02)*u33) + (q00 + q03)*(u01*u32 + - iu02*u32 - u03*u33))*u21 - u13*(-((q00 + q03)*(u03*u30 + (u01 + iu02) + *u31)) + (q01 + iq02)*(-(u01*u30) + iu02*u30 + u03*u31) + m*(u03*u32 + + (u01 + iu02)*u33))*u22 + (-u11 - iu12)*(-((q01 - iq02)*(u03*u30 + + (u01 + iu02)*u31)) + (q00 - q03)*(-(u01*u30) + iu02*u30 + u03*u31) + + m*(u01*u32 - iu02*u32 - u03*u33))*u22 + (-u11 + iu12)*(-((q00 + q03) + *(u03*u30 + (u01 + iu02)*u31)) + (q01 + iq02)*(-(u01*u30) + iu02*u30 + + u03*u31) + m*(u03*u32 + (u01 + iu02)*u33))*u23 + u13*(-((q01 - iq02) + *(u03*u30 + (u01 + iu02)*u31)) + (q00 - q03)*(-(u01*u30) + iu02*u30 + + u03*u31) + m*(u01*u32 - iu02*u32 - u03*u33))*u23)/tm + (u03 + *(-(m*(u13*u30 + (u11 + iu12)*u31)) + (q10 - q13)*(u13*u32 + (u11 + iu12) + *u33) - (q11 + iq12)*(u11*u32 - iu12*u32 - u13*u33))*u20 + (u01 + iu02) + *(m*(-(u11*u30) + iu12*u30 + u13*u31) - (q11 - iq12)*(u13*u32 + (u11 + + iu12)*u33) + (q10 + q13)*(u11*u32 - iu12*u32 - u13*u33))*u20 + (u01 + - iu02)*(-(m*(u13*u30 + (u11 + iu12)*u31)) + (q10 - q13)*(u13*u32 + (u11 + + iu12)*u33) - (q11 + iq12)*(u11*u32 - iu12*u32 - u13*u33))*u21 - u03 + *(m*(-(u11*u30) + iu12*u30 + u13*u31) - (q11 - iq12)*(u13*u32 + (u11 + + iu12)*u33) + (q10 + q13)*(u11*u32 - iu12*u32 - u13*u33))*u21 - u03 + *(-((q10 + q13)*(u13*u30 + (u11 + iu12)*u31)) + (q11 + iq12)*(-(u11*u30) + + iu12*u30 + u13*u31) + m*(u13*u32 + (u11 + iu12)*u33))*u22 + (-u01 + - iu02)*(-((q11 - iq12)*(u13*u30 + (u11 + iu12)*u31)) + (q10 - q13) + *(-(u11*u30) + iu12*u30 + u13*u31) + m*(u11*u32 - iu12*u32 - u13*u33)) + *u22 + (-u01 + iu02)*(-((q10 + q13)*(u13*u30 + (u11 + iu12)*u31)) + (q11 + + iq12)*(-(u11*u30) + iu12*u30 + u13*u31) + m*(u13*u32 + (u11 + iu12) + *u33))*u23 + u03*(-((q11 - iq12)*(u13*u30 + (u11 + iu12)*u31)) + (q10 + - q13)*(-(u11*u30) + iu12*u30 + u13*u31) + m*(u11*u32 - iu12*u32 - u13 + *u33))*u23)/um; + +} + +//========================================================================== + // Helicity matrix element for X -> two fermions. // Base class for the W, photon, and Z -> two fermions helicity matrix @@ -583,7 +689,7 @@ void HMEX2TwoFermions::initWaves(vector& p) { u.clear(); pMap.resize(4); - // Initialize W wave function. + // Initialize boson wave function. vector< Wave4 > u1; pMap[1] = 1; for (int h = 0; h < p[pMap[1]].spinStates(); h++) @@ -678,8 +784,8 @@ complex HMEGamma2TwoFermions::calculateME(vector h) { void HMEZ2TwoFermions::initConstants() { // Set the vector and axial couplings to the fermions. - p2CA = couplingsPtr->af(abs(pID[2])); - p2CV = couplingsPtr->vf(abs(pID[2])); + p2CA = coupSMPtr->af(abs(pID[2])); + p2CV = coupSMPtr->vf(abs(pID[2])); if (settingsPtr && abs(pID[0]) == 32) { p2CA = zpCoupling(abs(pID[2]), "a"); p2CV = zpCoupling(abs(pID[2]), "v"); @@ -877,7 +983,7 @@ double HMETauDecay::decayWeightMax(vector& p) { double on = real(p[0].rho[0][0]) > real(p[0].rho[1][1]) ? real(p[0].rho[0][0]) : real(p[0].rho[1][1]); // Determine the maximum off-diagonal element of rho. - double off = fabs(real(p[0].rho[0][1])) + fabs(imag(p[0].rho[0][1])); + double off = abs(real(p[0].rho[0][1])) + abs(imag(p[0].rho[0][1])); return DECAYWEIGHTMAX * (on + off); } diff --git a/src/TauDecays.cc b/src/TauDecays.cc index dc971902e6c..cdeaf6f9adb 100644 --- a/src/TauDecays.cc +++ b/src/TauDecays.cc @@ -1,5 +1,5 @@ // TauDecays.cc is a part of the PYTHIA event generator. -// Copyright (C) 2019 Philip Ilten, Torbjorn Sjostrand. +// Copyright (C) 2022 Philip Ilten, Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. @@ -36,60 +36,53 @@ const double TauDecays::WTCORRECTION[11] = { 1., 1., 1., // Additionally, the necessary matrix elements are initialized with the // Standard Model couplings, and particle data pointers. -void TauDecays::init(Info* infoPtrIn, Settings* settingsPtrIn, - ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, - Couplings* couplingsPtrIn) { - - // Set the pointers. - infoPtr = infoPtrIn; - settingsPtr = settingsPtrIn; - particleDataPtr = particleDataPtrIn; - rndmPtr = rndmPtrIn; - couplingsPtr = couplingsPtrIn; +void TauDecays::init() { // Initialize the hard matrix elements. hmeTwoFermions2W2TwoFermions - .initPointers(particleDataPtr, couplingsPtr, settingsPtr); + .initPointers(particleDataPtr, coupSMPtr, settingsPtr); hmeTwoFermions2GammaZ2TwoFermions - .initPointers(particleDataPtr, couplingsPtr, settingsPtr); + .initPointers(particleDataPtr, coupSMPtr, settingsPtr); + hmeTwoGammas2TwoFermions + .initPointers(particleDataPtr, coupSMPtr, settingsPtr); hmeW2TwoFermions - .initPointers(particleDataPtr, couplingsPtr, settingsPtr); + .initPointers(particleDataPtr, coupSMPtr, settingsPtr); hmeZ2TwoFermions - .initPointers(particleDataPtr, couplingsPtr, settingsPtr); + .initPointers(particleDataPtr, coupSMPtr, settingsPtr); hmeGamma2TwoFermions - .initPointers(particleDataPtr, couplingsPtr); + .initPointers(particleDataPtr, coupSMPtr); hmeHiggs2TwoFermions - .initPointers(particleDataPtr, couplingsPtr, settingsPtr); + .initPointers(particleDataPtr, coupSMPtr, settingsPtr); // Initialize the tau decay matrix elements. - hmeTau2Meson .initPointers(particleDataPtr, couplingsPtr); - hmeTau2TwoLeptons .initPointers(particleDataPtr, couplingsPtr); - hmeTau2TwoMesonsViaVector .initPointers(particleDataPtr, couplingsPtr); - hmeTau2TwoMesonsViaVectorScalar.initPointers(particleDataPtr, couplingsPtr); - hmeTau2ThreePions .initPointers(particleDataPtr, couplingsPtr); - hmeTau2ThreeMesonsWithKaons .initPointers(particleDataPtr, couplingsPtr); - hmeTau2ThreeMesonsGeneric .initPointers(particleDataPtr, couplingsPtr); - hmeTau2TwoPionsGamma .initPointers(particleDataPtr, couplingsPtr); - hmeTau2FourPions .initPointers(particleDataPtr, couplingsPtr); - hmeTau2FivePions .initPointers(particleDataPtr, couplingsPtr); - hmeTau2PhaseSpace .initPointers(particleDataPtr, couplingsPtr); + hmeTau2Meson .initPointers(particleDataPtr, coupSMPtr); + hmeTau2TwoLeptons .initPointers(particleDataPtr, coupSMPtr); + hmeTau2TwoMesonsViaVector .initPointers(particleDataPtr, coupSMPtr); + hmeTau2TwoMesonsViaVectorScalar.initPointers(particleDataPtr, coupSMPtr); + hmeTau2ThreePions .initPointers(particleDataPtr, coupSMPtr); + hmeTau2ThreeMesonsWithKaons .initPointers(particleDataPtr, coupSMPtr); + hmeTau2ThreeMesonsGeneric .initPointers(particleDataPtr, coupSMPtr); + hmeTau2TwoPionsGamma .initPointers(particleDataPtr, coupSMPtr); + hmeTau2FourPions .initPointers(particleDataPtr, coupSMPtr); + hmeTau2FivePions .initPointers(particleDataPtr, coupSMPtr); + hmeTau2PhaseSpace .initPointers(particleDataPtr, coupSMPtr); // User selected tau settings. - tauExt = settingsPtr->mode("TauDecays:externalMode"); - tauMode = settingsPtr->mode("TauDecays:mode"); - tauMother = settingsPtr->mode("TauDecays:tauMother"); - tauPol = settingsPtr->parm("TauDecays:tauPolarization"); + tauExt = mode("TauDecays:externalMode"); + tauMode = mode("TauDecays:mode"); + tauMother = mode("TauDecays:tauMother"); + tauPol = parm("TauDecays:tauPolarization"); // Parameters to determine if correlated partner should decay. - limitTau0 = settingsPtr->flag("ParticleDecays:limitTau0"); - tau0Max = settingsPtr->parm("ParticleDecays:tau0Max"); - limitTau = settingsPtr->flag("ParticleDecays:limitTau"); - tauMax = settingsPtr->parm("ParticleDecays:tauMax"); - limitRadius = settingsPtr->flag("ParticleDecays:limitRadius"); - rMax = settingsPtr->parm("ParticleDecays:rMax"); - limitCylinder = settingsPtr->flag("ParticleDecays:limitCylinder"); - xyMax = settingsPtr->parm("ParticleDecays:xyMax"); - zMax = settingsPtr->parm("ParticleDecays:zMax"); + limitTau0 = flag("ParticleDecays:limitTau0"); + tau0Max = parm("ParticleDecays:tau0Max"); + limitTau = flag("ParticleDecays:limitTau"); + tauMax = parm("ParticleDecays:tauMax"); + limitRadius = flag("ParticleDecays:limitRadius"); + rMax = parm("ParticleDecays:rMax"); + limitCylinder = flag("ParticleDecays:limitCylinder"); + xyMax = parm("ParticleDecays:xyMax"); + zMax = parm("ParticleDecays:zMax"); limitDecay = limitTau0 || limitTau || limitRadius || limitCylinder; } @@ -130,10 +123,15 @@ bool TauDecays::decay(int idxOut1, Event& event) { int idxOut2 = event[idxOut2Top].iBotCopyId(); out2 = HelicityParticle(event[idxOut2]); - // Set the mediator of the hard process. + // Set the mediator of the hard process (also handle no mediator). int idxMediator = event[idxOut1Top].mother1(); + if (event[idxOut1Top].mother2() > event[idxOut1Top].mother1() && + event[idxOut1Top].mother1() == event[idxOut2Top].mother1() && + event[idxOut1Top].mother2() == event[idxOut2Top].mother2()) + idxMediator = idxOut1Top; mediator = HelicityParticle(event[idxMediator]); mediator.direction = -1; + if (idxMediator == idxOut1Top) mediator.id(23); if (mediator.m() < out1.m() + out2.m()) { Vec4 p = out1.p() + out2.p(); mediator.p(p); @@ -172,6 +170,12 @@ bool TauDecays::decay(int idxOut1, Event& event) { // Check partner can decay. else if (!out2.canDecay()) correlated = false; else if (!out2.mayDecay()) correlated = false; + // Check partner not EW showered, set decay matrix otherwise. + else if (!out2.isFinal() && out2.statusAbs() > 40 + && out2.statusAbs() < 60) { + correlated = false; + out2.D = out2.rho; + } } // Set the production mechanism. @@ -192,7 +196,7 @@ bool TauDecays::decay(int idxOut1, Event& event) { } } - // Catch unknown production mechanims. + // Catch unknown production mechanisms. if (!known) { particles[1] = mediator; if (abs(mediator.id()) == 22) @@ -258,7 +262,7 @@ bool TauDecays::decay(int idxOut1, Event& event) { // Switch the taus. tau = &particles[idx]; // Calculate second tau's density matrix. - hardME->calculateRho(idx, particles); + if (hardME) hardME->calculateRho(idx, particles); // Decay the second tau. children.clear(); @@ -296,29 +300,30 @@ bool TauDecays::decay(int idxOut1, Event& event) { bool TauDecays::internalMechanism(Event&) { - // Flag if process is known. - bool known = true; - // Produced from a photon, Z, or Z'. if (abs(mediator.id()) == 22 || abs(mediator.id()) == 23 || abs(mediator.id()) == 32) { + // Produced from photons: t-channel. + if (in1.id() == 22 && in2.id() == 22) { + hardME = hmeTwoGammas2TwoFermions.initChannel(particles); // Produced from fermions: s-channel. - if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18 && in1.daughter2() == 0 && - in2.daughter2() == 0 && in1.daughter1() == in2.daughter1()) { + } else if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18 && + in1.daughter1() == in2.daughter1() && + in1.daughter2() == in2.daughter2()) { particles.push_back(mediator); hardME = hmeTwoFermions2GammaZ2TwoFermions.initChannel(particles); // Unknown photon production. - } else known = false; + } else return false; // Produced from a W or W'. } else if (abs(mediator.id()) == 24 || abs(mediator.id()) == 34) { // Produced from fermions: s-channel. if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18 && in1.daughter2() == 0 && - in2.daughter2() == 0 && in1.daughter1() == in2.daughter1()) { + in2.daughter2() == 0 && in1.daughter1() == in2.daughter1()) { particles.push_back(mediator); hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); - // Unknown W production. - } else known = false; + // Unknown W production. + } else return false; // Produced from a Higgs. } else if (abs(mediator.id()) == 25 || abs(mediator.id()) == 35 || @@ -357,41 +362,29 @@ bool TauDecays::internalMechanism(Event&) { hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); // Unknown production. - } else known = false; - return known; + } else return false; + return true; } //-------------------------------------------------------------------------- // Determine the tau polarization and tau decay correlation using the provided -// SPINUP digits. +// SPINUP digits interpreted as helicity states. bool TauDecays::externalMechanism(Event &event) { - // Flag if process is known. - bool known = true; - - // Uncorrelated, take directly from SPINUP if valid. + // Uncorrelated, take directly from tau SPINUP if valid. if (tauExt == 0) correlated = false; if (!correlated) { - double spinup = particles[2].pol(); - if (abs(spinup) > 1.001) spinup = event[particles[2].iTopCopyId()].pol(); - if (abs(spinup) > 1.001) known = false; - else { - particles[2].rho[0][0] = (1 - spinup) / 2; - particles[2].rho[1][1] = (1 + spinup) / 2; - } + if (particles[2].pol() == 9) + particles[2].pol(event[particles[2].iTopCopyId()].pol()); + if (particles[2].pol() == 9) return false; - // Correlated, try mother. + // Correlated, take from mother SPINUP if valid. } else if (tauExt == 1) { - double spinup = mediator.pol(); - if (abs(spinup) > 1.001) spinup = event[mediator.iTopCopyId()].pol(); - if (abs(spinup) > 1.001) spinup = 0; - if (mediator.rho.size() > 1) { - mediator.rho[0][0] = (1 - spinup) / mediator.spinStates(); - mediator.rho[1][1] = (1 + spinup) / mediator.spinStates(); - } + if (mediator.pol() == 9) mediator.pol(event[mediator.iTopCopyId()].pol()); + if (mediator.pol() == 9) return false; particles[1] = mediator; if (abs(mediator.id()) == 22) hardME = hmeGamma2TwoFermions.initChannel(particles); @@ -402,11 +395,11 @@ bool TauDecays::externalMechanism(Event &event) { else if (abs(mediator.id()) == 25 || abs(mediator.id()) == 35 || abs(mediator.id()) == 36 || abs(mediator.id()) == 37) hardME = hmeHiggs2TwoFermions.initChannel(particles); - else known = false; + else return false; - // Unknown mechanism. - } else known = false; - return known; + // Unknown mechanism. + } else return false; + return true; } @@ -615,24 +608,11 @@ void TauDecays::isotropicDecay(vector& children) { // Perform two-particle decays in the respective rest frame. vector pInv(decayMult + 1); for (int i = 1; i < decayMult; ++i) { - double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - children[i].m()) - * (mInv[i] + mInv[i+1] + children[i].m()) - * (mInv[i] + mInv[i+1] - children[i].m()) - * (mInv[i] - mInv[i+1] + children[i].m()) ) / mInv[i]; - - // Isotropic angles give three-momentum. - double cosTheta = 2. * rndmPtr->flat() - 1.; - double sinTheta = sqrt(1. - cosTheta*cosTheta); - double phi = 2. * M_PI * rndmPtr->flat(); - double pX = pAbs * sinTheta * cos(phi); - double pY = pAbs * sinTheta * sin(phi); - double pZ = pAbs * cosTheta; - - // Calculate energies, fill four-momenta. - double eHad = sqrt( children[i].m()*children[i].m() + pAbs*pAbs); - double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs); - children[i].p( pX, pY, pZ, eHad); - pInv[i+1].p( -pX, -pY, -pZ, eInv); + // Fill four-momenta + pair ps = + rndmPtr->phaseSpace2(mInv[i], mInv[i+1], children[i].m()); + pInv[i+1].p(ps.first); + children[i].p(ps.second); } // Boost decay products to the mother rest frame. From 48cc2ee58ceba648709eefebe9bbd1eed56dcb9b Mon Sep 17 00:00:00 2001 From: Saptaparna Date: Mon, 24 Oct 2022 20:40:57 -0500 Subject: [PATCH 3/3] Further check to see if patch is really applied --- src/HelicityMatrixElements.cc | 46 +++++++++++++++++------------------ src/TauDecays.cc | 39 ++++++++++++++++------------- 2 files changed, 45 insertions(+), 40 deletions(-) diff --git a/src/HelicityMatrixElements.cc b/src/HelicityMatrixElements.cc index 7bfadf66b44..8c11cea9b9b 100644 --- a/src/HelicityMatrixElements.cc +++ b/src/HelicityMatrixElements.cc @@ -688,15 +688,15 @@ complex HMETwoGammas2TwoFermions::calculateME(vector h) { void HMEX2TwoFermions::initWaves(vector& p) { u.clear(); - pMap.resize(4); + pMap.resize(3); // Initialize boson wave function. - vector< Wave4 > u1; - pMap[1] = 1; - for (int h = 0; h < p[pMap[1]].spinStates(); h++) - u1.push_back(p[pMap[1]].wave(h)); - u.push_back(u1); + vector< Wave4 > u0; + pMap[0] = 0; + for (int h = 0; h < p[pMap[0]].spinStates(); h++) + u0.push_back(p[pMap[0]].waveBar(h)); + u.push_back(u0); // Initialize fermion wave functions. - setFermionLine(2, p[2], p[3]); + setFermionLine(1, p[1], p[2]); } @@ -739,8 +739,8 @@ complex HMEW2TwoFermions::calculateME(vector h) { complex answer(0,0); for (int mu = 0; mu <= 3; mu++) { answer += - u[0][h[pMap[1]]](mu) * (u[2][h[pMap[3]]] * gamma[mu] - * (p2CV + p2CA * gamma[5]) * u[1][h[pMap[2]]]); + u[0][h[pMap[0]]](mu) * (u[2][h[pMap[2]]] * gamma[mu] + * (p2CV + p2CA * gamma[5]) * u[1][h[pMap[1]]]); } return answer; @@ -762,7 +762,7 @@ complex HMEGamma2TwoFermions::calculateME(vector h) { complex answer(0,0); for (int mu = 0; mu <= 3; mu++) { answer += - u[0][h[pMap[1]]](mu) * (u[2][h[pMap[3]]] * gamma[mu] * u[1][h[pMap[2]]]); + u[0][h[pMap[0]]](mu) * (u[2][h[pMap[2]]] * gamma[mu] * u[1][h[pMap[1]]]); } return answer; } @@ -802,8 +802,8 @@ complex HMEZ2TwoFermions::calculateME(vector h) { complex answer(0,0); for (int mu = 0; mu <= 3; mu++) { answer += - u[0][h[pMap[1]]](mu) * (u[2][h[pMap[3]]] * gamma[mu] - * (p2CV - p2CA * gamma[5]) * u[1][h[pMap[2]]]); + u[0][h[pMap[0]]](mu) * (u[2][h[pMap[2]]] * gamma[mu] + * (p2CV - p2CA * gamma[5]) * u[1][h[pMap[1]]]); } return answer; } @@ -870,15 +870,15 @@ void HMEHiggs2TwoFermions::initConstants() { // Set the H4 constants. p2CA = 0; p2CV = 0; - if (abs(pID[1]) == 37) { - p2CA = pID[1] == 37 ? 1 : -1; p2CV = 1; + if (abs(pID[0]) == 37) { + p2CA = pID[0] == 37 ? 1 : -1; p2CV = 1; // Neutral constants; settings available. } else if (settingsPtr) { int mode; double eta, phi; // Set the H1 mixing. - if (abs(pID[1]) == 25) { + if (abs(pID[0]) == 25) { mode = settingsPtr->mode("HiggsH1:parity"); eta = settingsPtr->parm("HiggsH1:etaParity"); phi = settingsPtr->parm("HiggsH1:phiParity"); @@ -887,7 +887,7 @@ void HMEHiggs2TwoFermions::initConstants() { else if (mode == 4) {p2CA = cos(phi); p2CV = complex(0, 1) * sin(phi);} else {p2CA = 0; p2CV = complex(0, 1);} // Set the H2 mixing. - } else if (abs(pID[1]) == 35) { + } else if (abs(pID[0]) == 35) { mode = settingsPtr->mode("HiggsH2:parity"); eta = settingsPtr->parm("HiggsH2:etaParity"); phi = settingsPtr->parm("HiggsH2:phiParity"); @@ -896,7 +896,7 @@ void HMEHiggs2TwoFermions::initConstants() { else if (mode == 4) {p2CA = cos(phi); p2CV = complex(0, 1) * sin(phi);} else {p2CA = 0; p2CV = complex(0, 1);} // Set the A3 mixing. - } else if (abs(pID[1]) == 36) { + } else if (abs(pID[0]) == 36) { mode = settingsPtr->mode("HiggsA3:parity"); eta = settingsPtr->parm("HiggsA3:etaParity"); phi = settingsPtr->parm("HiggsA3:phiParity"); @@ -908,9 +908,9 @@ void HMEHiggs2TwoFermions::initConstants() { // Neutral constants; default SM/MSSM. } else { - if (abs(pID[1]) == 25) {p2CA = 0; p2CV = complex(0, 1);} - else if (abs(pID[1]) == 35) {p2CA = 0; p2CV = complex(0, 1);} - else if (abs(pID[1]) == 36) {p2CA = 1; p2CV = 0;} + if (abs(pID[0]) == 25) {p2CA = 0; p2CV = complex(0, 1);} + else if (abs(pID[0]) == 35) {p2CA = 0; p2CV = complex(0, 1);} + else if (abs(pID[0]) == 36) {p2CA = 1; p2CV = 0;} } } @@ -921,8 +921,8 @@ void HMEHiggs2TwoFermions::initConstants() { void HMEHiggs2TwoFermions::initWaves(vector& p) { u.clear(); - pMap.resize(4); - setFermionLine(2, p[2], p[3]); + pMap.resize(3); + setFermionLine(1, p[1], p[2]); } @@ -932,7 +932,7 @@ void HMEHiggs2TwoFermions::initWaves(vector& p) { complex HMEHiggs2TwoFermions::calculateME(vector h) { - return (u[1][h[pMap[3]]] * (p2CV + p2CA * gamma[5]) * u[0][h[pMap[2]]]); + return (u[1][h[pMap[2]]] * (p2CV + p2CA * gamma[5]) * u[0][h[pMap[1]]]); } diff --git a/src/TauDecays.cc b/src/TauDecays.cc index cdeaf6f9adb..261ef0fc9d6 100644 --- a/src/TauDecays.cc +++ b/src/TauDecays.cc @@ -124,14 +124,16 @@ bool TauDecays::decay(int idxOut1, Event& event) { out2 = HelicityParticle(event[idxOut2]); // Set the mediator of the hard process (also handle no mediator). - int idxMediator = event[idxOut1Top].mother1(); + int idxMediator = event[idxOut1Top].mother1(); if (event[idxOut1Top].mother2() > event[idxOut1Top].mother1() && event[idxOut1Top].mother1() == event[idxOut2Top].mother1() && event[idxOut1Top].mother2() == event[idxOut2Top].mother2()) idxMediator = idxOut1Top; - mediator = HelicityParticle(event[idxMediator]); - mediator.direction = -1; - if (idxMediator == idxOut1Top) mediator.id(23); + mediator = HelicityParticle(event[idxMediator]); + if (idxMediator == idxOut1Top) { + if (out1.id() == -out2.id()) mediator.id(23); + else mediator.id(out1.id() < 0 ? 24 : -24); + } if (mediator.m() < out1.m() + out2.m()) { Vec4 p = out1.p() + out2.p(); mediator.p(p); @@ -143,10 +145,10 @@ bool TauDecays::decay(int idxOut1, Event& event) { int idxIn1 = event[idxMediatorTop].mother1(); int idxIn2 = event[idxMediatorTop].mother2(); in1 = HelicityParticle(event[idxIn1]); - in1.direction = -1; in2 = HelicityParticle(event[idxIn2]); + in1.direction = -1; in2.direction = -1; - + // Set the particles vector. particles.clear(); particles.push_back(in1); @@ -196,9 +198,11 @@ bool TauDecays::decay(int idxOut1, Event& event) { } } - // Catch unknown production mechanisms. + // Unknown production mechanisms, treat as vector boson decay. if (!known) { - particles[1] = mediator; + particles.erase(particles.begin()); + particles[0] = HelicityParticle(mediator.id(), 0, 0, 0, 0, 0, 0, 0, + mediator.p(), mediator.m(), 0, particleDataPtr); if (abs(mediator.id()) == 22) hardME = hmeGamma2TwoFermions.initChannel(particles); else if (abs(mediator.id()) == 23 || abs(mediator.id()) == 32) @@ -207,7 +211,7 @@ bool TauDecays::decay(int idxOut1, Event& event) { hardME = hmeW2TwoFermions.initChannel(particles); else if (correlated) { Vec4 p = out1.p() + out2.p(); - particles[1] = HelicityParticle(22, -22, idxIn1, idxIn2, idxOut1, + particles[0] = HelicityParticle(22, -22, idxIn1, idxIn2, idxOut1, idxOut2, 0, 0, p, p.mCalc(), 0, particleDataPtr); hardME = hmeGamma2TwoFermions.initChannel(particles); infoPtr->errorMsg("Warning in TauDecays::decay: unknown correlated " @@ -223,12 +227,13 @@ bool TauDecays::decay(int idxOut1, Event& event) { // Pick the first tau to decay. HelicityParticle* tau; - int idx = 2; - if (correlated) idx = (rndmPtr->flat() < 0.5) ? 2 : 3; - tau = &particles[idx]; + int idx1 = particles.size() == 3 ? 1 : 2; + int idx2 = idx1 + 1; + if (correlated && rndmPtr->flat() < 0.5) swap(idx1, idx2); + tau = &particles[idx1]; // Calculate the density matrix (if needed) and select channel. - if (hardME) hardME->calculateRho(idx, particles); + if (hardME) hardME->calculateRho(idx1, particles); vector children = createChildren(*tau); if (children.size() == 0) return false; @@ -254,15 +259,14 @@ bool TauDecays::decay(int idxOut1, Event& event) { // If a correlated second tau exists, decay that tau as well. if (correlated) { - idx = (idx == 2) ? 3 : 2; // Calculate the first tau decay matrix. decayME->calculateD(children); // Update the decay matrix for the tau. tau->D = children[0].D; // Switch the taus. - tau = &particles[idx]; + tau = &particles[idx2]; // Calculate second tau's density matrix. - if (hardME) hardME->calculateRho(idx, particles); + if (hardME) hardME->calculateRho(idx2, particles); // Decay the second tau. children.clear(); @@ -328,7 +332,8 @@ bool TauDecays::internalMechanism(Event&) { // Produced from a Higgs. } else if (abs(mediator.id()) == 25 || abs(mediator.id()) == 35 || abs(mediator.id()) == 36 || abs(mediator.id()) == 37) { - particles[1] = mediator; + particles.erase(particles.begin()); + particles[0] = mediator; hardME = hmeHiggs2TwoFermions.initChannel(particles); // Produced from a D or B hadron decay with a single tau.