Rivet analyses
title: MC_TAUPOL
Monte Carlo validation observables for tau polarisation
Experiment: ()
Status: VALIDATED
Authors: - George Marsden - Christian Gutschow
References: none listed
Beams: * *
Beam energies: ANY
Run details: - single- or di-tau production
Observables sensitive to the polarisation of the $\tau$ lepton.
Source code:MC_TAUPOL.cc
```c++ // -- C++ --
include "Rivet/Analysis.hh"
include "Rivet/Projections/FinalState.hh"
include "Rivet/Projections/TauFinder.hh"
namespace Rivet {
/// @brief Monte Carlo validation observables for tau polarisation class MC_TAUPOL : public Analysis { public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(MC_TAUPOL);
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
_mode = 0; // default is ditau
if (getOption("MODE") == "DITAU") {
_mode = 0;
}
if (getOption("MODE") == "SINGLE") {
_mode = 1;
}
_lothresh = getOption("MASSMIN", _mode ? 75 * GeV : 86. * GeV);
_hithresh = getOption("MASSMAX", _mode ? 85 * GeV : 96 * GeV);
// Initialise and register projections
declare(TauFinder(TauDecay::ANY), "taus");
// Book histograms
// tau invmass
book(_h["m_tau"], "h_m_tau", 50, 1., 3.);
book(_h["m_tautau"], "h_m_tautau", 50, 35., 420.);
book(_h["m_B"], "h_m_B", 50, 35., 420.);
book(_h["B_pT"], "h_B_pT", 100, 0., 500.);
// longnitudal spin correlations
book(_h["m_tau_vis"], "h_m_tau_vis", 50, 0., 5.);
book(_h["m_tautau_vis"], "h_m_tautau_vis", 50, 0., 100.);
book(_h["m_pipi_vis"], "h_m_pipi_vis", 60, 0., 100.);
book(_h["m_pipi_vis_lm"], "h_m_pipi_vis_lm", 60, 0., 100.);
book(_h["m_pipi_vis_hm"], "h_m_pipi_vis_hm", 60, 0., 100.);
// transverse spin correlations
book(_h["acoplanarity_plus"], "h_acoplanarity_plus", 30, 0., 6.2);
book(_h["acop_pt40"], "h_acop_pt40", 30, 0., 6.2);
book(_h["acop_pt100"], "h_acop_pt100", 30, 0., 6.2);
book(_h["acop_pt200"], "h_acop_pt200", 30, 0., 6.2);
book(_h["acop_pt1000"], "h_acop_pt1000", 30, 0., 6.2);
book(_h["acom_pt40"], "h_acom_pt40", 30, 0., 6.2);
book(_h["acom_pt100"], "h_acom_pt100", 30, 0., 6.2);
book(_h["acom_pt200"], "h_acom_pt200", 30, 0., 6.2);
book(_h["acom_pt1000"], "h_acom_pt1000", 30, 0., 6.2);
book(_h["acoplanarity_minus"], "h_acoplanarity_minus", 30, 0., 6.2);
// leptonic
book(_h["2B_xlep"], "h_2B_xlep", 50, 0.0, 1.0);
book(_h["2B_xel"], "h_2B_xel", 50, 0.0, 1.0);
book(_h["2B_xmu"], "h_2B_xmu", 50, 0.0, 1.0);
book(_h["xel_B"], "h_xel_B", 50, 0.0, 1.0);
book(_h["xmu_B"], "h_xmu_B", 50, 0.0, 1.0);
// hadronic
book(_h["1B_xpi"], "h_1B_xpi", 30, 0.0, 1.0);
book(_h["1B_xpi_lm"], "h_1B_xpi_lm", 30, 0.0, 1.0);
book(_h["1B_xpi_hm"], "h_1B_xpi_hm", 30, 0.0, 1.0);
book(_h["xpi_B"], "h_xpi_B", 30, 0.0, 1.0);
book(_h["2B_xrho"], "h_2B_xrho", 50, 0.0, 1.0);
// angles
book(_h["t_pi"], "h_t_pi", 40, -1.0, 1.0);
book(_h["t_rho"], "h_t_rho", 40, -1.0, 1.0);
book(_h["pipi_theta"], "h_pipi_theta", 40, -1.0, 1.0);
book(_ntaus, "h_n_taus", {0, 1, 2, 3, 4});
}
void findDecayProducts(const Particle& mother,
size_t& nstable,
Particles& ep,
Particles& em,
Particles& nu_e,
Particles& nu_ebar,
Particles& mup,
Particles& mum,
Particles& nu_mu,
Particles& nu_mubar,
Particles& pip,
Particles& pim,
Particles& pi0,
Particles& Kp,
Particles& Km,
Particles& K0S,
Particles& K0L,
Particles& eta,
Particles& gamma) {
for (const Particle& p : mother.children()) {
int id = p.pid();
if (id == PID::KPLUS) {
Kp.push_back(p);
++nstable;
}
else if (id == PID::KMINUS) {
Km.push_back(p);
++nstable;
}
else if (id == PID::PIPLUS) {
pip.push_back(p);
++nstable;
}
else if (id == PID::PIMINUS) {
pim.push_back(p);
++nstable;
}
else if (id == PID::EPLUS) {
ep.push_back(p);
++nstable;
}
else if (id == PID::EMINUS) {
em.push_back(p);
++nstable;
}
else if (id == PID::NU_E) {
nu_e.push_back(p);
++nstable;
}
else if (id == PID::NU_EBAR) {
nu_ebar.push_back(p);
++nstable;
}
else if (id == PID::NU_MU) {
nu_mu.push_back(p);
++nstable;
}
else if (id == PID::NU_MUBAR) {
nu_mubar.push_back(p);
++nstable;
}
else if (id == PID::ANTIMUON) {
mup.push_back(p);
++nstable;
}
else if (id == PID::MUON) {
mum.push_back(p);
++nstable;
}
else if (id == PID::PI0) {
pi0.push_back(p);
++nstable;
}
else if (id == PID::K0S) {
K0S.push_back(p);
++nstable;
}
else if (id == PID::K0L) {
K0L.push_back(p);
++nstable;
}
else if (id == PID::ETA) {
eta.push_back(p);
++nstable;
}
else if (id == PID::PHOTON) {
gamma.push_back(p);
++nstable;
}
else if (!p.children().empty()) {
findDecayProducts(p, nstable, ep, em, nu_e, nu_ebar, mup, mum, nu_mu, nu_mubar, pip, pim, pi0, Kp,
Km, K0S, K0L, eta, gamma);
}
else {
++nstable;
}
}
}
double angleRF(const FourMomentum& p1, const FourMomentum& p2, const FourMomentum& RF) const {
LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(RF.betaVec());
return cos(boost(p1).angle(boost(p2)));
}
/// Perform the per-event analysis
void analyze(const Event& event) {
Particles taus, t_nu, pions;
for (const auto& tau : apply<TauFinder>(event, "taus").particlesByPt()) {
//remove any from hadronic decay (only want prompt taus)
if (tau.parents()[0].isPrompt()) {
taus += tau;
_h["m_tau"]->fill(tau.mom().mass() / GeV);
}
}
//veto ditau events without exactly 2 prompt taus
if (taus.size() != 2 && _mode == 0) vetoEvent;
//veto single-tau events without exactly 1 prompt tau
if (taus.size() != 1 && _mode == 1) vetoEvent;
FourMomentum tautaumom, taunumom, phomom, Bmom;
if (taus.size() == 2 && _mode == 0) {
tautaumom = taus[0].mom() + taus[1].mom();
phomom = sum(scanAncestors(taus, 22), Kin::p4, FourMomentum());
Bmom = tautaumom + phomom;
_h["m_tautau"]->fill(tautaumom.mass());
}
if (taus.size() == 1 && _mode == 1) {
t_nu = scanAncestors(taus, 16);
phomom = sum(scanAncestors(taus, 22), Kin::p4, FourMomentum());
if (t_nu.size() != 1) vetoEvent;
if (taus[0].pid() * t_nu[0].pid() > 0) vetoEvent;
taunumom = taus[0].mom() + t_nu[0].mom();
Bmom = taunumom + phomom;
}
_ntaus->fill(int(taus.size()));
_h["m_B"]->fill(Bmom.mass());
_h["B_pT"]->fill(Bmom.pT());
// Check opposite charges
// Calculate invariant mass of tautau
FourMomentum tt_mom_vis;
for (const auto& t : taus) {
FourMomentum t_mom_vis;
for (const auto& c : t.children()) {
if (c.isVisible()) {
tt_mom_vis += c.mom();
t_mom_vis += c.mom();
}
}
_h["m_tau_vis"]->fill(t_mom_vis.mass());
}
_h["m_tautau_vis"]->fill(tt_mom_vis.mass());
size_t OS = 2;
if (Bmom.mass() < _lothresh)
OS = 0;
else if (Bmom.mass() < _hithresh)
OS = 1;
for (const Particle& tau : taus) {
size_t nstable(0);
Particles ep, em, nu_e, nu_ebar, mup, mum, nu_mu, nu_mubar;
Particles pip, pim, pi0, Kp, Km, K0S, K0L, eta, gamma;
findDecayProducts(tau, nstable, ep, em, nu_e, nu_ebar, mup, mum, nu_mu, nu_mubar, pip, pim, pi0, Kp,
Km, K0S, K0L, eta, gamma);
if (tau.pid() < 0) {
swap(pim, pip);
swap(Kp, Km);
swap(em, ep);
swap(mum, mup);
swap(nu_e, nu_ebar);
swap(nu_mu, nu_mubar);
}
// classify by number of stable decay products
FourMomentum tmom = tau.mom();
const LorentzTransform tboost = LorentzTransform::mkFrameTransformFromBeta(tmom.betaVec());
if (nstable == 2) {
if (pim.size() == 1) {
const double xpi = pim[0].mom().p() / tau.mom().p();
const FourMomentum pimom = pim[0].mom();
const double angle = cos(tboost(pimom).angle(tmom));
_h["xpi_B"]->fill(pim[0].mom().p() / Bmom.p());
_h["t_pi"]->fill(angle);
if (OS == 0)
_h["1B_xpi_lm"]->fill(xpi);
else if (OS == 1)
_h["1B_xpi"]->fill(xpi);
else
_h["1B_xpi_hm"]->fill(xpi);
pions.push_back(pim[0]);
}
}
else if (nstable == 3) {
if (em.size() == 1 && nu_ebar.size() == 1) {
const double xlep = em[0].mom().p() / tau.mom().p();
_h["2B_xlep"]->fill(xlep);
_h["2B_xel"]->fill(xlep);
_h["xel_B"]->fill(em[0].mom().p() / Bmom.p());
}
else if (mum.size() == 1 && nu_mubar.size() == 1) {
const double xlep = mum[0].mom().p() / tau.mom().p();
_h["2B_xlep"]->fill(xlep);
_h["2B_xmu"]->fill(xlep);
_h["xmu_B"]->fill(mum[0].mom().p() / Bmom.p());
}
else if (pim.size() == 1 && pi0.size() == 1) {
const FourMomentum ptot = pim[0].mom() + pi0[0].mom();
const double xrho = ptot.p() / tau.mom().p();
_h["2B_xrho"]->fill(xrho);
const double angle = cos(tboost(ptot).angle(tmom));
_h["t_rho"]->fill(angle);
}
} //end nstable = 3
} //end tau loop
// TAU LONGNITUDAL SPIN EFFECTS PART
//fill tautau->pipi inv mass
if (pions.size() == 2) {
FourMomentum pi1mom = pions[0].mom();
FourMomentum pi2mom = pions[1].mom();
FourMomentum pipi = pions[0].mom() + pions[1].mom();
if (OS == 0)
_h["m_pipi_vis_lm"]->fill(pipi.mass() / GeV);
else if (OS == 1)
_h["m_pipi_vis"]->fill(pipi.mass() / GeV);
else
_h["m_pipi_vis_hm"]->fill(pipi.mass() / GeV);
_h["pipi_theta"]->fill(angleRF(pi1mom, pi2mom, Bmom));
}
// TAU TRANSVERSE SPIN EFFECTS PART
if (pions.size() == 2 && _mode != 1 && OS == 1) {
FourMomentum tautau = taus[0].mom() + taus[1].mom();
// calculate angles between (tau, tau_pi) and (tau2, tau2_pi) planes
const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(tautau.betaVec());
taus[0] = Particle(PID::ANY, boost.transform(taus[0].mom()));
pions[0] = Particle(PID::ANY, boost.transform(pions[0].mom()));
taus[1] = Particle(PID::ANY, boost.transform(taus[1].mom()));
pions[1] = Particle(PID::ANY, boost.transform(pions[1].mom()));
// calculate the angle
Vector3 t0xp0 = taus[0].p3().cross(pions[0].p3()).unit();
Vector3 t1xp1 = taus[1].p3().cross(pions[1].p3()).unit();
double phi_star_sign = acos(t0xp0.dot(t1xp1));
// extend acoplanarity angle to range (0, 2pi)
Vector3 p0xp1 = pions[0].p3().cross(pions[1].p3()).unit();
double sign = p0xp1.dot(taus[0].p3().unit());
if (sign > 0) phi_star_sign = 2 * M_PI - phi_star_sign;
// alpha angles in the lab frame
Vector3 tlv3_Ez(0.0, 0.0, 1.0);
Vector3 cross1_pi = tlv3_Ez.cross(taus[0].p3()).unit();
Vector3 cross2_pi = pions[0].p3().cross(taus[0].p3()).unit();
const double alpha_pi = acos(abs(cross1_pi.dot(cross2_pi)));
if (alpha_pi > M_PI / 4.) {
_h["acoplanarity_plus"]->fill(phi_star_sign);
if (Bmom.pT() < 40 * GeV)
_h["acop_pt40"]->fill(phi_star_sign);
else if (Bmom.pT() < 100 * GeV)
_h["acop_pt100"]->fill(phi_star_sign);
else if (Bmom.pT() < 200 * GeV)
_h["acop_pt200"]->fill(phi_star_sign);
else
_h["acop_pt1000"]->fill(phi_star_sign);
}
if (alpha_pi < M_PI / 4.) {
_h["acoplanarity_minus"]->fill(phi_star_sign);
if (Bmom.pT() < 40 * GeV)
_h["acom_pt40"]->fill(phi_star_sign);
else if (Bmom.pT() < 100 * GeV)
_h["acom_pt100"]->fill(phi_star_sign);
else if (Bmom.pT() < 200 * GeV)
_h["acom_pt200"]->fill(phi_star_sign);
else
_h["acom_pt1000"]->fill(phi_star_sign);
}
} // end of tau spin corr
} //end void function
/// Normalise histograms etc., after the run
void finalize() {
normalize(_h);
normalize(_ntaus);
}
//@}
Particles scanAncestors(const Particles& taus, int id) {
Particles rtn;
for (const Particle& tau : taus) {
Particle thistau = tau;
while (true) {
// check the tau has parents first else end loop
if (thistau.parents().empty()) break;
//loop over taus siblings
for (const Particle& sibling : thistau.parents()[0].children()) {
// find ID
if (sibling.abspid() == id) {
// add particles to rtn vector
bool duplicate = true;
for (const auto& p : rtn) {
if (sibling.isSame(p)) duplicate = false;
}
if (duplicate) rtn += sibling;
}
}
// want to loop over taus parent's siblings next
thistau = thistau.parents()[0];
// if taus parent is not a tau end the loop
if (thistau.abspid() != PID::TAU) break;
}
}
return rtn;
}
/// @name Histograms
//@{
// histograms for leptonic decay
map<string, Histo1DPtr> _h;
BinnedHistoPtr<int> _ntaus;
size_t _mode;
double _lothresh, _hithresh;
//@}
};
RIVET_DECLARE_PLUGIN(MC_TAUPOL); } ```