Rivet analyses

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

namespace Rivet {

  /// Measurement of event shapes in minimum-bias events at 13 TeV
  class CMS_2025_I2924533 : public Analysis {
  public:
    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2025_I2924533);

    void init() {
      // Initialise and register projections

      // The basic final-state projection:
      // all final-state particles within
      // the given eta acceptance
      const ChargedFinalState cfs05(Cuts::abseta < 2.4 && Cuts::pT > 0.5 * GeV);
      declare(cfs05, "CFS05");
      // Book histograms
      // specify custom binning
      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
      book(_h["nch"], 1, 1, 1);
      book(_h["mass"], 2, 1, 1);
      book(_h["sphericity"], 3, 1, 1);
      book(_h["thrust"], 4, 1, 1);
      book(_h["broadening"], 5, 1, 1);
      book(_h["transverseThrust"], 7, 1, 1);
      book(_h["transverseSpherocity"], 8, 1, 1);
      book(_h2d["nch_mass"], 52, 1, 1);
      book(_h2d["nch_sphericity"], 53, 1, 1);
      book(_h2d["nch_thrust"], 54, 1, 1);
      book(_h2d["nch_broadening"], 55, 1, 1);
      book(_h2d["nch_transverseThrust"], 57, 1, 1);
      book(_h2d["nch_transverseSpherocity"], 58, 1, 1);
      vector<double> nch_bins = {3.0, 10.0, 20.0, 30.0, 40.0, 140.0};
      book(_b["nch_mass"], nch_bins);
      book(_b["nch_sphericity"], nch_bins);
      book(_b["nch_thrust"], nch_bins);
      book(_b["nch_broadening"], nch_bins);
      book(_b["nch_transverseThrust"], nch_bins);
      book(_b["nch_transverseSpherocity"], nch_bins);
      for (size_t i = 0; i < _b["nch_mass"]->numBins(); ++i) {
        book(_b["nch_mass"]->bin(i + 1), 17 + i, 1, 1);
        book(_b["nch_sphericity"]->bin(i + 1), 22 + i, 1, 1);
        book(_b["nch_thrust"]->bin(i + 1), 27 + i, 1, 1);
        book(_b["nch_broadening"]->bin(i + 1), 32 + i, 1, 1);
        book(_b["nch_transverseThrust"]->bin(i + 1), 42 + i, 1, 1);
        book(_b["nch_transverseSpherocity"]->bin(i + 1), 47 + i, 1, 1);
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const Particles& chargedParticles = apply<ChargedFinalState>(event, "CFS05").particlesByPt();
      vector<FourMomentum> momenta;
      double num = chargedParticles.size();
      double massSum = 0;
      FourMomentum fourMom(0, 0, 0, 0);
      for (const Particle& p : chargedParticles) {
        fourMom = fourMom + p.momentum();
        momenta.push_back(p.momentum());
      }
      massSum = fourMom.mass() / GeV;
      // Veto event if there are less than 3 charged particles within the abseta < 2.4 and pt > 0.5 GeV
      if (num < 3)
        vetoEvent;
      double t, s, ts, tt, b;
      // Boost the charged particle momenta to the center-of-mass frame for calculation of sphericity, thrust and broadening
      // The transverse spherocity and transverse thrust are calculated in the lab frame.
      LorentzTransform boostTrafo;
      boostTrafo.setBetaVec(-fourMom.betaVec());
      vector<FourMomentum> momentaBoost;
      vector<Vector3> momenta3Boost;
      vector<Vector3> momenta3Unboost;
      vector<Vector3> momentaTransUnboost;
      for (size_t ip = 0; ip < momenta.size(); ++ip) {
        momentaBoost.push_back(boostTrafo.transform(momenta[ip]));
        momenta3Boost.push_back(boostTrafo.transform(momenta[ip]).p3());
        momentaTransUnboost.push_back(momenta[ip].p3());
        momentaTransUnboost[ip].setZ(0);
        momenta3Unboost.push_back(momenta[ip].p3());
      }
      // Calculate the event shape variables
      Thrust thrust;
      thrust.calc(momenta3Boost);
      t = 1.0 - thrust.thrust();
      Thrust transverseThrust;
      transverseThrust.calc(momentaTransUnboost);
      tt = 1.0 - transverseThrust.thrust();
      Sphericity sphericity;
      sphericity.calc(momenta3Boost);
      s = sphericity.sphericity();
      Spherocity transverseSpherocity;
      transverseSpherocity.calc(momenta3Unboost);
      ts = transverseSpherocity.spherocity();
      Hemispheres hemiSpheres(thrust);
      hemiSpheres.calc(thrust.thrustAxis(), momentaBoost);
      b = hemiSpheres.Bsum();

      // Fill the 1-dimensional histograms
      _h["nch"]->fill(num);
      _h["mass"]->fill(massSum);
      _h["sphericity"]->fill(s);
      _h["thrust"]->fill(t);
      _h["broadening"]->fill(b);
      _h["transverseSpherocity"]->fill(ts);
      _h["transverseThrust"]->fill(tt);

      // Fill the 2-dimensional histograms
      _h2d["nch_mass"]->fill(num, massSum);
      _h2d["nch_sphericity"]->fill(num, s);
      _h2d["nch_thrust"]->fill(num, t);
      _h2d["nch_broadening"]->fill(num, b);
      _h2d["nch_transverseSpherocity"]->fill(num, ts);
      _h2d["nch_transverseThrust"]->fill(num, tt);

      // Fill the histograms in slices of charged particle multiplicity
      _b["nch_mass"]->fill(num, massSum);
      _b["nch_sphericity"]->fill(num, s);
      _b["nch_thrust"]->fill(num, t);
      _b["nch_broadening"]->fill(num, b);
      _b["nch_transverseSpherocity"]->fill(num, ts);
      _b["nch_transverseThrust"]->fill(num, tt);
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h);
      normalize(_h2d);
      normalizeGroup(_b);
    }

  private:
    /// @name Histograms
    map<string, Histo1DPtr> _h;
    map<string, Histo2DPtr> _h2d;
    map<string, Histo1DGroupPtr> _b;
  };

  RIVET_DECLARE_PLUGIN(CMS_2025_I2924533);

}  // namespace Rivet