Rivet analyses

title: CMS_2013_I1185414

Polarization of $\Upsilon(1,2,3)$ at 7 TeV

Experiment: CMS (LHC)

Inspire ID: 1185414

Status: VALIDATED

Authors: - Peter Richardson

References: - Phys.Rev.Lett. 110 (2013) 081802, 2013 - arXiv: 1209.2922 - Expt page: CMS-BPH-11-023

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - Upsilon production

Measurement of the polarization of $\Upsilon(1,2,3)$ at 7 TeV by CMS

Source code:CMS_2013_I1185414.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief Upsilon polarization at 7 TeV class CMS_2013_I1185414 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1185414);


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

/// Book histograms and initialise projections before the run
void init() {
  // projections
  declare(Beam(), "Beams");
  declare(UnstableParticles(), "UFS");
  // book profile hists for the moments
  // loop upsilon states
  for (unsigned int iups = 0; iups < 3; ++iups) {
    // rapidity intervals
    for (unsigned int iy = 0; iy < 2; ++iy) {
      // frame defns
      for (unsigned int iframe = 0; iframe < 3; ++iframe) {
        // 3 different moments
        for (unsigned int iobs = 0; iobs < 3; ++iobs) {
          const string name = "TMP/POL_" + toString(iups) + "_" + toString(iy) + "_" + toString(iframe)
              + "_" + toString(iobs);
          book(_p_Ups[iups][iy][iframe][iobs], name, refData(1, 1, 1));
        }
      }
    }
  }
}

void findDecayProducts(const Particle& mother, unsigned int& nstable, Particles& mup, Particles& mum) {
  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 (!p.children().empty()) {
      findDecayProducts(p, nstable, mup, mum);
    }
    else {
      ++nstable;
    }
  }
}


/// Perform the per-event analysis
void analyze(const Event& event) {
  // find the beams
  const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
  // 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 or Cuts::pid == 100553 or Cuts::pid == 200553)) {
    // check mu+mu- decay and find muons
    unsigned int nstable = 0;
    Particles mup, mum;
    findDecayProducts(p, nstable, mup, mum);
    if (mup.size() != 1 || mum.size() != 1 || nstable != 2) continue;
    // pT and rapidity
    const double absrap = p.absrap();
    const double xp = p.perp();
    if (absrap > 1.2 || xp < 10. || xp > 50.) continue;
    // type of upsilon
    unsigned int iups = p.pid() / 100000;
    // rapidity interval
    unsigned int iy = absrap > 0.6;
    // first the CS frame
    // first boost so upslion momentum =0 in z direction
    Vector3 beta = p.mom().betaVec();
    beta.setX(0.);
    beta.setY(0.);
    LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(beta);
    FourMomentum pp = boost.transform(p.mom());
    // and then transverse so pT=0
    beta = pp.betaVec();
    LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(beta);
    // get all the momenta in this frame
    Vector3 muDirn = boost2.transform(boost.transform(mup[0].mom())).p3().unit();
    FourMomentum p1 = boost2.transform(boost.transform(beams.first.mom()));
    FourMomentum p2 = boost2.transform(boost.transform(beams.second.mom()));
    if (beams.first.mom().z() < 0.) swap(p1, p2);
    if (p.rapidity() < 0.) swap(p1, p2);
    Vector3 axisy = (p1.p3().cross(p2.p3())).unit();
    Vector3 axisz(0., 0., 1.);
    Vector3 axisx = axisy.cross(axisz);
    double cTheta = axisz.dot(muDirn);
    double cPhi = axisx.dot(muDirn);
    // fill the moments
    _p_Ups[iups][iy][0][0]->fill(xp, 1.25 * (3. * sqr(cTheta) - 1.));
    _p_Ups[iups][iy][0][1]->fill(xp, 1.25 * (1. - sqr(cTheta)) * (2. * sqr(cPhi) - 1.));
    _p_Ups[iups][iy][0][2]->fill(xp, 2.5 * cTheta * sqrt(1. - sqr(cTheta)) * cPhi);
    // now for the HX frame
    beta = p.mom().betaVec();
    boost = LorentzTransform::mkFrameTransformFromBeta(beta);
    axisz = pp.p3().unit();
    axisx = axisy.cross(axisz);
    cTheta = axisz.dot(muDirn);
    cPhi = axisx.dot(muDirn);
    // fill the moments
    _p_Ups[iups][iy][1][0]->fill(xp, 1.25 * (3. * sqr(cTheta) - 1.));
    _p_Ups[iups][iy][1][1]->fill(xp, 1.25 * (1. - sqr(cTheta)) * (2. * sqr(cPhi) - 1.));
    _p_Ups[iups][iy][1][2]->fill(xp, 2.5 * cTheta * sqrt(1. - sqr(cTheta)) * cPhi);
    // then PX
    axisz = (p1.p3().unit() + p2.p3().unit()).unit();
    axisx = axisy.cross(axisz);
    cTheta = axisz.dot(muDirn);
    cPhi = axisx.dot(muDirn);
    // fill the moments
    _p_Ups[iups][iy][2][0]->fill(xp, 1.25 * (3. * sqr(cTheta) - 1.));
    _p_Ups[iups][iy][2][1]->fill(xp, 1.25 * (1. - sqr(cTheta)) * (2. * sqr(cPhi) - 1.));
    _p_Ups[iups][iy][2][2]->fill(xp, 2.5 * cTheta * sqrt(1. - sqr(cTheta)) * cPhi);
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  // Loop over states
  for (unsigned int iups = 0; iups < 3; ++iups) {
    // Loop over rapidity ranges
    for (unsigned int iy = 0; iy < 2; ++iy) {
      // Loop over frame definition
      for (unsigned int iframe = 0; iframe < 3; ++iframe) {
        // base no for the ihistos in rivet
        unsigned int ibase = 24 * iups + iy + 8 * iframe;
        // book scatters
        Estimate1DPtr lTheta, lPhi, lThetaPhi, lTilde;
        book(lTheta, ibase + 1, 1, 1);
        book(lPhi, ibase + 3, 1, 1);
        book(lThetaPhi, ibase + 5, 1, 1);
        book(lTilde, ibase + 7, 1, 1);
        // histos for the moments
        Profile1DPtr moment[3];
        for (unsigned int ix = 0; ix < 3; ++ix) {
          moment[ix] = _p_Ups[iups][iy][iframe][ix];
        }
        // loop over bins
        for (unsigned int ibin = 1; ibin <= moment[0]->bins().size(); ++ibin) {
          // extract moments and errors
          double val[3], err[3];
          // m1 = lTheta/(3+lTheta), m2 = lPhi/(3+lTheta), m3 = lThetaPhi/(3+lTheta)
          for (unsigned int ix = 0; ix < 3; ++ix) {
            bool has0 = moment[ix]->bin(ibin).numEntries() > 0
                && moment[ix]->bin(ibin).effNumEntries() > 0;
            bool has1 = moment[ix]->bin(ibin).numEntries() > 1
                && moment[ix]->bin(ibin).effNumEntries() > 1;
            val[ix] = has0 ? moment[ix]->bin(ibin).mean(2) : 0.;
            err[ix] = has1 ? moment[ix]->bin(ibin).stdErr(2) : 0.;
          }
          // values of the lambdas and their errors
          double l1 = 3. * val[0] / (1. - val[0]);
          double l2 = (3. + l1) * val[1];
          // fill the scatters
          lTheta->bin(ibin).setVal(l1);
          lTheta->bin(ibin).setErr(3. / sqr(1. - val[0]) * err[0]);
          lPhi->bin(ibin).setVal(l2);
          lPhi->bin(ibin).setErr(3. / sqr(1. - val[0])
                                 * sqrt(sqr(err[0] * val[1]) + sqr(err[1] * (1. - val[0]))));
          lThetaPhi->bin(ibin).setVal((3. + l1) * val[2]);
          lThetaPhi->bin(ibin).setErr(3. / sqr(1. - val[0])
                                      * sqrt(sqr(err[0] * val[1]) + sqr(err[1] * (1. - val[0]))));
          lTilde->bin(ibin).setVal((l1 + 3. * l2) / (1. - l2));
          lTilde->bin(ibin).setErr(3. / sqr(1. - val[0] - 3 * val[1])
                                   * sqrt(sqr(err[0]) + 9. * sqr(err[1])));
        }
      }
    }
  }
}

/// @}


/// @name Histograms
/// @{
Profile1DPtr _p_Ups[3][2][3][3];
/// @}

};

RIVET_DECLARE_PLUGIN(CMS_2013_I1185414);

} ```