Rivet analyses

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

namespace Rivet {


  /// @brief ARGUS pi+/-, K+/- and proton/antiproton spectrum at Upsilon(4S)
  ///
  /// @author Peter Richardson
  class ARGUS_1993_I340894 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1993_I340894);


    void analyze(const Event& e) {
      // Find the upsilons
      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
        _weightSum->fill();
        Particles pionsA,pionsB,protonsA,protonsB,kaons;
        // Find the decay products we want
        findDecayProducts(p, pionsA, pionsB, protonsA, protonsB, kaons);
        LorentzTransform cms_boost;
        if (p.p3().mod() > 1*MeV)
          cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
        for (size_t ix = 0; ix < pionsA.size(); ++ix) {
          FourMomentum ptemp(pionsA[ix].momentum());
          FourMomentum p2 = cms_boost.transform(ptemp);
          double pcm = cms_boost.transform(ptemp).vector3().mod();
          _histPiA->fill(pcm);
        }
        _multPiA->fill(Ecm,double(pionsA.size()));
        for (size_t ix = 0; ix < pionsB.size(); ++ix) {
          double pcm = cms_boost.transform(pionsB[ix].momentum()).vector3().mod();
          _histPiB->fill(pcm);
        }
        _multPiB->fill(Ecm,double(pionsB.size()));
        for (size_t ix = 0; ix < protonsA.size(); ++ix) {
          double pcm = cms_boost.transform(protonsA[ix].momentum()).vector3().mod();
          _histpA->fill(pcm);
        }
        _multpA->fill(Ecm,double(protonsA.size()));
        for (size_t ix = 0; ix < protonsB.size(); ++ix) {
          double pcm = cms_boost.transform(protonsB[ix].momentum()).vector3().mod();
          _histpB->fill(pcm);
        }
        _multpB->fill(Ecm,double(protonsB.size()));
        for (size_t ix = 0 ;ix < kaons.size(); ++ix) {
          double pcm = cms_boost.transform(kaons[ix].momentum()).vector3().mod();
          _histKA->fill(pcm);
          _histKB->fill(pcm);
        }
        _multK->fill(Ecm,double(kaons.size()));
      }
    }


    void finalize() {
      if (_weightSum->val() > 0.) {
        scale(_histPiA, 1. / *_weightSum);
        scale(_histPiB, 1. / *_weightSum);
        scale(_histKA , 1. / *_weightSum);
        scale(_histKB , 1. / *_weightSum);
        scale(_histpA , 1. / *_weightSum);
        scale(_histpB , 1. / *_weightSum);
        scale(_multPiA, 1. / *_weightSum);
        scale(_multPiB, 1. / *_weightSum);
        scale(_multK  , 1. / *_weightSum);
        scale(_multpA , 1. / *_weightSum);
        scale(_multpB , 1. / *_weightSum);
      }
    }


    void init() {
      declare(UnstableParticles(), "UFS");

      // spectra
      book(_histPiA ,1, 1, 1);
      book(_histPiB ,2, 1, 1);
      book(_histKA  ,3, 1, 1);
      book(_histKB  ,6, 1, 1);
      book(_histpA  ,4, 1, 1);
      book(_histpB  ,5, 1, 1);
      // multiplicities
      book(_multPiA , 7, 1, 1);
      book(_multPiB , 8, 1, 1);
      book(_multK   , 9, 1, 1);
      book(_multpA  ,10, 1, 1);
      book(_multpB  ,11, 1, 1);

      book(_weightSum, "TMP/weightSum");
    }


  private:

    /// @{
    /// Count of weights
    CounterPtr _weightSum;
    /// Spectra
    Histo1DPtr _histPiA, _histPiB, _histKA, _histKB, _histpA, _histpB;
    /// Multiplicities
    BinnedHistoPtr<string> _multPiA, _multPiB, _multK, _multpA, _multpB;
    const string Ecm = "10.58";
    /// @}


    void findDecayProducts(Particle parent, Particles & pionsA, Particles & pionsB,
                           Particles & protonsA, Particles & protonsB, Particles & kaons) {
      int parentId = parent.pid();
      for (const Particle & p : parent.children()) {
        int id = abs(p.pid());
        if (id == PID::PIPLUS) {
          if (parentId != PID::LAMBDA && parentId != PID::K0S) {
            pionsA.push_back(p);
            pionsB.push_back(p);
          }
          else
            pionsB.push_back(p);
        }
        else if (id == PID::PROTON) {
          if (parentId != PID::LAMBDA && parentId != PID::K0S) {
            protonsA.push_back(p);
            protonsB.push_back(p);
          }
          else
            protonsB.push_back(p);
        }
        else if (id == PID::KPLUS) {
          kaons.push_back(p);
        }
        else if (!p.children().empty())
          findDecayProducts(p, pionsA, pionsB, protonsA, protonsB, kaons);
      }
    }


  };



  RIVET_DECLARE_ALIASED_PLUGIN(ARGUS_1993_I340894, ARGUS_1993_S2653028);

}