Rivet analyses
K*0 meson production measured by OPAL at LEP 1.
Experiment: OPAL (LEP 1)
Inspire ID: 447146
Status: VALIDATED
Authors: - Peter Richardson
References: - Phys.Lett.B412:210-224,1997 - hep-ex/9708022
Beams: e+ e-
Beam energies: (45.6, 45.6)GeV
Run details: - Hadronic Z decay events generated on the Z pole ($\sqrt{s} = 91.2$ GeV)
The K*0 fragmentation function has been measured in hadronic Z0 decays. In addition the helicity density matrix elements for inclusive K*(892)0 mesons from hadronic Z0 decays have been measured over the full range of K*0 momentum using data taken with the OPAL experiment at LEP. Both the fragmentation function and polarization measurements are implemented.
Source
code:OPAL_1997_I447146.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief OPAL K*0 fragmentation function paper
///
/// @author Peter Richardson
class OPAL_1997_I447146 : public Analysis {
public:
RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1997_I447146);
/// @name Analysis methods
/// @{
void init() {
declare(Beam(), "Beams");
declare(ChargedFinalState(), "FS");
declare(UnstableParticles(), "UFS");
book(_histXeK0, 1, 1, 1);
const vector<double> edges{0., 0.01, 0.03, 0.1, 0.125, 0.14, 0.16, 0.2, 0.3, 0.4, 0.5, 0.7, 1.0};
book(_h_ctheta, edges);
book(_h_alpha, edges);
for (size_t ix=0; ix < _h_ctheta->numBins(); ++ix) {
string suff = to_string(ix);
if (suff.length() == 1) suff = "0"+suff;
book(_h_ctheta->bin(ix+1), "ctheta_"+suff, 20, -1.0, 1.0);
book(_h_alpha->bin(ix+1), "alpha_"+suff, 20, 0.0, 0.5*M_PI);
}
book(_h_ctheta_large,"ctheta_large", 20,-1.0, 1.0);
book(_h_alpha_low, {0.3, 0.4, 0.5, 0.7, 1.0});
book(_h_alpha_high, {0.3, 0.4, 0.5, 0.7, 1.0});
for (size_t ix=0; ix < _h_alpha_low->numBins(); ++ix) {
string suff = std::to_string(ix);
book(_h_alpha_low->bin(ix+1), "alpha_low_"+suff, 20, 0.0, 0.5*M_PI);
book(_h_alpha_high->bin(ix+1), "alpha_high_"+suff, 20, 0.0, 0.5*M_PI);
}
book(_h_alpha_large ,"_h_alpha_large" ,20,0.,0.5*M_PI);
book(_h_alpha_large_low ,"_h_alpha_large_low" ,20,0.,0.5*M_PI);
book(_h_alpha_large_high,"_h_alpha_large_high",20,0.,0.5*M_PI);
}
pair<double,double> calcRho(Histo1DPtr hist,unsigned int imode) {
if (hist->numEntries()==0.) return make_pair(0.,0.);
double sum1(0.),sum2(0.);
for (const auto& bin : hist->bins() ) {
double Oi = bin.sumW();
if(Oi==0.) continue;
double ai,bi;
if(imode==0) {
ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
}
else {
ai = 2.*(bin.xMax()-bin.xMin())/M_PI;
bi = 2.*(sin(2.*bin.xMax())-sin(2.*bin.xMin()))/M_PI;
}
double Ei = bin.errW();
sum1 += sqr(bi/Ei);
sum2 += bi/sqr(Ei)*(Oi-ai);
}
return make_pair(sum2/sum1,sqrt(1./sum1));
}
void analyze(const Event& e) {
// First, veto on leptonic events by requiring at least 4 charged FS particles
const FinalState& fs = apply<FinalState>(e, "FS");
const size_t numParticles = fs.particles().size();
// Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
if (numParticles < 2) {
MSG_DEBUG("Failed leptonic event cut");
vetoEvent;
}
MSG_DEBUG("Passed leptonic event cut");
// Get beams and average beam momentum
const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
const double meanBeamMom = 0.5*(beams.first.p3().mod() + beams.second.p3().mod());
MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
Vector3 axis;
if (beams.first.pid()>0) {
axis = beams.first .momentum().p3().unit();
}
else {
axis = beams.second.momentum().p3().unit();
}
// Final state of unstable particles to get particle spectra
const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
for (const Particle& p : ufs.particles(Cuts::abspid==313)) {
double xp = p.p3().mod()/meanBeamMom;
_histXeK0->fill(xp);
int sign = p.pid()/313;
if(p.children().size()!=2) continue;
Particle kaon;
if (p.children()[0].pid()==sign*321 && p.children()[1].pid()==-sign*211) {
kaon = p.children()[0];
}
else if (p.children()[1].pid()==sign*321 && p.children()[0].pid()==-sign*211) {
kaon = p.children()[1];
}
else continue;
// spin axes
Vector3 e1z = p.momentum().p3().unit();
Vector3 e1y = e1z.cross(axis).unit();
Vector3 e1x = e1y.cross(e1z).unit();
LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
Vector3 axis1 = boost.transform(kaon.momentum()).p3().unit();
double ctheta = e1z.dot(axis1);
double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
double alpha = abs(abs(phi)-0.5*M_PI);
_h_ctheta->fill(xp,ctheta);
_h_alpha->fill(xp,alpha );
double cBeam = axis.dot(e1z);
(cBeam < 0.5 ? _h_alpha_low : _h_alpha_high)->fill(xp,alpha);
if (xp > 0.3) {
_h_ctheta_large->fill(ctheta);
_h_alpha_large->fill(alpha);
(cBeam < 0.5 ? _h_alpha_large_low : _h_alpha_large_high)->fill(alpha);
}
}
}
/// Finalize
void finalize() {
scale(_histXeK0, 1./sumOfWeights());
const vector<double> x = {0., 0.01 ,0.03 ,0.10 ,0.125,0.14 ,0.16 ,0.20 ,0.30 ,0.40 ,0.50 ,0.70 ,1.00};
Estimate1DPtr h_rho00;
book(h_rho00,2,1,1);
Estimate1DPtr h_rho_off;
book(h_rho_off,2,1,2);
Estimate1DPtr h_ratio;
book(h_ratio,3,1,1);
Estimate1DPtr h_off_low;
book(h_off_low,4,1,1);
Estimate1DPtr h_off_high;
book(h_off_high,4,1,2);
for (size_t ix=0; ix<_h_ctheta->numBins(); ++ix) {
// extract the rho00 component
normalize(_h_ctheta->bin(ix+1));
pair<double,double> rho00 = calcRho(_h_ctheta->bin(ix+1), 0);
h_rho00->bin(ix+1).set(rho00.first, rho00.second);
// extract the rho 1 -1 component
normalize(_h_alpha->bin(ix+1));
pair<double,double> rho_off = calcRho(_h_alpha->bin(ix+1), 1);
h_rho_off->bin(ix+1).set(rho_off.first, rho_off.second);
if (ix<8) continue;
// at large xp also the ratio
size_t iy = ix-8;
double ratio = rho_off.first/(1.-rho00.first);
double dr = ((rho_off.second - rho00.first*rho_off.second + rho_off.first*rho00.second))/sqr(1.-rho00.first);
h_ratio->bin(iy+1).set(ratio, dr);
// rho 1 -1 for cos theta <0.5
normalize(_h_alpha_low->bin(iy+1));
rho_off = calcRho(_h_alpha_low->bin(iy+1), 1);
h_off_low->bin(iy+1).set(rho_off.first, rho_off.second);
// rho 1 -1 for cos theta >0.5
normalize(_h_alpha_high->bin(iy+1));
rho_off = calcRho(_h_alpha_high->bin(iy+1), 1);
h_off_high->bin(iy+1).set(rho_off.first, rho_off.second);
}
// ratio for xp 0.3
normalize(_h_ctheta_large);
pair<double,double> rho00 = calcRho(_h_ctheta_large,0);
normalize(_h_alpha_large);
pair<double,double> rho_off = calcRho(_h_alpha_large,1);
double ratio = rho_off.first/(1.-rho00.first);
double dr = ((rho_off.second - rho00.first*rho_off.second + rho_off.first*rho00.second))/sqr(1.-rho00.first);
Estimate1DPtr h_ratio_large;
book(h_ratio_large,3,2,1);
h_ratio_large->bin(1).set(ratio, dr);
// rho 1 -1 for xp >0.3 and cos theta < 0.5
normalize(_h_alpha_large_low );
rho_off = calcRho(_h_alpha_large_low,1);
book(h_off_low,4,2,1);
h_off_low->bin(1).set(rho_off.first, rho_off.second);
// rho 1 -1 for xp >0.3 and cos theta > 0.5
normalize(_h_alpha_large_high);
rho_off = calcRho(_h_alpha_large_high,1);
book(h_off_high,4,2,2);
h_off_high->bin(1).set(rho_off.first, rho_off.second);
}
/// @}
private:
/// @{
Histo1DPtr _histXeK0;
Histo1DGroupPtr _h_ctheta, _h_alpha, _h_alpha_low, _h_alpha_high;
Histo1DPtr _h_ctheta_large, _h_alpha_large, _h_alpha_large_low, _h_alpha_large_high;
/// @}
};
RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1997_I447146, OPAL_1997_S3608263);
}Aliases: - OPAL_1997_S3608263