Rivet analyses

title: ATLAS_2016_I1419652

Track-based minimum bias at 13 TeV in ATLAS

Experiment: ATLAS (LHC)

Inspire ID: 1419652

Status: VALIDATED

Authors: - Roman Lysak

References: - Expt page: ATLAS-STDM-2015-02 - arXiv: 1602.01633 - DOI: 10.1016/j.physletb.2016.04.050 - Physics Letters B (2016), Vol. 758, pp. 67-88

Beams: p+ p+

Beam energies: (6500.0, 6500.0)GeV

Run details: - $pp$ QCD interactions at 13 TeV. Diffractive events should be included. Multiple kinematic cuts should not be required.

Measurements from proton-proton collisions at centre-of-mass energy of $\sqrt{s} = 13$ TeV recorded with the ATLAS detector at the LHC. Events were collected using a single-arm minimum-bias trigger. The charged-particle multiplicity, its dependence on transverse momentum and pseudorapidity and the relationship between the mean transverse momentum and charged-particle multiplicity are measured. The observed distributions are corrected to well-defined phase-space regions ($p_\text{T} > 500$ MeV and $|\eta| < 2.5$ of the particles), using model-independent corrections.

Source code:ATLAS_2016_I1419652.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {

class ATLAS_2016_I1419652 : public Analysis { public:

/// Particle types included
enum PartTypes { k_NoStrange, k_AllCharged, kNPartTypes };

/// Phase space regions
enum RegionID { k_pt500_nch1_eta25, k_pt500_nch1_eta08, kNregions };

/// Nch cut for each region
const static int nchCut[kNregions];


/// Default constructor
ATLAS_2016_I1419652()
    : Analysis("ATLAS_2016_I1419652") { }


/// Initialization, called once before running
void init() {

  // Projections
  const ChargedFinalState cfs500_25((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >= 500.0 * MeV));
  declare(cfs500_25, "CFS500_25");

  const ChargedFinalState cfs500_08((Cuts::etaIn(-0.8, 0.8) && Cuts::pT >= 500.0 * MeV));
  declare(cfs500_08, "CFS500_08");

  for (int iT = 0; iT < kNPartTypes; ++iT) {
    for (int iR = 0; iR < kNregions; ++iR) {
      size_t offset = 8 * iR + 4 * iT;
      book(_sumW[iT][iR], "_sumW" + to_str(iT) + to_str(iR));
      book(_hist_eta[iT][iR], offset + 3, 1, 1);
      book(_hist_pt[iT][iR], offset + 4, 1, 1);
      book(_hist_nch[iT][iR], offset + 5, 1, 1);
      book(_hist_ptnch[iT][iR], offset + 6, 1, 1);
    }
  }
}


void analyze(const Event& event) {
  string fsName;
  for (int iR = 0; iR < kNregions; ++iR) {
    switch (iR) {
      case k_pt500_nch1_eta25: fsName = "CFS500_25"; break;
      case k_pt500_nch1_eta08: fsName = "CFS500_08"; break;
    }

    const ChargedFinalState& cfs = apply<ChargedFinalState>(event, fsName);

    /// What's the benefit in separating this code which is only called from one place?!
    fillPtEtaNch(cfs, iR);
  }
}


void finalize() {
  // Standard histograms
  for (int iT = 0; iT < kNPartTypes; ++iT) {
    for (int iR = 0; iR < kNregions; ++iR) {

      double etaRangeSize = -999.0; //intentionally crazy
      switch (iR) {
        case k_pt500_nch1_eta25: etaRangeSize = 5.0; break;
        case k_pt500_nch1_eta08: etaRangeSize = 1.6; break;
        default: etaRangeSize = -999.0; break; //intentionally crazy
      }

      if (_sumW[iT][iR]->val() > 0) {
        scale(_hist_nch[iT][iR], 1.0 / *_sumW[iT][iR]);
        scale(_hist_pt[iT][iR], 1.0 / dbl(*_sumW[iT][iR]) / TWOPI / etaRangeSize);
        scale(_hist_eta[iT][iR], 1.0 / *_sumW[iT][iR]);
      }
      else {
        MSG_WARNING("Sum of weights is zero (!) in type/region: " << iT << " " << iR);
      }
    }
  }
}


/// Helper for collectively filling Nch, pT, eta, and pT vs. Nch histograms
void fillPtEtaNch(const ChargedFinalState& cfs, int iRegion) {

  // Get number of particles
  int nch[kNPartTypes];
  int nch_noStrange = 0;
  for (const Particle& p : cfs.particles()) {
    PdgId pdg = p.abspid();
    if (pdg == 3112 || // Sigma-
        pdg == 3222 || // Sigma+
        pdg == 3312 || // Xi-
        pdg == 3334)   // Omega-
      continue;
    nch_noStrange++;
  }
  nch[k_AllCharged] = cfs.size();
  nch[k_NoStrange] = nch_noStrange;

  // Skip if event fails cut for all charged (noStrange will always be less)
  if (nch[k_AllCharged] < nchCut[iRegion]) return;

  // Fill event weight info
  _sumW[k_AllCharged][iRegion]->fill();
  if (nch[k_NoStrange] >= nchCut[iRegion]) {
    _sumW[k_NoStrange][iRegion]->fill();
  }

  // Fill nch
  _hist_nch[k_AllCharged][iRegion]->fill(nch[k_AllCharged]);
  if (nch[k_NoStrange] >= nchCut[iRegion]) {
    _hist_nch[k_NoStrange][iRegion]->fill(nch[k_NoStrange]);
  }

  // Loop over particles, fill pT, eta and ptnch
  for (const Particle& p : cfs.particles()) {
    const double pt = p.pT() / GeV;
    const double eta = p.eta();
    _hist_pt[k_AllCharged][iRegion]->fill(pt, 1.0 / pt);
    _hist_eta[k_AllCharged][iRegion]->fill(eta);
    _hist_ptnch[k_AllCharged][iRegion]->fill(nch[k_AllCharged], pt);

    // Make sure nch cut is passed for nonStrange particles!
    if (nch[k_NoStrange] >= nchCut[iRegion]) {
      PdgId pdg = p.abspid();
      if (pdg == 3112 || // Sigma-
          pdg == 3222 || // Sigma+
          pdg == 3312 || // Xi-
          pdg == 3334)   // Omega-
        continue;
      // Here we don't have strange particles anymore
      _hist_pt[k_NoStrange][iRegion]->fill(pt, 1.0 / pt);
      _hist_eta[k_NoStrange][iRegion]->fill(eta);
      _hist_ptnch[k_NoStrange][iRegion]->fill(nch[k_NoStrange], pt);
    }
  }
}

private:

CounterPtr _sumW[kNPartTypes][kNregions];

Histo1DPtr _hist_nch[kNPartTypes][kNregions];
Histo1DPtr _hist_pt[kNPartTypes][kNregions];
Histo1DPtr _hist_eta[kNPartTypes][kNregions];
Profile1DPtr _hist_ptnch[kNPartTypes][kNregions];

};

// Constants: pT & eta regions const int ATLAS_2016_I1419652::nchCut[] = {1, 1};

RIVET_DECLARE_PLUGIN(ATLAS_2016_I1419652);

} ```