Skip to content

Commit

Permalink
Merge pull request #29 from Saptaparna/cms/306
Browse files Browse the repository at this point in the history 
Pythia Update to fix tau decays
  • Loading branch information
@smuzaffar
smuzaffar authored Oct 25, 2022
2 parents e7384d4 + 48cc2ee commit 2859daf
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 42 deletions.
48 changes: 24 additions & 24 deletions src/HelicityMatrixElements.cc
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
// HelicityMatrixElements.cc is a part of the PYTHIA event generator.
// Copyright (C) 2021 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.

Expand Down Expand Up @@ -688,15 +688,15 @@ complex HMETwoGammas2TwoFermions::calculateME(vector<int> h) {
void HMEX2TwoFermions::initWaves(vector<HelicityParticle>& 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]);

}

Expand Down Expand Up @@ -739,8 +739,8 @@ complex HMEW2TwoFermions::calculateME(vector<int> 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;

Expand All @@ -762,7 +762,7 @@ complex HMEGamma2TwoFermions::calculateME(vector<int> 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;
}
Expand Down Expand Up @@ -802,8 +802,8 @@ complex HMEZ2TwoFermions::calculateME(vector<int> 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;
}
Expand Down Expand Up @@ -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");
Expand All @@ -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");
Expand All @@ -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");
Expand All @@ -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;}
}
}

Expand All @@ -921,8 +921,8 @@ void HMEHiggs2TwoFermions::initConstants() {
void HMEHiggs2TwoFermions::initWaves(vector<HelicityParticle>& p) {

u.clear();
pMap.resize(4);
setFermionLine(2, p[2], p[3]);
pMap.resize(3);
setFermionLine(1, p[1], p[2]);

}

Expand All @@ -932,7 +932,7 @@ void HMEHiggs2TwoFermions::initWaves(vector<HelicityParticle>& p) {

complex HMEHiggs2TwoFermions::calculateME(vector<int> 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]]]);

}

Expand Down
41 changes: 23 additions & 18 deletions src/TauDecays.cc
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
// TauDecays.cc is a part of the PYTHIA event generator.
// Copyright (C) 2021 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.

Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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)
Expand All @@ -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 "
Expand All @@ -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<HelicityParticle> children = createChildren(*tau);
if (children.size() == 0) return false;

Expand All @@ -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();
Expand Down Expand Up @@ -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.
Expand Down

0 comments on commit 2859daf

Please sign in to comment.