Rivet analyses

title: ALEPH_2000_I507531

Measurement of $\pi^0$, $\eta$, $\eta^\prime$, $K^0_S$ and $\Lambda^0$ spectra in two and three jet events

Experiment: ALEPH (LEP)

Inspire ID: 507531

Status: VALIDATED

Authors: - Peter Richardson

References: - Eur.Phys.J. C16 (2000) 613

Beams: e- e+

Beam energies: (45.6, 45.6)GeV

Run details: - e+e- -> hadrons

Measurement of $\pi^0$, $\eta$, $\eta^\prime$, $K^0_S$ and $\Lambda^0$ spectra in two and three jet events. In addition to the normal inclusive spectra the spectra in individual jets are measured for three jet events.

Source code:ALEPH_2000_I507531.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief pi, eta, eta', K0, lambda spectra class ALEPH_2000_I507531 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_2000_I507531);


/// @name Analysis methods
/// @{

/// Book histograms and initialise projections before the run
void init() {

  // Projections
  declare(Beam(), "Beams");
  declare(UnstableParticles(), "UFS");
  declare(FinalState(), "FS");

  // Histograms
  // incl
  book(_h_pi0, 1, 1, 1);
  book(_h_eta, 2, 1, 1);
  book(_h_etaP, 3, 1, 1);
  book(_h_K0, 16, 1, 1);
  book(_h_lam, 17, 1, 1);
  // two jet
  book(_h_2_pi0, 4, 1, 1);
  book(_h_2_eta, 5, 1, 1);
  book(_h_2_etaP, 6, 1, 1);
  book(_h_2_K0, 18, 1, 1);
  book(_h_2_lam, 19, 1, 1);
  // three jet
  book(_h_3_pi0[0], 7, 1, 1);
  book(_h_3_pi0[1], 8, 1, 1);
  book(_h_3_pi0[2], 9, 1, 1);
  book(_h_3_eta[0], 10, 1, 1);
  book(_h_3_eta[1], 11, 1, 1);
  book(_h_3_eta[2], 12, 1, 1);
  book(_h_3_etaP[0], 13, 1, 1);
  book(_h_3_etaP[1], 14, 1, 1);
  book(_h_3_etaP[2], 15, 1, 1);
  book(_h_3_K0[0], 20, 1, 1);
  book(_h_3_K0[1], 21, 1, 1);
  book(_h_3_K0[2], 22, 1, 1);
  book(_h_3_lam[0], 23, 1, 1);
  book(_h_3_lam[1], 24, 1, 1);
  book(_h_3_lam[2], 25, 1, 1);
  book(_w2, "/TMP/W2");
  book(_w3, "/TMP/W3");
}

void findDecayProducts(const Particle& parent, Particles& decay) {
  for (const Particle& child : parent.children()) {
    if (child.children().empty()) {
      decay.push_back(child);
    }
    else {
      findDecayProducts(child, decay);
    }
  }
}

/// Perform the per-event analysis
void analyze(const Event& event) {
  // Get beams and average beam momentum
  const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
  const double meanBeamMom = (beams.first.p3().mod() + beams.second.p3().mod()) / 2.0;

  Particles decay, fs;
  // unstable particles
  const UnstableParticles ufs = apply<UnstableParticles>(event, "UFS");
  for (const Particle& part : ufs.particles(Cuts::pid == 111 or Cuts::pid == 221 or Cuts::pid == 331
                                            or Cuts::pid == 310 or Cuts::abspid == 3122)) {
    fs.push_back(part);
    findDecayProducts(part, decay);
  }
  // FS particles
  for (const Particle& part : apply<FinalState>(event, "FS").particles()) {
    bool skip = false;
    for (const Particle& dec : decay) {
      if (dec.genParticle() == part.genParticle()) {
        skip = true;
        break;
      }
    }
    if (skip) continue;
    fs.push_back(part);
  }
  // Definition of the Durham algorithm
  fastjet::JetDefinition durham_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best);
  // pseudojets
  vector<fastjet::PseudoJet> input_particles;
  // Pseudo-jets from the non photons
  unsigned int ix = 0;
  for (const Particle& p : fs) {
    const FourMomentum p4 = p.momentum();
    input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
    input_particles.back().set_user_index(ix);
    ++ix;
  }
  // cluster the jets
  fastjet::ClusterSequence clust_seq(input_particles, durham_def);
  PseudoJets jets = fastjet::sorted_by_E(clust_seq.exclusive_jets_ycut(0.01));
  if (jets.size() == 2)
    _w2->fill();
  else if (jets.size() == 3)
    _w3->fill();
  ix = 0;
  for (const Particle& part : fs) {
    double xE = part.momentum().E() / meanBeamMom;
    double xP = part.momentum().p3().mod() / meanBeamMom;
    int ijet = jets.size() != 3 ? -1 : findJet(ix, jets);
    if (part.pid() == 111) {
      _h_pi0->fill(xE);
      if (jets.size() == 2) {
        _h_2_pi0->fill(xE);
      }
      else if (jets.size() == 3) {
        _h_3_pi0[ijet]->fill(xE);
      }
    }
    else if (part.pid() == 221) {
      _h_eta->fill(xE);
      if (jets.size() == 2) {
        _h_2_eta->fill(xE);
      }
      else if (jets.size() == 3) {
        _h_3_eta[ijet]->fill(xE);
      }
    }
    else if (part.pid() == 331) {
      _h_etaP->fill(xE);
      if (jets.size() == 2) {
        _h_2_etaP->fill(xE);
      }
      else if (jets.size() == 3) {
        _h_3_etaP[ijet]->fill(xE);
      }
    }
    else if (part.pid() == 310) {
      double xi = -log(xP);
      _h_K0->fill(xi);
      if (jets.size() == 2) {
        _h_2_K0->fill(xi);
      }
      else if (jets.size() == 3) {
        _h_3_K0[ijet]->fill(xi);
      }
    }
    else if (part.abspid() == 3122) {
      double xi = -log(xP);
      _h_lam->fill(xi);
      if (jets.size() == 2) {
        _h_2_lam->fill(xi);
      }
      else if (jets.size() == 3) {
        _h_3_lam[ijet]->fill(xi);
      }
    }
    else {
      break;
    }
    ix += 1;
  }
}

int findJet(int id, const PseudoJets& jets) {
  for (unsigned int ijet = 0; ijet < jets.size(); ++ijet) {
    for (const PseudoJet& con : jets[ijet].constituents()) {
      if (con.user_index() == id) return ijet;
    }
  }
  return -1;
}


/// Normalise histograms etc., after the run
void finalize() {
  scale(_h_pi0, 1. / sumOfWeights());
  scale(_h_eta, 1. / sumOfWeights());
  scale(_h_etaP, 1. / sumOfWeights());
  scale(_h_K0, 1. / sumOfWeights());
  scale(_h_lam, 1. / sumOfWeights());
  scale(_h_2_pi0, 1. / *_w2);
  scale(_h_2_eta, 1. / *_w2);
  scale(_h_2_etaP, 1. / *_w2);
  scale(_h_2_K0, 1. / *_w2);
  scale(_h_2_lam, 1. / *_w2);
  for (unsigned int ix = 0; ix < 3; ++ix) {
    scale(_h_3_pi0[ix], 1. / *_w3);
    scale(_h_3_eta[ix], 1. / *_w3);
    scale(_h_3_etaP[ix], 1. / *_w3);
    scale(_h_3_K0[ix], 1. / *_w3);
    scale(_h_3_lam[ix], 1. / *_w3);
  }
}

/// @}


/// @name Histograms
/// @{
Histo1DPtr _h_pi0, _h_eta, _h_etaP, _h_K0, _h_lam;
Histo1DPtr _h_2_pi0, _h_2_eta, _h_2_etaP, _h_2_K0, _h_2_lam;
Histo1DPtr _h_3_pi0[3], _h_3_eta[3], _h_3_etaP[3], _h_3_K0[3], _h_3_lam[3];
CounterPtr _w2, _w3;
/// @}

};

RIVET_DECLARE_PLUGIN(ALEPH_2000_I507531);

} ```