Rivet analyses

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

namespace Rivet {

  /// @brief Event shapes at 55.2
  class AMY_1990_I283337 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(AMY_1990_I283337);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(Beam(), "Beams");
      const FinalState fs;
      declare(fs, "FS");
      const Sphericity sphere(fs);
      declare(sphere, "Sphericity");
      const Thrust thrust(fs);
      declare(thrust, "Thrust");
      declare(Hemispheres(sphere), "Hemispheres");
      // histograms
      book(_histRapidityT , 1, 1, 1);
      book(_histMajor     ,13, 1, 1);

      book(_histScaledMom , 2, 1, 1);
      book(_histPl        , 3, 1, 1);
      book(_histPt        , 4, 1, 1);
      book(_histPt2       , 5, 1, 1);
      book(_histPtIn      , 6, 1, 1);
      book(_histPtOut     , 7, 1, 1);
      book(_histMeanPtIn2 , 8, 1, 1);
      book(_histMeanPtOut2, 9, 1, 1);
      book(_histNtheta    ,10, 1, 1);
      book(_histEtheta    ,11, 1, 1);
      book(_histThrust    ,12, 1, 1);
      book(_histMinor     ,14, 1, 1);
      book(_histOblateness,15, 1, 1);
      book(_histSphericity,16, 1, 1);
      book(_histAplanarity,17, 1, 1);
      book(_histQx        ,18, 1, 1);
      book(_histQ21       ,19, 1, 1);
      book(_histRhoLight  ,20, 1, 1);
      book(_histRhoHeavy  ,21, 1, 1);
      book(_histRhoDiff   ,22, 1, 1);
      book(_wSum,"TMP/wSum");
    }


    /// Perform the per-event analysis
    void analyze(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");
      _wSum->fill();

      // 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;
      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);

      // Thrusts
      MSG_DEBUG("Calculating thrust");
      const Thrust& thrust = apply<Thrust>(event, "Thrust");
      _histThrust    ->fill(thrust.thrust()     );
      _histMajor     ->fill(thrust.thrustMajor());
      _histMinor     ->fill(thrust.thrustMinor());
      _histOblateness->fill(thrust.oblateness() );
      // Sphericities
      MSG_DEBUG("Calculating sphericity");
      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
      _histSphericity->fill(sphericity.sphericity());
      _histAplanarity->fill(sphericity.aplanarity());
      _histQx        ->fill(sqrt(1./3.)*(sphericity.lambda1()-sphericity.lambda2()));
      _histQ21       ->fill(sphericity.lambda2()-sphericity.lambda3());
      // Hemispheres
      MSG_DEBUG("Calculating hemisphere variables");
      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
      _histRhoHeavy->fill(hemi.scaledM2high());
      _histRhoLight->fill(hemi.scaledM2low() );
      _histRhoDiff ->fill(hemi.scaledM2diff());
      // single particle distributions
      double pTIn2(0.),pTOut2(0.);
      unsigned int nCharged(0);
      for (const Particle& p : fs.particles()) {
        // Get momentum and energy of each particle.
        const Vector3 mom3 = p.p3();
        const double energy = p.E();
        const double mom = mom3.mod();
        const double scaledMom = mom/meanBeamMom;
        const double momT = dot(thrust.thrustAxis(), mom3);
        const double momS = dot(sphericity.sphericityAxis(), mom3);
        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
        const double pT = sqrt(pow(pTinS, 2) + pow(pToutS, 2));

        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
        double angle = sphericity.sphericityAxis().angle(p.p3())/M_PI*180.;
        if(angle>90.) angle=180.-angle;
        if(PID::isCharged(p.pid())) {
          _histScaledMom->fill(scaledMom);
          _histRapidityT->fill(fabs(rapidityT));
          _histPl       ->fill(fabs(momS)     );
          _histPt       ->fill(pT             );
          _histPt2      ->fill(sqr(pT)        );
          _histPtIn     ->fill(fabs(pTinS)    );
          _histPtOut    ->fill(fabs(pToutS)   );
          pTIn2  += sqr(pTinS);
          pTOut2 += sqr(pToutS);
          _histNtheta->fill(angle);
          ++nCharged;
        }
        _histEtheta->fill(angle,energy); 
      }
      if (nCharged) {
        _histMeanPtIn2 ->fill( pTIn2/nCharged);
        _histMeanPtOut2->fill(pTOut2/nCharged);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // histograms
      scale(_histRapidityT , 1./ *_wSum);
      scale(_histScaledMom , 1./ *_wSum);
      scale(_histPl        , 1./ *_wSum);
      scale(_histPt        , 1./ *_wSum);
      scale(_histPt2       , 1./ *_wSum);
      scale(_histPtIn      , 1./ *_wSum);
      scale(_histPtOut     , 1./ *_wSum);
      scale(_histMeanPtIn2 , 1./ *_wSum);
      scale(_histMeanPtOut2, 1./ *_wSum);
      scale(_histNtheta    , 1./ *_wSum);
      scale(_histEtheta    , 1./ *_wSum);
      scale(_histThrust    , 1./ *_wSum);
      scale(_histMajor     , 1./ *_wSum);
      scale(_histMinor     , 1./ *_wSum);
      scale(_histOblateness, 1./ *_wSum);
      scale(_histSphericity, 1./ *_wSum);
      scale(_histAplanarity, 1./ *_wSum);
      scale(_histQx        , 1./ *_wSum);
      scale(_histQ21       , 1./ *_wSum);
      scale(_histRhoLight  , 1./ *_wSum);
      scale(_histRhoHeavy  , 1./ *_wSum);
      scale(_histRhoDiff   , 1./ *_wSum);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _histRapidityT, _histScaledMom, _histPl, _histPt, _histPt2, _histPtIn, _histPtOut,
      _histMeanPtIn2, _histMeanPtOut2, _histNtheta, _histEtheta, _histThrust, _histMajor, _histMinor,
      _histOblateness, _histSphericity, _histAplanarity, _histQx, _histQ21, _histRhoLight,
      _histRhoHeavy, _histRhoDiff;
    CounterPtr _wSum;
    /// @}


  };


  RIVET_DECLARE_PLUGIN(AMY_1990_I283337);


}