Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"

namespace Rivet {


  /// UA1 minbias track multiplicities, \f$ p_\perp \f$ and \f$ E_\perp \f$
  class UA1_1990_I280412 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(UA1_1990_I280412);


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

    /// Book projections and histograms
    void init() {
      declare(ChargedFinalState(Cuts::abseta < 5.5), "TriggerFS");
      declare(ChargedFinalState(Cuts::abseta < 2.5), "TrackFS");
      const FinalState trkcalofs(Cuts::abseta < 2.5);
      declare(MissingMomentum(trkcalofs), "MET25");
      const FinalState calofs(Cuts::abseta < 6.0);
      declare(MissingMomentum(calofs), "MET60");

      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
        if (en == "63"s) book(_p[en+"Pt"], 8, 1, 1);
        else if (en == "200"s) {
          book(_h[en+"Nch"], 1, 1, 1);
          book(_h[en+"Esigd3p"], 2, 1, 1);
          book(_p[en+"Pt"], 6, 1, 1);
          book(_h[en+"Et"], 9, 1, 1);
          book(_p[en+"Etavg"], 12, 1, 1);
        }
        else if (en == "500"s) {
          book(_h[en+"Nch"], 1, 1, 2);
          book(_h[en+"Esigd3p"], 2, 1, 2);
          book(_h[en+"Et"], 10, 1, 1);
          book(_p[en+"Etavg"], 12, 1, 2);
        }
        else if (en == "900"s) {
          book(_h[en+"Nch"], 1, 1, 3);
          book(_h[en+"Esigd3p"], 2, 1, 3);
          book(_p[en+"Pt"], 7, 1, 1);
          book(_h[en+"Et"], 11, 1, 1);
          book(_p[en+"Etavg"], 12, 1, 3);
          book(_h[en+"Esigd3p08"], 3, 1, 1);
          book(_h[en+"Esigd3p40"], 4, 1, 1);
          book(_h[en+"Esigd3p80"], 5, 1, 1);
        }
        book(_c[en], "TMP/sumwTrig"+en);
      }
      raiseBeamErrorIf(_sqs.empty());
    }


    void analyze(const Event& event) {
      // Trigger
      const FinalState& trigfs = apply<FinalState>(event, "TriggerFS");
      unsigned int n_minus(0), n_plus(0);
      for (const Particle& p : trigfs.particles()) {
        const double eta = p.eta();
        if (inRange(eta, -5.5, -1.5)) n_minus++;
        else if (inRange(eta, 1.5, 5.5)) n_plus++;
      }
      if (n_plus == 0 || n_minus == 0) vetoEvent;
      _c[_sqs]->fill();

      // Use good central detector tracks
      const FinalState& cfs = apply<FinalState>(event, "TrackFS");
      const double Et25 = apply<MissingMomentum>(event, "MET25").scalarEt();
      const double Et60 = apply<MissingMomentum>(event, "MET60").scalarEt();
      const unsigned int nch = cfs.size();

      // Event level histos
      if (_sqs != "63"s) {
        _h[_sqs+"Nch"]->fill(nch);
        _h[_sqs+"Et"]->fill(Et60/GeV);
        _p[_sqs+"Etavg"]->fill(nch, Et25/GeV);
      }

      // Particle/track level histos
      const double deta = 2 * 5.0;
      const double dphi = TWOPI;
      const double dnch_deta = nch/deta;
      for (const Particle& p : cfs.particles()) {
        const double pt = p.pT();
        const double scaled_weight = 1.0/(deta*dphi*pt/GeV);
        if (_sqs != "500"s)  _p[_sqs+"Pt"]->fill(nch, pt/GeV);
        if (_sqs != "63"s)   _h[_sqs+"Esigd3p"]->fill(pt/GeV, scaled_weight);
        // Also fill for specific dn/deta ranges at 900 GeV
        if (_sqs == "900"s) {
          if (inRange(dnch_deta, 0.8, 4.0)) {
            _h[_sqs+"Esigd3p08"]->fill(pt/GeV, scaled_weight);
          } else if (inRange(dnch_deta, 4.0, 8.0)) {
            _h[_sqs+"Esigd3p40"]->fill(pt/GeV, scaled_weight);
          } else if (dnch_deta > 8.0) {
            _h[_sqs+"Esigd3p80"]->fill(pt/GeV, scaled_weight);
          }
        }
      }
    }


    void finalize() {
      for (const auto& item : _c) {
        if (isZero(item.second->val())) {
          MSG_DEBUG("No events passed the trigger!");
          continue;
        }
        const double sf = crossSectionPerEvent()/millibarn;
        if (item.first != "63"s) {
          scale(_h[item.first+"Nch"], 2*sf); ///< Factor of 2 for Nch bin widths?
          scale(_h[item.first+"Esigd3p"], sf);
          scale(_h[item.first+"Et"], sf);
        }
        if (item.first == "900"s) {
          // NB. Ref data is normalised to a fixed value not reproducible from MC.
          double scale08 = 0., scale40 = 0., scale80 = 0.;
          if (_h[item.first+"Esigd3p08"]->bin(1).sumW()) {
            scale08 = 0.933e5*_h[item.first+"Esigd3p08"]->bin(1).dVol()/_h[item.first+"Esigd3p08"]->bin(1).sumW();
          }
          scale(_h[item.first+"Esigd3p08"], scale08);
          if (_h[item.first+"Esigd3p40"]->bin(1).sumW()) {
            scale40 = 1.369e5*_h[item.first+"Esigd3p40"]->bin(1).dVol()/_h[item.first+"Esigd3p40"]->bin(1).sumW();
          }
          scale(_h[item.first+"Esigd3p40"], scale40);
          if (_h[item.first+"Esigd3p80"]->bin(1).sumW()) {
            scale80 = 1.657e5*_h[item.first+"Esigd3p80"]->bin(1).dVol()/_h[item.first+"Esigd3p80"]->bin(1).sumW();
          }
          scale(_h[item.first+"Esigd3p80"], scale80);
        }
      }
    }

    /// @}


  private:

    /// @name Weight counters
    /// @{
    map<string,CounterPtr> _c;
    /// @}

    /// @name Histogram collections
    /// @{
    map<string,Histo1DPtr> _h;
    map<string,Profile1DPtr> _p;
    string _sqs = "";
    /// @}

  };

  RIVET_DECLARE_ALIASED_PLUGIN(UA1_1990_I280412, UA1_1990_S2044935);
}