Rivet analyses

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

namespace Rivet {


  /// @brief Charged-particle multiplicities, R, and event shapes between 12 and 43 GeV
  class TASSO_1984_I195333 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1984_I195333);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // Initialise and register projections
      const FinalState fs;
      declare(fs, "FS");
      declare(Sphericity(fs), "Sphericity");
      declare(Thrust(fs), "Thrust");

      // Counters for R
      book(_h_hadrons, "_sigma_hadrons", refData<YODA::BinnedEstimate<string>>(1,1,1));
      book(_h_muons,   "_sigma_muons"  , refData<YODA::BinnedEstimate<string>>(1,1,1));
      for (const string& en : _h_hadrons.binning().edges<0>()) {
        const size_t idx = en.find("-");
        if (idx!=std::string::npos) {
          const double emin = std::stod(en.substr(0,idx));
          const double emax = std::stod(en.substr(idx+1,string::npos));
          if (inRange(sqrtS()/GeV, emin, emax)) {
            _sqs[0] = en; break;
          }
        }
        else {
          const double eval = std::stod(en)*GeV;
          if (isCompatibleWithSqrtS(eval)) {
            _sqs[0] = en; break;
          }
        }
      }
      size_t ih = 1;
      for (double eVal : {14*GeV, 22*GeV, 34*GeV}) {
        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal)) _sqs[2] = en;
        book(_c[en], "_sumW_"+en);
        book(_d[en+"mult"], 3, 1, ih);
        book(_h[en+"p"],    5, 1, ih);
        book(_h[en+"xp"],   6, 1, ih);
        book(_h[en+"pl"],   7, 1, ih);
        book(_h[en+"pt"],   8, 1, ih);
        book(_h[en+"pt2"],  9, 1, ih);
        book(_h[en+"xl"],  10, 1, ih);
        book(_h[en+"xT"],  11, 1, ih);
        book(_h[en+"S"],   12, 1, ih);
        book(_h[en+"T"],   13, 1, ih);
        book(_h[en+"y"],   14, 1, ih);
        ++ih;
      }
      // average event shapes
      book(_p["charged"], 4, 1, 1);
      book(_p["total"], "TMP/ntotal",refData<YODA::BinnedEstimate<string>>(4,1,2));
      book(_p["sphericity"], 4, 1,  4);
      book(_p["thrust"],     4, 1,  5);
      book(_p["p_total"],    4, 1,  6);
      book(_p["p_l"],        4, 1,  7);
      book(_p["pt"],         4, 1,  8);
      book(_p["pt2"],        4, 1,  9);
      book(_p["pt2_in"],     4, 1, 10);
      book(_p["pt2_out"],    4, 1, 11);
      for (const string& en : _p["charged"].binning().edges<0>()) {
        const size_t idx = en.find("-");
        if (idx != string::npos) {
          const double emin = stod(en.substr(0,idx));
          const double emax = stod(en.substr(idx+1,string::npos));
          if (inRange(sqrtS()/GeV, emin, emax)) {
            _sqs[1] = en; break;
          }
        }
        else {
          const double eval = stod(en)*GeV;
          if (isCompatibleWithSqrtS(eval)) {
            _sqs[1] = en; break;
          }
        }
      }
      raiseBeamErrorIf(_sqs[0].empty() && _sqs[1].empty());
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const FinalState& fs = apply<FinalState>(event, "FS");

      map<long,int> nCount;
      double ntotal(0);
      unsigned int nCharged(0);
      for (const Particle& p : fs.particles()) {
        nCount[p.pid()] += 1;
        if ((p.pid()!=PID::PHOTON && p.abspid()!=PID::ELECTRON) ||
            p.parents().empty() || p.parents()[0].pid()!=PID::PI0) {
          ntotal += 1.;
        }
        else if (p.pid()==PID::PHOTON) {
          ntotal += 0.5 ;
        }
        else {
          ntotal += 0.25;
          nCharged -=1;
        }
        if (PID::isCharged(p.pid())) {
          ++nCharged;
        }
      }
      // mu+mu- + photons
      if (nCount[-13]==1 && nCount[13]==1 && int(fs.particles().size())==2+nCount[22]) {
        _h_muons->fill(_sqs[0]); vetoEvent;
      }
      // everything else
      _h_hadrons->fill(_sqs[0]);
      _p["charged"]->fill(_sqs[1], nCharged);
      _p["total"]->fill(_sqs[1], ntotal);
      // thrust
      const Thrust& thrust = apply<Thrust>(event, "Thrust");
      _p["thrust"]->fill(_sqs[1],thrust.thrust());
      // sphericity
      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
      _p["sphericity"]->fill(_sqs[1],sphericity.sphericity());
      // global distributions
      if (!_sqs[2].empty()) {
        _c[_sqs[2]]->fill();
        _d[_sqs[2]+"mult"]->fill(nCharged);
        _h[_sqs[2]+"S"]->fill(sphericity.sphericity());
        _h[_sqs[2]+"T"]->fill(thrust.thrust());
      }
      // single particle distributions
      for (const Particle& p : fs.particles()) {
        if (!PID::isCharged(p.pid())) continue;
        const Vector3 mom3 = p.p3();
        double pp = mom3.mod();
        const double mom = dot(sphericity.sphericityAxis(), mom3);
        const double pTin = dot(mom3, sphericity.sphericityMajorAxis());
        const double pTout = dot(mom3, sphericity.sphericityMinorAxis());
        double pT2 = sqr(pTin) + sqr(pTout);
        double pT = sqrt(pT2);
        _p["p_total"]->fill(_sqs[1],pp);
        _p["p_l"]->fill(_sqs[1],fabs(mom));
        _p["pt2_in"]->fill(_sqs[1],sqr(pTin));
        _p["pt2_out"]->fill(_sqs[1],sqr(pTout));
        _p["pt2"]->fill(_sqs[1],pT2);
        _p["pt"]->fill(_sqs[1],pT);
        if (!_sqs[2].empty()) {
          const double rap = 0.5 * log((p.E() + mom) / (p.E() - mom));
          _h[_sqs[2]+"p"]->fill(pp);
          _h[_sqs[2]+"xp"]->fill(2.*pp/sqrtS());
          _h[_sqs[2]+"pl"]->fill(fabs(mom));
          _h[_sqs[2]+"xl"]->fill(2.*fabs(mom)/sqrtS());
          _h[_sqs[2]+"pt2"]->fill(pT2);
          _h[_sqs[2]+"pt"]->fill(pT);
          _h[_sqs[2]+"xT"]->fill(2.*pT/sqrtS());
          _h[_sqs[2]+"y"]->fill(fabs(rap));
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      BinnedEstimatePtr<string> mult;
      book(mult, 1, 1, 1);
      divide(_h_hadrons, _h_muons, mult);
      // charged particle multiplicity distribution
      normalize(_d);
      // charged fraction
      book(mult, 4, 1, 2);
      divide(_p["charged"], _p["total"], mult);
      book(mult, 4, 1, 3);
      const auto& b = _p["charged"]->binAt(_sqs[1]);
      mult->binAt(_sqs[1]).set(b.stdDev(2),b.stdErr(2));
      // scale the distributions
      scale(_c, crossSectionPerEvent());
      scale(_h, crossSectionPerEvent());
      for (auto& item : _h) {
        const double w = _c[item.first.substr(0,2)]->sumW();
        if (!isZero(w))  scale(item.second, 1.0/w);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    map<string,CounterPtr> _c;
    map<string,Histo1DPtr> _h;
    map<string,BinnedHistoPtr<int>> _d;
    map<string,BinnedProfilePtr<string>> _p;
    BinnedHistoPtr<string> _h_hadrons, _h_muons;
    string _sqs[3];
    /// @}

  };

  RIVET_DECLARE_PLUGIN(TASSO_1984_I195333);

}