Rivet analyses

title: CDF_2010_I849042

CDF Run 2 underlying event in Drell-Yan

Experiment: CDF (Tevatron Run 2)

Inspire ID: 849042

Status: VALIDATED

Authors: - Hendrik Hoeth

References: - Phys.Rev.D82:034001,2010

Beams: p- p+

Beam energies: (980.0, 980.0)GeV

Run details: - ppbar collisions at 1960 GeV. * Drell-Yan events with $Z/\gamma -> e e$ and $Z/\gamma -> \mu\mu$. * A mass cut $m_{ll} > 70~\text{GeV}$ can be applied on generator level. * Particles with $c \tau > 10~\text{mm}$ should be set stable.

Deepak Kar and Rick Field's measurement of the underlying event in Drell-Yan events. $Z -> ee$ and $Z -> \mu\mu$ events are selected using a $Z$ mass window cut between 70 and 110~GeV. Toward'',away'' and ``transverse'' regions are defined in the same way as in the original (2001) CDF underlying event analysis. The reconstructed $Z$ defines the $\phi$ direction of the toward region. The leptons are ignored after the $Z$ has been reconstructed. Thus the region most sensitive to the underlying event is the toward region (the recoil jet is boosted into the away region).

Source code:CDF_2010_I849042.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/ChargedFinalState.hh"

include "Rivet/Projections/ChargedLeptons.hh"

include "Rivet/Projections/FastJets.hh"

namespace Rivet {

/// @brief CDF Run II underlying event in Drell-Yan /// /// @author Hendrik Hoeth /// /// Measurement of the underlying event in Drell-Yan /// \f$ Z/\gamma^ \to e^+ e^- \f$ and /// \f$ Z/\gamma^ \to \mu^+ \mu^- \f$ events. The reconstructed /// Z defines the \f$ \phi \f$ orientation. A Z mass window cut is applied. /// /// @par Run conditions /// /// @arg \f$ \sqrt{s} = \f$ 1960 GeV /// @arg produce Drell-Yan events /// @arg Set particles with c*tau > 10 mm stable /// @arg Z decay mode: Z -> e+e- and Z -> mu+mu- /// @arg gamma decay mode: gamma -> e+e- and gamma -> mu+mu- /// @arg minimum invariant mass of the fermion pair coming from the Z/gamma: 70 GeV class CDF_2010_I849042 : public Analysis { public:

RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2010_I849042);


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

void init() {

  _mode = 0;
  if (getOption("MODE") == "DY")
    _mode = 1;
  else if (getOption("MODE") == "QCD")
    _mode = 2;

  // Set up projections
  const ChargedFinalState cfs(Cuts::abseta < 1.0 && Cuts::pT >= 0.5 * GeV);
  const ChargedFinalState clfs(Cuts::abseta < 1.0 && Cuts::pT >= 20 * GeV);
  declare(cfs, "CFS");
  declare(ChargedLeptons(clfs), "CL");

  // Final state for the jet finding
  const FinalState fsj(Cuts::abseta < 4.0);
  declare(fsj, "FSJ");
  declare(FastJets(fsj, JetAlg::CDFMIDPOINT, 0.7), "MidpointJets");

  // Book histograms
  if (_mode == 0 || _mode == 1) {
    string mode("Z");
    book(_p[mode + "tnchg"], 1, 1, 1);
    book(_p[mode + "pnchg"], 1, 1, 2);
    book(_p[mode + "anchg"], 1, 1, 3);
    book(_p[mode + "pmaxnchg"], 2, 1, 1);
    book(_p[mode + "pminnchg"], 2, 1, 2);
    book(_p[mode + "pdifnchg"], 2, 1, 3);
    book(_p[mode + "tcptsum"], 3, 1, 1);
    book(_p[mode + "pcptsum"], 3, 1, 2);
    book(_p[mode + "acptsum"], 3, 1, 3);
    book(_p[mode + "pmaxcptsum"], 4, 1, 1);
    book(_p[mode + "pmincptsum"], 4, 1, 2);
    book(_p[mode + "pdifcptsum"], 4, 1, 3);
    book(_p[mode + "tcptave"], 5, 1, 1);
    book(_p[mode + "pcptave"], 5, 1, 2);
    book(_p[mode + "tcptmax"], 6, 1, 1);
    book(_p[mode + "pcptmax"], 6, 1, 2);
    book(_p[mode + "zptvsnchg"], 7, 1, 1);
    book(_p[mode + "cptavevsnchg"], 8, 1, 1);
    book(_p[mode + "cptavevsnchgsmallzpt"], 9, 1, 1);
  }

  if (_mode == 0 || _mode == 2) {
    string mode("QCD");
    book(_p[mode + "tnchg"], 10, 1, 1);
    book(_p[mode + "pnchg"], 10, 1, 2);
    book(_p[mode + "anchg"], 10, 1, 3);
    book(_p[mode + "pmaxnchg"], 11, 1, 1);
    book(_p[mode + "pminnchg"], 11, 1, 2);
    book(_p[mode + "pdifnchg"], 11, 1, 3);
    book(_p[mode + "tcptsum"], 12, 1, 1);
    book(_p[mode + "pcptsum"], 12, 1, 2);
    book(_p[mode + "acptsum"], 12, 1, 3);
    book(_p[mode + "pmaxcptsum"], 13, 1, 1);
    book(_p[mode + "pmincptsum"], 13, 1, 2);
    book(_p[mode + "pdifcptsum"], 13, 1, 3);
    book(_p[mode + "pcptave"], 14, 1, 1);
    book(_p[mode + "pcptmax"], 15, 1, 1);
  }
}


/// Do the analysis
void analyze(const Event& event) {

  if (_mode == 0 || _mode == 1) doDYanalysis(event);
  if (_mode == 0 || _mode == 2) doQCDanalysis(event);
}


void doDYanalysis(const Event& e) {

  const string pre("Z");

  const FinalState& fs = apply<FinalState>(e, "CFS");
  const size_t numParticles = fs.particles().size();

  // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
  if (numParticles < 1) {
    MSG_DEBUG("Failed multiplicity cut");
    vetoEvent;
  }

  // Get the leptons
  const Particles& leptons = apply<ChargedLeptons>(e, "CL").chargedLeptons();

  // We want exactly two leptons of the same flavour.
  MSG_DEBUG("lepton multiplicity = " << leptons.size());
  if (leptons.size() != 2 || leptons[0].pid() != -leptons[1].pid()) vetoEvent;

  // Lepton pT > 20 GeV
  if (leptons[0].pT() / GeV <= 20 || leptons[1].pT() / GeV <= 20) vetoEvent;

  // Lepton pair should have an invariant mass between 70 and 110 and |eta| < 6
  const FourMomentum dilepton = leptons[0].momentum() + leptons[1].momentum();
  if (!inRange(dilepton.mass() / GeV, 70., 110.) || fabs(dilepton.eta()) >= 6) vetoEvent;
  MSG_DEBUG("Dilepton mass = " << dilepton.mass() / GeV << " GeV");
  MSG_DEBUG("Dilepton pT   = " << dilepton.pT() / GeV << " GeV");

  // Calculate the observables
  size_t numToward(0), numAway(0);
  long int numTrans1(0), numTrans2(0);
  double ptSumToward(0.0), ptSumTrans1(0.0), ptSumTrans2(0.0), ptSumAway(0.0);
  double ptMaxToward(0.0), ptMaxTrans1(0.0), ptMaxTrans2(0.0), ptMaxAway(0.0);
  const double phiZ = dilepton.azimuthalAngle();
  const double pTZ = dilepton.pT();
  /// @todo Replace with for
  for (Particles::const_iterator p = fs.particles().begin(); p != fs.particles().end(); ++p) {
    // Don't use the leptons
    /// @todo Replace with PID::isLepton
    if (abs(p->pid()) < 20) continue;

    const double dPhi = deltaPhi(p->momentum().phi(), phiZ);
    const double pT = p->pT();
    double rotatedphi = p->momentum().phi() - phiZ;
    while (rotatedphi < 0) rotatedphi += 2 * PI;

    if (dPhi < PI / 3.0) {
      ptSumToward += pT;
      ++numToward;
      if (pT > ptMaxToward) ptMaxToward = pT;
    }
    else if (dPhi < 2 * PI / 3.0) {
      if (rotatedphi <= PI) {
        ptSumTrans1 += pT;
        ++numTrans1;
        if (pT > ptMaxTrans1) ptMaxTrans1 = pT;
      }
      else {
        ptSumTrans2 += pT;
        ++numTrans2;
        if (pT > ptMaxTrans2) ptMaxTrans2 = pT;
      }
    }
    else {
      ptSumAway += pT;
      ++numAway;
      if (pT > ptMaxAway) ptMaxAway = pT;
    }
    // We need to subtract the two leptons from the number of particles to get the correct multiplicity
    _p[pre + "cptavevsnchg"]->fill(numParticles - 2, pT);
    if (pTZ < 10) _p[pre + "cptavevsnchgsmallzpt"]->fill(numParticles - 2, pT);
  }

  // Fill the histograms
  _p[pre + "tnchg"]->fill(pTZ, numToward / (4 * PI / 3));
  _p[pre + "pnchg"]->fill(pTZ, (numTrans1 + numTrans2) / (4 * PI / 3));
  _p[pre + "pmaxnchg"]->fill(pTZ, (numTrans1 > numTrans2 ? numTrans1 : numTrans2) / (2 * PI / 3));
  _p[pre + "pminnchg"]->fill(pTZ, (numTrans1 < numTrans2 ? numTrans1 : numTrans2) / (2 * PI / 3));
  _p[pre + "pdifnchg"]->fill(pTZ, abs(numTrans1 - numTrans2) / (2 * PI / 3));
  _p[pre + "anchg"]->fill(pTZ, numAway / (4 * PI / 3));

  _p[pre + "tcptsum"]->fill(pTZ, ptSumToward / (4 * PI / 3));
  _p[pre + "pcptsum"]->fill(pTZ, (ptSumTrans1 + ptSumTrans2) / (4 * PI / 3));
  _p[pre + "pmaxcptsum"]->fill(pTZ,
                               (ptSumTrans1 > ptSumTrans2 ? ptSumTrans1 : ptSumTrans2) / (2 * PI / 3));
  _p[pre + "pmincptsum"]->fill(pTZ,
                               (ptSumTrans1 < ptSumTrans2 ? ptSumTrans1 : ptSumTrans2) / (2 * PI / 3));
  _p[pre + "pdifcptsum"]->fill(pTZ, fabs(ptSumTrans1 - ptSumTrans2) / (2 * PI / 3));
  _p[pre + "acptsum"]->fill(pTZ, ptSumAway / (4 * PI / 3));

  if (numToward > 0) {
    _p[pre + "tcptave"]->fill(pTZ, ptSumToward / numToward);
    _p[pre + "tcptmax"]->fill(pTZ, ptMaxToward);
  }
  if ((numTrans1 + numTrans2) > 0) {
    _p[pre + "pcptave"]->fill(pTZ, (ptSumTrans1 + ptSumTrans2) / (numTrans1 + numTrans2));
    _p[pre + "pcptmax"]->fill(pTZ, (ptMaxTrans1 > ptMaxTrans2 ? ptMaxTrans1 : ptMaxTrans2));
  }

  // We need to subtract the two leptons from the number of particles to get the correct multiplicity
  _p[pre + "zptvsnchg"]->fill(numParticles - 2, pTZ);
}


void doQCDanalysis(const Event& e) {

  const string pre("QCD");

  const FinalState& fsj = apply<FinalState>(e, "FSJ");
  if (fsj.particles().size() < 1) {
    MSG_DEBUG("Failed multiplicity cut");
    vetoEvent;
  }

  const Jets& jets = apply<FastJets>(e, "MidpointJets").jetsByPt();
  MSG_DEBUG("Jet multiplicity = " << jets.size());

  // We require the leading jet to be within |eta|<2
  if (jets.size() < 1 || fabs(jets[0].eta()) >= 2) {
    MSG_DEBUG("Failed leading jet cut");
    vetoEvent;
  }

  const double jetphi = jets[0].phi();
  const double jeteta = jets[0].eta();
  const double jetpT = jets[0].pT();
  MSG_DEBUG("Leading jet: pT = " << jetpT << ", eta = " << jeteta << ", phi = " << jetphi);

  // Get the final states to work with for filling the distributions
  const FinalState& cfs = apply<ChargedFinalState>(e, "CFS");

  size_t numToward(0), numAway(0);
  long int numTrans1(0), numTrans2(0);
  double ptSumToward(0.0), ptSumTrans1(0.0), ptSumTrans2(0.0), ptSumAway(0.0);
  double ptMaxOverall(0.0), ptMaxToward(0.0), ptMaxTrans1(0.0), ptMaxTrans2(0.0), ptMaxAway(0.0);

  // Calculate all the charged stuff
  for (const Particle& p : cfs.particles()) {
    const double dPhi = deltaPhi(p.phi(), jetphi);
    const double pT = p.pT();
    const double phi = p.phi();
    double rotatedphi = phi - jetphi;
    while (rotatedphi < 0) rotatedphi += 2 * PI;

    if (pT > ptMaxOverall) {
      ptMaxOverall = pT;
    }

    if (dPhi < PI / 3.0) {
      ptSumToward += pT;
      ++numToward;
      if (pT > ptMaxToward) ptMaxToward = pT;
    }
    else if (dPhi < 2 * PI / 3.0) {
      if (rotatedphi <= PI) {
        ptSumTrans1 += pT;
        ++numTrans1;
        if (pT > ptMaxTrans1) ptMaxTrans1 = pT;
      }
      else {
        ptSumTrans2 += pT;
        ++numTrans2;
        if (pT > ptMaxTrans2) ptMaxTrans2 = pT;
      }
    }
    else {
      ptSumAway += pT;
      ++numAway;
      if (pT > ptMaxAway) ptMaxAway = pT;
    }
  } // end charged particle loop

  // Fill the histograms
  _p[pre + "tnchg"]->fill(jetpT / GeV, numToward / (4 * PI / 3));
  _p[pre + "pnchg"]->fill(jetpT / GeV, (numTrans1 + numTrans2) / (4 * PI / 3));
  _p[pre + "pmaxnchg"]->fill(jetpT / GeV, (numTrans1 > numTrans2 ? numTrans1 : numTrans2) / (2 * PI / 3));
  _p[pre + "pminnchg"]->fill(jetpT / GeV, (numTrans1 < numTrans2 ? numTrans1 : numTrans2) / (2 * PI / 3));
  _p[pre + "pdifnchg"]->fill(jetpT / GeV, abs(numTrans1 - numTrans2) / (2 * PI / 3));
  _p[pre + "anchg"]->fill(jetpT / GeV, numAway / (4 * PI / 3));

  _p[pre + "tcptsum"]->fill(jetpT / GeV, ptSumToward / GeV / (4 * PI / 3));
  _p[pre + "pcptsum"]->fill(jetpT / GeV, (ptSumTrans1 + ptSumTrans2) / GeV / (4 * PI / 3));
  _p[pre + "pmaxcptsum"]->fill(
      jetpT / GeV, (ptSumTrans1 > ptSumTrans2 ? ptSumTrans1 : ptSumTrans2) / GeV / (2 * PI / 3));
  _p[pre + "pmincptsum"]->fill(
      jetpT / GeV, (ptSumTrans1 < ptSumTrans2 ? ptSumTrans1 : ptSumTrans2) / GeV / (2 * PI / 3));
  _p[pre + "pdifcptsum"]->fill(jetpT / GeV, fabs(ptSumTrans1 - ptSumTrans2) / GeV / (2 * PI / 3));
  _p[pre + "acptsum"]->fill(jetpT / GeV, ptSumAway / GeV / (4 * PI / 3));

  if ((numTrans1 + numTrans2) > 0) {
    _p[pre + "pcptave"]->fill(jetpT / GeV, (ptSumTrans1 + ptSumTrans2) / GeV / (numTrans1 + numTrans2));
    _p[pre + "pcptmax"]->fill(jetpT / GeV, (ptMaxTrans1 > ptMaxTrans2 ? ptMaxTrans1 : ptMaxTrans2) / GeV);
  }
}


// void finalize() {    }

/// @}

private:

size_t _mode;
map<string, Profile1DPtr> _p;

};

RIVET_DECLARE_ALIASED_PLUGIN(CDF_2010_I849042, CDF_2010_S8591881);

} ```

Aliases: - CDF_2010_S8591881