Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/Hemispheres.hh"
#include "Rivet/Projections/Sphericity.hh"
#include "Rivet/Projections/Thrust.hh"

namespace Rivet {


  /// @brief Event shapes at 29 GeV
  class MARKII_1988_I246184 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(MARKII_1988_I246184);


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

    /// Book histograms and initialise projections before the run
    void init() {
      const FinalState fs;
      declare(fs, "FS");
      const ChargedFinalState cfs;
      declare(cfs, "CFS");
      Sphericity sphere(fs);
      declare(sphere, "Sphericity");
      declare(Thrust(fs), "Thrust");
      declare(Hemispheres(sphere), "Hemispheres");
      // histograms
      unsigned int ioff = 18;
      for (unsigned int ix = 0; ix < 3; ++ix) {
        book(_histAplanarity[ix], 1 + ioff * ix, 1, 1);
        book(_histQx[ix], 2 + ioff * ix, 1, 1);
        book(_histQ2Q1[ix], 3 + ioff * ix, 1, 1);
        book(_histSphericity[ix], 4 + ioff * ix, 1, 1);
        book(_histThrust[ix], 5 + ioff * ix, 1, 1);
        book(_histMinor[ix], 6 + ioff * ix, 1, 1);
        book(_histOblateness[ix], 7 + ioff * ix, 1, 1);
        book(_histMJetBroad[ix], 8 + ioff * ix, 1, 1);
        book(_histMJetSlim[ix], 9 + ioff * ix, 1, 1);
        book(_histMJetDiff[ix], 10 + ioff * ix, 1, 1);
        book(_histScaledMom[ix], 15 + ioff * ix, 1, 1);
        book(_histPt2S[ix], 11 + ioff * ix, 1, 1);
        book(_histPtS[ix], 12 + ioff * ix, 1, 1);
        book(_histPtSIn[ix], 14 + ioff * ix, 1, 1);
        book(_histPtSOut[ix], 13 + ioff * ix, 1, 1);
        book(_histRapidityS[ix], 16 + ioff * ix, 1, 1);
        book(_histTheta[ix], 17 + ioff * ix, 1, 1);
        book(_histETheta[ix], 18 + ioff * ix, 1, 1);
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Sphericity related
      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
      for (unsigned int ix = 0; ix < 3; ++ix) {
        _histSphericity[ix]->fill(sphericity.sphericity());
        _histAplanarity[ix]->fill(sphericity.aplanarity());
        _histQx[ix]->fill((sphericity.lambda1() - sphericity.lambda2()) / sqrt(3.));
        _histQ2Q1[ix]->fill(sphericity.lambda2() - sphericity.lambda3());
      }
      // thrust related
      const Thrust& thrust = apply<Thrust>(event, "Thrust");
      for (unsigned int ix = 0; ix < 3; ++ix) {
        _histThrust[ix]->fill(thrust.thrust());
        _histMinor[ix]->fill(thrust.thrustMinor());
        _histOblateness[ix]->fill(thrust.oblateness());
      }
      // hemisphere related
      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
      double mWide = hemi.scaledM2high(), mNarrow = hemi.scaledM2low();
      if (!hemi.massMatchesBroadening()) swap(mWide, mNarrow);
      for (unsigned int ix = 0; ix < 3; ++ix) {
        _histMJetBroad[ix]->fill(mWide);
        _histMJetSlim[ix]->fill(mNarrow);
        _histMJetDiff[ix]->fill(hemi.scaledM2diff());
      }
      // dists w.r.t sphericity axis
      const FinalState& cfs = apply<FinalState>(event, "CFS");
      for (const Particle& p : cfs.particles()) {
        // Get momentum and energy of each particle.
        const Vector3 mom3 = p.p3();
        const double energy = p.E();
        // Scaled momenta.
        const double mom = mom3.mod();
        const double scaledMom = 2. * mom / sqrtS();
        // Get momenta components w.r.t. thrust and sphericity.
        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
        double pT2 = sqr(pTinS) + sqr(pToutS);
        double pT = sqrt(pT2);
        const double momS = dot(sphericity.sphericityAxis(), mom3);
        const double rapidityS = 0.5 * std::log((energy + momS) / (energy - momS));
        // angle
        double theta = sphericity.sphericityAxis().angle(mom3) / M_PI * 180.;
        if (theta > 90.) theta = 180. - theta;
        // fill histos
        for (unsigned int ix = 0; ix < 3; ++ix) {
          _histScaledMom[ix]->fill(scaledMom);
          _histPt2S[ix]->fill(fabs(pT2 / GeV));
          _histPtS[ix]->fill(fabs(pT / GeV));
          _histPtSIn[ix]->fill(fabs(pTinS / GeV));
          _histPtSOut[ix]->fill(fabs(pToutS / GeV));
          _histRapidityS[ix]->fill(fabs(rapidityS));
          _histTheta[ix]->fill(theta);
        }
      }
      // energy flow includes neutral w.r.t sphericity axis
      const FinalState& fs = apply<FinalState>(event, "FS");
      for (const Particle& p : fs.particles()) {
        // Get momentum and energy of each particle.
        const Vector3 mom3 = p.p3();
        const double energy = p.E();
        // angle
        double theta = sphericity.sphericityAxis().angle(mom3) / M_PI * 180.;
        if (theta > 90.) theta = 180. - theta;
        // fill histos
        for (unsigned int ix = 0; ix < 3; ++ix) {
          _histETheta[ix]->fill(theta, energy);
        }
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int ix = 0; ix < 3; ++ix) {
        scale(_histAplanarity[ix], 1. / sumOfWeights());
        scale(_histQx[ix], 1. / sumOfWeights());
        scale(_histQ2Q1[ix], 1. / sumOfWeights());
        scale(_histSphericity[ix], 1. / sumOfWeights());
        scale(_histThrust[ix], 1. / sumOfWeights());
        scale(_histMinor[ix], 1. / sumOfWeights());
        scale(_histOblateness[ix], 1. / sumOfWeights());
        scale(_histMJetBroad[ix], 1. / sumOfWeights());
        scale(_histMJetSlim[ix], 1. / sumOfWeights());
        scale(_histMJetDiff[ix], 1. / sumOfWeights());
        scale(_histScaledMom[ix], 1. / sumOfWeights());
        scale(_histPt2S[ix], 1. / sumOfWeights());
        scale(_histPtS[ix], 1. / sumOfWeights());
        scale(_histPtSIn[ix], 1. / sumOfWeights());
        scale(_histPtSOut[ix], 1. / sumOfWeights());
        scale(_histRapidityS[ix], 1. / sumOfWeights());
        scale(_histTheta[ix], 1. / sumOfWeights());
        scale(_histETheta[ix], 1. / sumOfWeights());
      }
    }

    ///@}


    /// @name Histograms
    ///@{
    Histo1DPtr _histAplanarity[3], _histQx[3], _histQ2Q1[3], _histSphericity[3];
    Histo1DPtr _histThrust[3], _histMinor[3], _histOblateness[3];
    Histo1DPtr _histMJetBroad[3], _histMJetSlim[3], _histMJetDiff[3];
    Histo1DPtr _histScaledMom[3], _histPt2S[3], _histPtS[3], _histPtSIn[3], _histPtSOut[3], _histRapidityS[3];
    Histo1DPtr _histTheta[3], _histETheta[3];
    ///@}
  };


  RIVET_DECLARE_PLUGIN(MARKII_1988_I246184);

}