Rivet analyses

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

namespace Rivet {


  /// @brief Thrust like variable at Upsilon(1S,2S)
  class LENA_1981_I164397 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LENA_1981_I164397);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(), "UFS");
      declare(ChargedFinalState(), "FS");
      const double eval = sqrtS()/GeV;
      if (inRange(eval, 7.35, 7.49)) {
        _sqs = "7.35 - 7.49"s;
      }
      else if (inRange(eval, 8.629, 9.142)) {
        _sqs = "8.629 - 9.142"s;
      }
      else if (inRange(eval, 9.15, 9.41)) {
        _sqs = "9.15 - 9.41"s;
      }
      else if(isCompatibleWithSqrtS(9.5149, 1e-2)) {
        _sqs = "9.5149";
      }
      else if(isCompatibleWithSqrtS(9.9903, 1e-2)) {
        _sqs = "9.9903";
      }
      else if(isCompatibleWithSqrtS(9.4624, 1e-2)) {
        _sqs = "9.4624";
      }
      else if(isCompatibleWithSqrtS(10.0148, 1e-2)) {
        _sqs = "10.0148";
      }
      else raiseBeamErrorIf(true);
      book(_mult,         3, 1, 1);
      book(_h["9.5149"],  4, 1, 1);
      book(_h["9.9903"],  4, 1, 2);
      book(_h["9.4624"],  4, 1, 3);
      book(_h["10.0148"], 4, 1, 4);
      book(_c["9.5149"],  "TMP/weightSum_cont1");
      book(_c["9.9903"],  "TMP/weightSum_cont2");
      book(_c["9.4624"],  "TMP/weightSum_Ups1");
      book(_c["10.0148"], "TMP/weightSum_Ups2");
    }


    /// Recursively walk the decay tree to find the charged decay products of @a p
    void findDecayProducts(const Particle& mother, Particles& charged) const {
      for (const Particle& p: mother.children()) {
        const int id = p.pid();
        if (!p.children().empty())    findDecayProducts(p, charged);
        else if (PID::isCharged(id))  charged += p;
      }
    }

    // defn of thrust in paper used just the direction
    double thrustPrime(const LorentzTransform& boost, const Particles& particles) const {
      vector<Vector3> vecs; vecs.reserve(particles.size());
      for (const Particle& p : particles) {
        vecs.push_back(boost.transform(p.mom()).p3());
      }
      Thrust thrust;
      thrust.calc(vecs);
      return thrust.thrust();
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Find the Upsilons among the unstables
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      const Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
      if (upsilons.empty()) {
        MSG_DEBUG("No Upsilons found => continuum event");
        const Particles cfs = apply<ChargedFinalState>(event, "FS").particles();
        _mult->fill(_sqs,cfs.size());
        if (_sqs == "9.5149"s or _sqs == "9.9903"s) {
          _c[_sqs]->fill();
          LorentzTransform boost;
          _h[_sqs]->fill(thrustPrime(boost,cfs));
        }
      }
      // Upsilon(s) found
      else {
        for (const Particle& ups : upsilons) {
          const int parentId = ups.pid();
          Particles charged;
          // boost to rest frame (if required)
          LorentzTransform boost;
          if (ups.p3().mod() > 1*MeV)
            boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
          // Find the decay products we want
          findDecayProducts(ups, charged);
          if (parentId==553) {
            _c["9.4624"s]->fill();
            _mult->fill("9.4624"s, charged.size());
            _h["9.4624"s]->fill(thrustPrime(boost,charged));
          }
          else if (parentId == 100553) {
            _c["10.0148"s]->fill();
            _mult->fill("10.0148"s, charged.size());
            _h["10.0148"s]->fill(thrustPrime(boost,charged));
          }
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // charged particle multiplicity
      scale(_h, crossSectionPerEvent());
      scale(_c, crossSectionPerEvent());
      for (const auto& item : _c) {
        const double w = item.second->val();
        if (isZero(w))  continue;
        scale(_h[item.first], 1.0/w);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    BinnedProfilePtr<string> _mult;
    map<string,Histo1DPtr> _h;
    map<string,CounterPtr> _c;
    string _sqs = "";
    /// @}

  };


  RIVET_DECLARE_PLUGIN(LENA_1981_I164397);

}