Rivet analyses

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

namespace Rivet {


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

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2013_I1204994);


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

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

      declare(UnstableParticles(Cuts::pid==553 || Cuts::pid==100553 || Cuts::pid==200553), "UFS");

      // branching ratios
      _br = {0.0248,0.0193,0.0218};

      // total cross section
      book(_h_total,1,1,1);

      // double differential
      const vector<double> ybins = {0.,1.2,2.25};
      for(unsigned int ix=0;ix<3;++ix) {
        book(_h_Upsilon_fid_y[ix],5,1,ix+1);
        book(_h_Upsilon_y    [ix],9,1,ix+1);
        book(_h_Upsilon_y_r  [ix],"TMP/Ups"+toString(ix),refData(12,1,1));
        book(_h_Upsilon_fid[ix], ybins);
        book(_h_Upsilon[ix],     ybins);
        book(_h_Upsilon_r[ix],   ybins);
        for (size_t iy=1; iy < ybins.size(); ++iy) {
          book(_h_Upsilon_fid[ix]->bin(iy), ix+2,1,iy);
          book(_h_Upsilon[ix]->bin(iy), ix+6,1,iy);
          book(_h_Upsilon_r[ix]->bin(iy), "TMP/Ups"+toString(ix)+"_"+toString(iy), refData(11,1,iy));
        }
      }
    }

    void findDecayProducts(const Particle & mother, unsigned int & nstable,
                           Particles & mup, Particles & mum, unsigned int & ngamma) {
      for (const Particle & p : mother.children()) {
        int id = p.pid();
        if (id == PID::MUON ) {
          ++nstable;
          mum.push_back(p);
        }
        else if (id == PID::ANTIMUON) {
          ++nstable;
          mup.push_back(p);
        }
        else if (id == PID::PI0 || id == PID::K0S || id == PID::K0L ) {
          ++nstable;
        }
        else if (id == PID::GAMMA && p.children().empty() ) {
          ++ngamma;
          ++nstable;
        }
        else if ( !p.children().empty() ) {
          findDecayProducts(p, nstable, mup, mum,ngamma);
        }
        else
          ++nstable;
      }
    }

    /// 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");

      for (const Particle& p : ufs.particles()) {
        // check if in fiducal region
        unsigned int nstable=0, ngamma=0;
        Particles mup,mum;
        findDecayProducts(p,nstable,mup,mum,ngamma);
        bool fiducial = false;
        if (mup.size()==1 && mum.size()==1 && nstable==ngamma+2) {
          fiducial = true;
          if(mup[0].perp()<4. || mup[0].abseta()>2.4) fiducial = false;
          if(mum[0].perp()<4. || mum[0].abseta()>2.4) fiducial = false;
        }
        // pT and rapidity
        double absrap = p.absrap();
        double xp = p.perp();
        // type of upsilon
        unsigned int iups=p.pid()/100000;
        if (fiducial) {
          if (xp<70.) _h_Upsilon_fid_y[iups]->fill(absrap);
          _h_Upsilon_fid[iups]->fill(absrap,xp);
        }
        if (xp<70.) {
          _h_Upsilon_y[iups]->fill(absrap);
          _h_Upsilon_y_r[iups]->fill(absrap);
        }
        _h_Upsilon[iups]->fill(absrap,xp);
        _h_Upsilon_r[iups]->fill(absrap,xp);
        if (absrap<2.25 && xp<70) _h_total->fill(_edges[iups],_br[iups]);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {

      scale(_h_total,crossSection()/nanobarn/sumOfWeights());

      for (unsigned int ix=0; ix<3; ++ix) {
        scale(_h_Upsilon_fid_y[ix], 0.5*crossSection()/picobarn/sumOfWeights());
        scale(_h_Upsilon_y[ix],     0.5*_br[ix]*crossSection()/picobarn/sumOfWeights());
        scale(_h_Upsilon_y_r[ix],   0.5*_br[ix]*crossSection()/picobarn/sumOfWeights());

        scale(_h_Upsilon_fid[ix], 0.5*crossSection()/femtobarn/sumOfWeights());
        scale(_h_Upsilon[ix],     0.5*_br[ix]*crossSection()/femtobarn/sumOfWeights());
        scale(_h_Upsilon_r[ix],   0.5*_br[ix]*crossSection()/femtobarn/sumOfWeights());
        divByGroupWidth(_h_Upsilon_fid[ix]);
        divByGroupWidth(_h_Upsilon[ix]);
        divByGroupWidth(_h_Upsilon_r[ix]);
      }

      // ratios
      for (size_t iy=1; iy<_h_Upsilon[0]->numBins()+1; ++iy) {
        Estimate1DPtr tmp;
        book(tmp,10,1,iy);
        divide(_h_Upsilon_r[1]->bin(iy), _h_Upsilon_r[0]->bin(iy), tmp);
        book(tmp,11,1,iy);
        divide(_h_Upsilon_r[2]->bin(iy), _h_Upsilon_r[0]->bin(iy), tmp);
      }
      for (unsigned int iy=0; iy<2; ++iy) {
        Estimate1DPtr tmp;
        book(tmp,12,1,iy+1);
        divide(_h_Upsilon_y_r[iy+1], _h_Upsilon_y_r[0], tmp);
      }
    }
    /// @}


    /// @name Histograms
    /// @{
    BinnedHistoPtr<string> _h_total;
    Histo1DPtr _h_Upsilon_fid_y[3],_h_Upsilon_y[3],_h_Upsilon_y_r[3];
    Histo1DGroupPtr _h_Upsilon[3],_h_Upsilon_fid[3],_h_Upsilon_r[3];
    vector<double> _br;
    vector<string> _edges = { "P P --> UPSILON(1S) < MU+ MU- > X",
                              "P P --> UPSILON(2S) < MU+ MU- > X",
                              "P P --> UPSILON(3S) < MU+ MU- > X" };
    /// @}


  };


  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1204994);

}