Rivet analyses
MC analyis of V → Pℓ+ℓ− and P → Vℓ+ℓ− decays
Experiment: ()
Status: VALIDATED
Authors: - Peter Richardson
References: none listed
Beams: * *
Beam energies: ANY
Run details: - Any process producing vector or scalar meson decays
A Monte Carlo analysis for the decay of vector mesons to a pseudoscalar meson, or a pseudoscalar meson to a vector meson, and an e+e− or μ+μ− pair. All such decays are automatically identified. Based on old Herwig++ internal analysis.
Source
code:MC_DECAY_MESON_MESON_LEPTONS.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/UnstableParticles.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
namespace Rivet {
/// @brief MC decay M -> M l+ l-
class MC_DECAY_MESON_MESON_LEPTONS : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_MESON_MESON_LEPTONS);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
declare(UnstableParticles(),"UFS");
// book the Histos
// pi0 dalitz
bookHistos(111,22,11,0.140);
// eta dalitz
bookHistos(221,22,11,0.55);
bookHistos(221,22,13,0.55);
// eta' dalitz
bookHistos(331,22,11,0.96);
bookHistos(331,22,13,0.96);
// omega -> pi0
bookHistos(223,111,11,0.8);
bookHistos(223,111,13,0.8);
// phi -> pi0
bookHistos(333,111,11,1.1);
bookHistos(333,111,13,1.1);
// phi -> eta
bookHistos(333,221,11,1.1);
bookHistos(333,221,13,1.1);
// J/psi dalitz
bookHistos(443,22,11,3.1);
bookHistos(443,22,13,3.1);
// B -> s gamma
bookHistos(511,313,11,5.3);
bookHistos(511,313,13,5.3);
}
void bookHistos(int id1, int id2, int il, double dM) {
if(abs(id2)%10==3 || id2==22) {
_incoming_P.push_back(id1);
_outgoingV.push_back(id2);
_outgoingf_P.push_back(il);
std::ostringstream title;
title << "h2_" << abs(id1);
if(id1>0) title << "p";
else title << "m";
title << "_" << abs(id2);
if(id2>0) title << "p";
else title << "m";
title << "_" << il << "_";
_mff_P .push_back(Histo1DPtr());
book(_mff_P.back(), title.str()+"mff" , 100, 0., dM);
_mVf .push_back(Histo1DPtr());
book(_mVf.back(), title.str()+"mVf" , 100, 0., dM);
_mVfbar.push_back(Histo1DPtr());
book(_mVfbar.back(), title.str()+"mVfbar", 100, 0., dM);
}
else {
_incomingV.push_back(id1);
_outgoingP.push_back(id2);
_outgoingf_V.push_back(il);
std::ostringstream title;
title << "h_" << abs(id1);
if(id1>0) title << "p";
else title << "m";
title << "_" << abs(id2);
if(id2>0) title << "p";
else title << "m";
title << "_" << il << "_";
_mff_V .push_back(Histo1DPtr());
book(_mff_V.back(), title.str()+"mff" , 100, 0., dM);
_mPf .push_back(Histo1DPtr());
book(_mPf.back(), title.str()+"mPf" , 100, 0., dM);
_mPfbar.push_back(Histo1DPtr());
book(_mPfbar.back(), title.str()+"mPfbar", 100, 0., dM);
}
}
void findDecayProducts(const Particle & mother,
unsigned int & nstable,
Particles& lp, Particles& lm,
Particles& scalar,
Particles& vector) {
for(const Particle & p : mother.children()) {
int id = p.pid();
if ( id == PID::EMINUS || id == PID::MUON ) {
lm.push_back(p);
++nstable;
}
else if (id == PID::EPLUS || id == PID::ANTIMUON) {
lp.push_back(p);
++nstable;
}
else if (abs(id)%10==1 && PID::isMeson(id)) {
scalar.push_back(p);
++nstable;
}
else if ((abs(id)%10==3 && PID::isMeson(id)) ||
id==PID::PHOTON ) {
vector.push_back(p);
++nstable;
}
else if ( !p.children().empty() ) {
findDecayProducts(p,nstable,lp,lm,scalar,vector);
}
else
++nstable;
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// loop over unstable particles
for(const Particle& iMeson : apply<UnstableParticles>(event, "UFS").particles()) {
// only consider scalar/vector mesons
long pid = iMeson.pid();
if(!PID::isMeson(pid)) continue;
if(abs(pid)%10!=3 and abs(pid)%10!=1 ) continue;
Particles lp,lm,scalar,vector;
unsigned int nstable(0);
findDecayProducts(iMeson,nstable,lp,lm,scalar,vector);
if(nstable!=3 || lp.size()!=1 || lm.size()!=1 || lp[0].pid()!=-lm[0].pid()) continue;
if(scalar.size()==1) {
// check if we already have this decay
unsigned int ix=0; bool found(false);
while(!found&&ix<_incomingV.size()) {
if(_incomingV[ix]==pid && _outgoingP[ix]==scalar[0].pid() &&
_outgoingf_V[ix]==lm[0].pid()) {
found=true;
}
else {
++ix;
}
}
// create a new graph if needed
if(!found) {
MSG_WARNING("MC_DECAY_MESON_MESON_LEPTONS S" << iMeson.pid() << " " << scalar[0].pid() << " "
<< iMeson.mass() << "\n");
continue;
}
// add the results to the histogram
_mff_V [ix]->fill((lm [0].momentum()+lp[0].momentum()).mass());
_mPf [ix]->fill((scalar[0].momentum()+lm[0].momentum()).mass());
_mPfbar[ix]->fill((scalar[0].momentum()+lp[0].momentum()).mass());
}
else if(vector.size()==1) {
// check if we already have this decay
unsigned int ix=0; bool found(false);
while(!found&&ix<_incoming_P.size()) {
if(_incoming_P[ix]==pid && _outgoingV[ix]==vector[0].pid() &&
_outgoingf_P[ix]==lm[0].pid()) {
found=true;
}
else {
++ix;
}
}
// create a new graph if needed
if(!found) {
MSG_WARNING("MC_DECAY_MESON_MESON_LEPTONS V" << iMeson.pid() << " " << vector[0].pid() << " "
<< iMeson.mass() << "\n");
continue;
}
// add the results to the histogram
_mff_P [ix]->fill((lm [0].momentum()+lp[0].momentum()).mass());
_mVf [ix]->fill((vector[0].momentum()+lm[0].momentum()).mass());
_mVfbar[ix]->fill((vector[0].momentum()+lp[0].momentum()).mass());
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
// normalize to unity V->P
for(unsigned int ix=0;ix<_mff_V.size();++ix) {
normalize(_mff_V);
normalize(_mPf);
normalize(_mPfbar);
}
// normalize to unity P->V
for(unsigned int ix=0;ix<_mff_P.size();++ix) {
normalize(_mff_P);
normalize(_mVf);
normalize(_mVfbar);
}
}
/// @}
/// @name Histograms for V -> P
/// @{
/**
* PDG codes of the incoming particles
*/
vector<long> _incomingV;
/**
* PDG codes of the outgoing pseudoscalar mesons
*/
vector<long> _outgoingP;
/**
* PDG codes of the outgoing fermion
*/
vector<long> _outgoingf_V;
/**
* Histograms for the mass of the fermion-antifermion pair
*/
vector<Histo1DPtr> _mff_V;
/**
* Histograms for the masses of the pseudoscalar and the fermion
*/
vector<Histo1DPtr> _mPf;
/**
* Histograms for the masses of the pseudoscalar and the antifermion
*/
vector<Histo1DPtr> _mPfbar;
/// @}
/// @name Histograms P->V
/// @{
/**
* PDG codes of the incoming_P particles
*/
vector<long> _incoming_P;
/**
* PDG codes of the outgoing vector mesons
*/
vector<long> _outgoingV;
/**
* PDG codes of the outgoing fermion
*/
vector<long> _outgoingf_P;
/**
* Histograms for the mass of the fermion-antifermion pair
*/
vector<Histo1DPtr> _mff_P;
/**
* Histograms for the masses of the vector and the fermion
*/
vector<Histo1DPtr> _mVf;
/**
* Histograms for the masses of the vector and the antifermion
*/
vector<Histo1DPtr> _mVfbar;
/// @}
};
RIVET_DECLARE_PLUGIN(MC_DECAY_MESON_MESON_LEPTONS);
}