Rivet analyses

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

namespace Rivet {


  /// @brief chi_b production at 7,8 TeV
  class LHCB_2014_I1308738 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2014_I1308738);


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

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

      declare(UnstableParticles(), "UFS");

      for (double eVal : allowedEnergies()) {
        const int en = int(eVal+0.5);
        if (isCompatibleWithSqrtS(eVal))  _sqs = en;

        int ih = int(en==8000)+1;

        // decays to 1S
        for (size_t ix=0; ix<3; ++ix) {
          book(_h_pT_1S[en+ix], "TMP/Ups_1_"+toString(en+ix), refData(3*ih-2,1,1+ix));
        }
        for (size_t ix=0; ix<2; ++ix) {
          book(_h_pT_1S[en+ix+3], "TMP/Ups_1_"+toString(en+ix+3), refData(3*ih-2,1,2+ix));
        }
        for (size_t ix=0; ix<2; ++ix) { // decays to 2S
          for (size_t iy=0; iy<2; ++iy) {
            book(_h_pT_2S[en+ix+2*iy], "TMP/Ups_2_"+toString(en+ix+2*iy), refData(3*ih-1,1,2+ix));
          }
        }
        for (size_t ix=0; ix<2; ++ix) { // decays to 3S
          book(_h_pT_3S[en+ix], "TMP/Ups_3_"+toString(en+ix), refData(3*ih,1,3));
        }
      }
      raiseBeamErrorIf(_sqs==0);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      // Final state of unstable particles to get particle spectra
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");

      // first the upsilon states for the denominator
      for (const Particle& p : ufs.particles(Cuts::pid==553 || Cuts::pid==100553 || Cuts::pid==200553)) {
        // rapidity cut
        double absrap=p.absrap();
        if (absrap<2. || absrap>4.5) continue;
        double xp = p.perp();
        if (p.pid()==553) {
          for (size_t ix=0; ix<2; ++ix) {
            _h_pT_1S[_sqs+3+ix]->fill(xp);
          }
        }
        else if (p.pid()==100553) {
          for (size_t ix=0; ix<2; ++ix) {
            _h_pT_2S[_sqs+2+ix]->fill(xp);
          }
        }
        else if (p.pid()==200553) {
          _h_pT_3S[_sqs+1]->fill(xp);
        }
      }
      // P states
      for (const Particle& p : ufs.particles(Cuts::pid==10551 || Cuts::pid==110551 || Cuts::pid==210551 ||
                                             Cuts::pid==20553 || Cuts::pid==120553 || Cuts::pid==220553 ||
                                             Cuts::pid==  555 || Cuts::pid==100555 || Cuts::pid==200555 )) {
        Particle Upsilon;
        if (p.children()[0].pid()==22 && p.children()[1].abspid()%100000==553) {
          Upsilon=p.children()[1];
        }
        else if (p.children()[1].pid()==22 && p.children()[0].abspid()%100000==553) {
          Upsilon=p.children()[0];
        }
        else {
          continue;
        }
        double absrap = Upsilon.absrap();
        if (absrap<2. || absrap>4.5) continue;
        unsigned int iups = p.children()[0].pid()/100000;
        unsigned int ichi = p.pid()/100000;
        double xp  = Upsilon.perp();
        if (iups==0)                 _h_pT_1S[_sqs+ichi]->fill(xp);
        else if (iups==1 && ichi>0)  _h_pT_2S[_sqs+ichi-1]->fill(xp);
        else if (iups==2 && ichi==2) _h_pT_3S[_sqs+0]->fill(xp);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (double eVal : allowedEnergies()) {
        const int en = int(eVal+0.5);
        int ih = int(en==8000)+1;
        // calculate the ratios
        // for 1S
        for (size_t iy=0; iy<3; ++iy) {
          Estimate1DPtr tmp;
          book(tmp,3*ih-2,1,iy+1);
          if (iy<2)  efficiency(_h_pT_1S[en+iy], _h_pT_1S[en+3], tmp);
          else       efficiency(_h_pT_1S[en+iy], _h_pT_1S[en+4], tmp);
          tmp->scale(100.);
        }
        Estimate1DPtr tmp;
        book(tmp, 3*ih-1, 1, 2);
        efficiency(_h_pT_2S[en+0], _h_pT_2S[en+2], tmp);
        tmp->scale(100.);
        book(tmp, 3*ih-1, 1, 3);
        efficiency(_h_pT_2S[en+1], _h_pT_2S[en+3], tmp);
        tmp->scale(100.);
        book(tmp, 3*ih, 1, 3);
        efficiency(_h_pT_3S[en+0], _h_pT_3S[en+1], tmp);
        tmp->scale(100.);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    map<size_t,Histo1DPtr> _h_pT_1S, _h_pT_2S, _h_pT_3S;
    int _sqs = 0;
    /// @}


  };


  RIVET_DECLARE_PLUGIN(LHCB_2014_I1308738);

}