Rivet analyses

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

namespace Rivet {


  /// @brief Azimuthal ordering of charged hadrons
  class ATLAS_2012_I1091481 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1091481);

    /// Book histograms and initialise projections before the run
    void init() {

      ChargedFinalState cfs100(Cuts::abseta < 2.5 && Cuts::pT > 0.1*GeV);
      declare(cfs100,"CFS100");
      ChargedFinalState cfs500(Cuts::abseta < 2.5 && Cuts::pT > 0.5*GeV);
      declare(cfs500,"CFS500");

      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
        size_t ih = bool(en == "900") + 1;

        book(_h[en+"E_10_100"], ih, 1, 1);
        book(_h[en+"E_1_100"],  ih, 1, 2);
        book(_h[en+"E_10_500"], ih, 1, 3);

        book(_h[en+"eta_10_100"], ih, 2, 1);
        book(_h[en+"eta_1_100"],  ih, 2, 2);
        book(_h[en+"eta_10_500"], ih, 2, 3);

        book(_c[en+"inclusive"], "_norm_inclusive_"+en);
        book(_c[en+"lowPt"],     "_norm_lowPt_"+en);
        book(_c[en+"pt500"],     "_norm_pt500_"+en);
      }
      raiseBeamErrorIf(_sqs.empty());
    }


    // Recalculate particle energy assuming pion mass
    double getPionEnergy(const Particle& p) {
      double m_pi = 0.1396*GeV;
      double p2 = p.p3().mod2()/(GeV*GeV);
      return sqrt(sqr(m_pi) + p2);
    }


    // S_eta core for one event
    //
    //  -1 + 1/Nch * |sum_j^Nch exp[i*(xi eta_j - Phi_j)]|^2
    //
    double getSeta(const Particles& part, double xi) {
      std::complex<double> c_eta (0.0, 0.0);
      for (const Particle& p : part) {
        double eta = p.eta();
        double phi = p.phi();
        double arg = xi*eta-phi;
         std::complex<double> temp(cos(arg), sin(arg));
         c_eta += temp;
      }
      return std::norm(c_eta)/part.size() - 1.0;
    }


    // S_E core for one event
    //
    //  -1 + 1/Nch * |sum_j^Nch exp[i*(omega X_j - Phi_j)]|^2
    //
    double getSE(const Particles& part, double omega) {
      double Xj = 0.0;
      std::complex<double> c_E (0.0, 0.0);
      for (unsigned int i=0; i < part.size(); ++i) {
        Xj += 0.5*getPionEnergy(part[i]);
        double phi = part[i].phi();
        double arg = omega*Xj - phi;
        std::complex<double> temp(cos(arg), sin(arg));
        c_E += temp;
        Xj += 0.5*getPionEnergy(part[i]);
      }
      return std::norm(c_E)/part.size() - 1.0;
    }


    // Convenient fill function
    void fillS(Histo1DPtr h, const Particles& part, bool SE=true) {
      // Loop over bins, take bin centers as parameter values
      // @todo use barchart()
      for (auto& b : h->bins()) {
        double x = b.xMid();
        double width = b.xWidth();
        double y;
        if (SE)  y = getSE(part,   x);
        else     y = getSeta(part, x);
        h->fill(x, y * width);
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Charged fs
      const ChargedFinalState& cfs100  = apply<ChargedFinalState>(event, "CFS100");
      const Particles          part100 = cfs100.particles(cmpMomByEta);
      const ChargedFinalState& cfs500  = apply<ChargedFinalState>(event, "CFS500");
      const Particles&         part500 = cfs500.particles(cmpMomByEta);

      // Veto event if the most inclusive phase space has less than 10 particles and the max pT is > 10 GeV
      if (part100.size() < 11) vetoEvent;
      double ptmax = cfs100.particlesByPt()[0].pT()/GeV;
      if (ptmax > 10.0) vetoEvent;

      // Fill the pt>100, pTmax<10 GeV histos
      fillS(_h[_sqs+"E_10_100"],   part100, true);
      fillS(_h[_sqs+"eta_10_100"], part100, false);
      _c[_sqs+"inclusive"]->fill();

      // Fill the pt>100, pTmax<1 GeV histos
      if (ptmax < 1.0) {
        fillS(_h[_sqs+"E_1_100"],   part100, true);
        fillS(_h[_sqs+"eta_1_100"], part100, false);
        _c[_sqs+"lowPt"]->fill();
      }

      // Fill the pt>500, pTmax<10 GeV histos
      if (part500.size() > 10) {
        fillS(_h[_sqs+"E_10_500"],   part500, true );
        fillS(_h[_sqs+"eta_10_500"], part500, false);
        _c[_sqs+"pt500"]->fill();
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      // The scaling takes the multiple fills per event into account
      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));

        scale(_h[en+"E_10_100"], 1.0/ *_c[en+"inclusive"]);
        scale(_h[en+ "E_1_100"], 1.0/ *_c[en+"lowPt"]);
        scale(_h[en+"E_10_500"], 1.0/ *_c[en+"pt500"]);

        scale(_h[en+"eta_10_100"], 1.0/ *_c[en+"inclusive"]);
        scale(_h[en+ "eta_1_100"], 1.0/ *_c[en+"lowPt"]);
        scale(_h[en+"eta_10_500"], 1.0/ *_c[en+"pt500"]);
      }
    }


  private:

    map<string,Histo1DPtr> _h;
    map<string,CounterPtr> _c;

    string _sqs = "";
  };


  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1091481);

}