// -*- 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"#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
#include"Rivet/Projections/InitialQuarks.hh"namespace Rivet {
/// @brief pi0 spectrum
classDELPHI_1996_I401100:public Analysis {
public:/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1996_I401100);
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
voidinit() {
declare(Beam(), "Beams");
declare(ChargedFinalState(), "FS");
declare(UnstableParticles(), "UFS");
declare(InitialQuarks(), "IQF");
// Book histograms
book(_h_pi_all, 1, 1, 1);
book(_h_pi_bot, 3, 1, 1);
book(_wAll,"TMP/wAll");
book(_wBot,"TMP/wBot");
}
/// Perform the per-event analysis
voidanalyze(const Event& event) {
// First, veto on leptonic events by requiring at least 4 charged FS particles
const FinalState& fs = apply<FinalState>(event, "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");
int flavour =0;
const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
// If we only have two quarks (qqbar), just take the flavour.
// If we have more than two quarks, look for the highest energetic q-qbar pair.
if (iqf.particles().size() ==2) {
flavour = iqf.particles().front().abspid();
}
else {
map<int, double> quarkmap;
for (const Particle&p : iqf.particles()) {
if (quarkmap[p.pid()] < p.E()) {
quarkmap[p.pid()] = p.E();
}
}
double maxenergy =0.;
for (int i =1; i <=5; ++i) {
if (quarkmap[i]+quarkmap[-i] > maxenergy) {
flavour = i;
}
}
}
_wAll->fill();
if(flavour==5) _wBot->fill();
// Get beams and average beam momentum
const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
constdouble meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) /2.0;
MSG_DEBUG("Avg beam momentum = "<< meanBeamMom);
// Final state of unstable particles to get particle spectra
const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
for (const Particle&p : ufs.particles(Cuts::pid==PID::PI0)) {
double xp = p.p3().mod()/meanBeamMom;
_h_pi_all->fill(xp);
if(flavour==5) _h_pi_bot->fill(xp);
}
}
/// Normalise histograms etc., after the run
voidfinalize() {
scale(_h_pi_all, 1./*_wAll);
scale(_h_pi_bot, 1./*_wBot);
}
//@}
/// @name Histograms
//@{
Histo1DPtr _h_pi_all, _h_pi_bot;
CounterPtr _wAll,_wBot;
//@}
};
// The hook for the plugin system
RIVET_DECLARE_PLUGIN(DELPHI_1996_I401100);
}