Rivet analyses

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

namespace Rivet {


  /// @brief Upsilon(1,2,3S) at 7 TeV
  class CMS_2013_I1225274 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1225274);

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

    /// Book histograms and initialise projections before the run
    void init() {
      // projection
      declare(UnstableParticles(), "UFS");
      // histograms
      for (unsigned int ix=0; ix<3; ++ix) {
        if (ix<2)  book(_h_total[ix],1,1,1+ix);
        for (unsigned int iy=0; iy<2; ++iy) {
          book(_h_pT_y[ix][iy], {0.,0.4,0.8,1.2,1.6,2.0,2.4});
          for(unsigned int iz=0; iz<6; ++iz) {
            book(_h_pT_y[ix][iy]->bin(iz+1), 5+3*iz+ix, 1, iy+1);
          }
          book(_h_pT[ix][iy], 2 +ix, 1, iy+1);
          book(_h_y [ix][iy], 23+ix, 1, iy+1);
          book(_h_r [ix][iy], "TMP/h_r_"+toString(ix)+"_"+toString(iy), refData(26,1,iy+1));
        }
        book(_h_pT_acc[ix], 29+ix, 1, 1);
      }
    }

    void findChildren(const Particle& p, Particles& mum, Particles& mup, unsigned int& nstable) {
      for (const Particle & child : p.children()) {
        if (child.pid()==PID::MUON) {
          mum.push_back(child);
          ++nstable;
        }
        else if (child.pid()==PID::ANTIMUON) {
          mup.push_back(child);
          ++nstable;
        }
        else if (child.pid()==PID::PHOTON) {
          continue;
        }
        else if (child.children().empty()) {
          ++nstable;
        }
        else {
          findChildren(child,mum,mup,nstable);
        }
      }
    }

    // from eqn 1 of paper
    bool acceptMuon(const Particle & p) const {
      const double abseta = p.abseta();
      const double xp = p.perp();
      if (abseta<0.8)      return xp>3.75;
      else if (abseta<1.6) return xp>3.5;
      else if (abseta<2.4) return xp>3.0;
      else return false;
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty()) _edges = _h_total[0]->xEdges();
      // Final state of unstable particles to get particle spectra
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");

      for (const Particle& p : ufs.particles(Cuts::pid==553 || Cuts::pid==100553 || Cuts::pid==200553)) {
        const double absrap = p.absrap();
        // rapidity cut
        if (absrap>2.4) continue;
        double xp = p.perp();
        unsigned int iloc=0;
        if      (p.pid()==   553) iloc=0;
        else if (p.pid()==100553) iloc=1;
        else if (p.pid()==200553) iloc=2;
        // acceptance corrected only hist
        if (absrap<1.2) {
          _h_pT_acc[iloc]->fill(xp);
        }
        // check if children muons and within acceptance
        unsigned int imin=1;
        // find the children
        Particles mum,mup;
        unsigned int nstable(0);
        findChildren(p,mum,mup,nstable);
        if (mup.size()==1 && mup.size()==1 && nstable==2) {
          if (acceptMuon(mup[0]) && acceptMuon(mum[0])) imin=0;
        }
        // fill the histos
        for (unsigned int ix=imin;ix<2;++ix) {
          _h_pT_y[iloc][ix] ->fill(absrap,xp);
          _h_pT  [iloc][ix]->fill(xp);
          _h_r   [iloc][ix]->fill(xp);
          if (xp<50.) _h_y[iloc][ix]->fill(absrap);
          if (ix==0) {
            _h_total[ix]->fill(_edges[iloc]);
          }
          else {
            _h_total[ix]->fill(_edges[iloc], _br[iloc]);
          }
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      double factor = crossSection() / nanobarn/ sumOfWeights();
      for (unsigned int ix=0; ix<2; ++ix) {
        // total cross sections, just the factor
        scale(_h_total[ix], factor);
        for (unsigned int iy=0; iy<3; ++iy) {
          double factor2 = factor;
          if (ix==1) factor2*=_br[iy];
          // pT integrated over y, just the factor
          scale(_h_pT[iy][ix], factor2);
          scale(_h_r[iy][ix], factor2);
          if (ix==1) scale(_h_pT_acc[iy],factor2);
          // not integrated over y, alsso undo y +/- folding
          scale(_h_y[iy][ix], 0.5*factor2);
          scale(_h_pT_y[iy][ix], 0.5*factor2);
          divByGroupWidth(_h_pT_y[iy][ix]);
        }
        // ratios
        Estimate1DPtr tmp;
        // ups 3/ ups1
        book(tmp, 26, 1, ix+1);
        divide(_h_r[2][ix], _h_r[0][ix], tmp);
        // ups 2/ ups1
        book(tmp, 27, 1, ix+1);
        divide(_h_r[1][ix], _h_r[0][ix], tmp);
        // ups 3/ ups2
        book(tmp, 28, 1, ix+1);
        divide(_h_r[2][ix], _h_r[1][ix], tmp);
      }
    }

    ///@}


    /// @name Histograms
    ///@{
    Histo1DGroupPtr _h_pT_y[3][2];
    BinnedHistoPtr<string> _h_total[2];
    Histo1DPtr _h_pT[3][2],_h_pT_acc[3],_h_y[3][2],_h_r[3][2];
    vector<double> _br{0.0248,0.0193,0.0218};
    vector<string> _edges;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(CMS_2013_I1225274);

}