Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {


  /// @brief K+- Lambda asymmetries
  class DELPHI_1995_I382285 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1995_I382285);


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

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

      // Initialise and register projections
      declare(Beam(), "Beams");
      declare(FinalState(), "FS");
      declare(UnstableParticles(), "UFS");

      // Book histograms
      book(_h_Kp, "/TMP/cos_Kp", 20, -1., 1.);
      book(_h_Km, "/TMP/cos_Km", 20, -1., 1.);
      book(_h_lm, "/TMP/cos_lm", 20, -1., 1.);
      book(_h_lb, "/TMP/cos_lb", 20, -1., 1.);
    }


    /// 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;
      Vector3 beamAxis;
      if (beams.first.pid() == 11) {
        beamAxis = beams.first.momentum().p3().unit();
      }
      else {
        beamAxis = beams.second.momentum().p3().unit();
      }
      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      for (const Particle& p : ufs.particles(Cuts::abspid == 3122 || Cuts::abspid == 321)) {
        double modp = p.momentum().p3().mod();
        if (p.abspid() == 321) {
          if (modp < 10. || modp > 18.) continue;
          const double cK = beamAxis.dot(p.momentum().p3().unit());
          if (p.pid() > 0)
            _h_Kp->fill(cK);
          else
            _h_Km->fill(cK);
        }
        else {
          if (modp < 11.41 || modp > 22.82) continue;
          const double cLam = beamAxis.dot(p.momentum().p3().unit());
          if (p.pid() > 0)
            _h_lm->fill(cLam);
          else
            _h_lb->fill(cLam);
        }
      }
    }

    pair<double, double> calcAsymmetry(const Estimate1DPtr& hist) {
      double sum1(0.), sum2(0.);
      for (const auto& bin : hist->bins()) {
        double Oi = bin.val();
        if (Oi == 0.) continue;
        const double bi = 4. * (bin.xMax() + bin.xMin())
            / (3. + sqr(bin.xMax()) + bin.xMax() * bin.xMin() + sqr(bin.xMin()));
        const double Ei = bin.errAvg();
        sum1 += sqr(bi / Ei);
        sum2 += bi / sqr(Ei) * Oi;
      }
      return make_pair(sum2 / sum1, sqrt(1. / sum1));
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h_Kp);
      normalize(_h_Km);
      Estimate1DPtr sK;
      book(sK, "a_K", 20, -1., 1.);
      asymm(_h_Kp, _h_Km, sK);
      pair<double, double> alpha = calcAsymmetry(sK);
      BinnedEstimatePtr<string> h_K;
      book(h_K, 1, 1, 1);
      h_K->bin(1).set(-alpha.first, alpha.second);

      normalize(_h_lm);
      normalize(_h_lb);
      Estimate1DPtr sLam;
      book(sLam, "a_Lam", 20, -1., 1.);
      asymm(_h_lm, _h_lb, sLam);
      alpha = calcAsymmetry(sLam);
      BinnedEstimatePtr<string> h_lam;
      book(h_lam, 1, 1, 2);
      h_lam->bin(1).set(alpha.first, alpha.second);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_Kp, _h_Km, _h_lm, _h_lb;
    const string Ecm = "91.2";
    /// @}
  };


  RIVET_DECLARE_PLUGIN(DELPHI_1995_I382285);


}