Rivet analyses
Charged particle asymmetries in e+e− → μ+μ−γ and e+e− → π+π−γ
Experiment: BABAR (PEP-II)
Inspire ID: 1388182
Status: VALIDATED NOHEPDATA SINGLEWEIGHT
Authors: - Peter Richardson
References: - Phys.Rev.D 92 (2015) 7, 072015
Beams: e+ e-
Beam energies: (3.5, 8.0); (5.3, 5.3)GeV
Run details: - e+e-> mu+ mu- gamma or pi+pi-gamma
Measurement of Charged particle asymmetries in e+e− → μ+μ−γ and e+e− → π+π−γ which are senistive to interference between the initial- and final-state QED radiation
Source
code:BABAR_2015_I1388182.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
namespace Rivet {
/// @brief e+ e- > mu+ mu- gamma or pi+ pi- gamma
class BABAR_2015_I1388182 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2015_I1388182);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
declare(Beam(), "Beams");
declare(FinalState(), "FS");
// histograms
Profile1DPtr tmp;
book(tmp, "TMP/pi0", refData(3, 1, 1));
_h_pipi.push_back(tmp);
for (unsigned int ix = 0; ix < 15; ++ix) {
book(tmp, 3, 1, 1 + ix);
_h_pipi.push_back(tmp);
if (ix > 13) continue;
book(tmp, 1, 1, 1 + ix);
_h_mumu.push_back(tmp);
}
book(tmp, "TMP/pi16", refData(3, 1, 15));
_h_pipi.push_back(tmp);
book(tmp, "TMP/pi17", refData(3, 1, 15));
_h_pipi.push_back(tmp);
book(_h_nopsi, "TMP/nopsi", refData(1, 1, 7));
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// get the axis, direction of incoming positron
const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
bool CMF = fuzzyEquals(-beams.first.mom().z() / beams.second.mom().z(), 1);
Vector3 axis1 = beams.first.mom().p3().unit();
Vector3 axis2 = beams.second.mom().p3().unit();
if (beams.first.pid() < 0) swap(axis1, axis2);
// find the final state final state particles
Particle mup, mum, pip, pim, gamma;
Particles fs = apply<FinalState>(event, "FS").particles();
// boost to CMF frame
LorentzTransform boost;
if (!CMF) {
boost = LorentzTransform::mkFrameTransformFromBeta(
(beams.first.mom() + beams.second.mom()).betaVec());
}
else {
double E1 = 3.5, E2 = 0.25 * sqr(sqrtS()) / E1;
FourMomentum pnew(E1 + E2, 0, 0, E1 - E2);
boost = LorentzTransform::mkFrameTransformFromBeta(-pnew.betaVec());
}
FourMomentum pGamma;
for (const Particle& p : fs) {
FourMomentum pLab, pCMF;
if (CMF) {
pCMF = p.mom();
pLab = boost.transform(p.mom());
}
else {
pCMF = boost.transform(p.mom());
pLab = p.mom();
}
double theta = acos(pLab.p3().unit().dot(axis1));
if (p.isCharged()) {
if (theta < .4 || theta > 2.45) continue;
if (pLab.p3().mod() < 1.) continue;
if (p.pid() == PID::MUON && mum.pid() != PID::MUON) {
mum = p;
}
else if (p.pid() == PID::ANTIMUON && mup.pid() != PID::ANTIMUON) {
mup = p;
}
else if (p.pid() == PID::PIPLUS && pip.pid() != PID::PIPLUS) {
pip = p;
}
else if (p.pid() == PID::PIMINUS && pim.pid() != PID::PIMINUS) {
pim = p;
}
}
else if (p.pid() == PID::GAMMA) {
// angle cut on the photon
if (theta < .35 || theta > 2.4) continue;
if (gamma.pid() != PID::GAMMA) {
gamma = p;
pGamma = pCMF;
}
else {
if (pCMF.E() > pGamma.E()) {
gamma = p;
pGamma = pCMF;
}
}
}
else {
vetoEvent;
}
}
if (gamma.pid() != PID::GAMMA) vetoEvent;
if (!(((pip.pid() == PID::PIPLUS && pim.pid() == PID::PIMINUS)
|| (mum.pid() == PID::MUON && mup.pid() == PID::ANTIMUON)))) {
vetoEvent;
}
if (pip.pid() == PID::PIPLUS && mum.pid() == PID::MUON) vetoEvent;
if (pGamma.E() < 3.) vetoEvent;
Vector3 axisZ = pGamma.p3().unit();
Vector3 axisX = (axis2 - axisZ.dot(axis2) * axisZ).unit();
Vector3 axisY = axisZ.cross(axisX);
FourMomentum pMinus, pPlus;
if (CMF) {
if (mum.pid() == PID::MUON) {
pMinus = mum.mom();
pPlus = mup.mom();
}
else {
pMinus = pim.mom();
pPlus = pip.mom();
}
}
else {
if (mum.pid() == PID::MUON) {
pMinus = boost.transform(mum.mom());
pPlus = boost.transform(mup.mom());
}
else {
pMinus = boost.transform(pim.mom());
pPlus = boost.transform(pip.mom());
}
}
double phiM = atan2(pMinus.p3().dot(axisY), pMinus.p3().dot(axisX));
if (phiM < 0.) phiM += 2. * M_PI;
double phiP = atan2(pPlus.p3().dot(axisY), pPlus.p3().dot(axisX));
if (phiP < 0.) phiP += 2. * M_PI;
double mass = (pMinus + pPlus).mass();
if (mum.pid() == PID::MUON) {
if (mass > 0.2 && mass < 7.) {
unsigned int imass = int(mass / .5);
if (phiM < M_PI)
_h_mumu[imass]->fill(cos(phiM), 1.);
else
_h_mumu[imass]->fill(cos(phiP), -1.);
if (imass == 6 && mass > 3.2) {
if (phiM < M_PI)
_h_nopsi->fill(cos(phiM), 1.);
else
_h_nopsi->fill(cos(phiP), -1.);
}
}
}
else if (pip.pid() == PID::PIPLUS) {
if (mass > 0.2 && mass < 2.0) {
unsigned int imass = int((mass - 0.2) / .1);
if (phiM < M_PI)
_h_pipi[imass]->fill(cos(phiM), 1.);
else
_h_pipi[imass]->fill(cos(phiP), -1.);
}
}
}
pair<double, double> calcA0(Profile1DPtr hist) {
if (hist->numEntries() == 0.) return make_pair(0., 0.);
double sum1(0.), sum2(0.);
for (const auto& bin : hist->bins()) {
if (bin.numEntries() < 2) continue;
double Oi = bin.mean(2);
double Ei = bin.stdErr(2);
if (Ei == 0.) continue;
double xi = 0.5 * (bin.xMin() + bin.xMax());
sum1 += sqr(xi / Ei);
sum2 += Oi * xi / sqr(Ei);
}
return make_pair(sum2 / sum1, sqrt(1. / sum1));
}
/// Normalise histograms etc., after the run
void finalize() {
// muon assymetries
Estimate1DPtr _A_mumu;
book(_A_mumu, 2, 1, 1);
for (unsigned int ix = 0; ix < 14; ++ix) {
pair<double, double> A0 = calcA0(ix != 6 ? _h_mumu[ix] : _h_nopsi);
_A_mumu->bin(ix + 1).set(A0.first, A0.second);
}
// pion assymetries
Estimate1DPtr _A_pipi;
book(_A_pipi, 4, 1, 1);
for (unsigned int ix = 0; ix < _h_pipi.size(); ++ix) {
pair<double, double> A0 = calcA0(_h_pipi[ix]);
_A_pipi->bin(ix + 1).set(A0.first, A0.second);
}
}
/// @}
/// @name Histograms
/// @{
vector<Profile1DPtr> _h_mumu, _h_pipi;
Profile1DPtr _h_nopsi;
/// @}
};
RIVET_DECLARE_PLUGIN(BABAR_2015_I1388182);
}