Rivet analyses

title: ATLAS_2013_I1204994

$\Upsilon(1,2,3S)$ production at 7 TeV

Experiment: ATLAS (LHC)

Inspire ID: 1204994

Status: VALIDATED

Authors: - Peter Richardsonm

References: - Expt page: ATLAS-BPHY-2011-06 - arXiv: 1211.7255 - Phys.Rev.D 87 (2013) 052004, 2013

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - Upsilon(1,2,3S) production

Measurement of the double differential cross section for $\Upsilon(1,2,3S)$ production at 7 TeV by the ATLAS collaboration

Source code:ATLAS_2013_I1204994.cc

```c++ // -- 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);

} ```