rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2018_I1708620

Energy density vs pseudorapidity in proton-proton collisions at $\sqrt{s} = 13$ TeV
Experiment: CMS (LHC)
Inspire ID: 1708620
Status: VALIDATED
Authors:
  • Deniz Sunar Cerci
  • Sercan Sen
  • Salim Cerci
  • Caglar Zorbilmez
  • Ilknur Hos
References:
  • CMS-FSQ-15-006
  • CERN-EP-2018-308
  • arXiv: 1812.04095
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • Inelastic, single-diffractive and non-single diffractive events at $\sqrt{s} = 13$ TeV.

A measurement of the energy density in proton-proton collisions at a centre-of-mass energy of $\sqrt{s} = 13$ TeV is presented. The data have been recorded with the CMS experiment at the LHC during low luminosity operations in 2015. The energy density is studied as a function of pseudorapidity in the ranges $-6.6 < \eta < -5.2$ and $3.15 < |\eta| < 5.20$. The results are compared with the predictions of several models. All the models considered suggest a different shape of the pseudorapidity dependence compared to that observed in the data. A comparison with LHC proton-proton collision data at $\sqrt{s} = 0.9$ and 7 TeV confirms the compatibility of the data with the hypothesis of limiting fragmentation.

Source code: CMS_2018_I1708620.cc
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
/// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {


  /// Forward energy flow at 13 TeV with CMS
  ///
  /// @note Rivet 3 conversion by A Buckley
  class CMS_2018_I1708620 : public Analysis {
  public:

    /// Constructor
    DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2018_I1708620);


    /// Initialise
    void init() {

      declare(FinalState(), "FS");
      const ChargedFinalState cfsBSCplus(Cuts::eta > 3.9 && Cuts::eta < 4.4);
      declare(cfsBSCplus, "cfsBSCplus");
      const ChargedFinalState cfsBSCminus(Cuts::eta > -4.4 && Cuts::eta < -3.9);
      declare(cfsBSCminus, "cfsBSCminus");

      book(_noe_inel, "TMP/Ninel");
      book(_noe_nsd, "TMP/Nnsd");
      book(_noe_bsc, "TMP/Nbsc");
      book(_noe_sd, "TMP/Nsd");
      book(_noe_nsd_sd, "TMP/Nnsdsd");

      book(_h_inel, 1, 1, 1);
      book(_h_nsd , 2, 1, 1);
      book(_h_et  , 3, 1, 1);
      book(_h_sd  , 4, 1, 1);
    }


    void analyze(const Event& event) {

      const ChargedFinalState& cfsBSCplus = applyProjection<ChargedFinalState>(event, "cfsBSCplus");
      const ChargedFinalState& cfsBSCminus = applyProjection<ChargedFinalState>(event, "cfsBSCminus");
      const bool bscplus = !cfsBSCplus.empty();
      const bool bscminus = !cfsBSCminus.empty();

      // Find final-state particles
      const FinalState& fs = applyProjection<FinalState>(event, "FS");
      // const Particles particlesByRapidity = fs.particlesByPt();
      // sortBy(particlesByRapidity, cmpMomByRap);
      const Particles particlesByRapidity = fs.particles(cmpMomByRap);
      const size_t num_particles = particlesByRapidity.size();

      // Find gaps, and choose the middle one as the event gap centre
      vector<double> gaps, midpoints;
      for (size_t ip = 1; ip < num_particles; ++ip) {
        const Particle& p1 = particlesByRapidity[ip-1];
        const Particle& p2 = particlesByRapidity[ip];
        const double gap = p2.rapidity() - p1.rapidity();
        const double mid = (p2.rapidity() + p1.rapidity()) / 2;
        gaps.push_back(gap);
        midpoints.push_back(mid);
      }
      const size_t imid = std::distance(gaps.begin(), max_element(gaps.begin(), gaps.end()));
      const double gapcenter = midpoints[imid];

      // Assign particles to sides and compute Mx, My
      FourMomentum v4mx, v4my;
      for (const Particle& p : particlesByRapidity) {
        (p.rapidity() > gapcenter ? v4mx : v4my) += p.momentum();
      }
      const double mx = v4mx.mass();
      const double my = v4my.mass();

      // Compute xi variables
      const double xix = sqr(mx/(sqrtS()/GeV));
      const double xiy = sqr(my/(sqrtS()/GeV));
      const double xi_sd = max(xix, xiy);

      // Compute if inelastic, and other variables
      const bool inel = (xi_sd > 1e-6);
      if (inel) _noe_inel->fill();
      const bool bsc = (bscplus && bscminus);
      if (bsc) _noe_bsc->fill(); ///< @todo Not re-entry safe: FIX

      // Count/histogram backward and forward Et
      static const double YBEAM = 9.54;
      int nplus  = 0, nminus = 0;
      for (const Particle& p : particlesByRapidity) {
        const double eta = p.eta();
        if (!p.isVisible()) continue;
        if (inRange(eta, 2.866, 5.205)) nplus += 1;
        if (inRange(eta, -5.205, -2.866)) nminus += 1;
        if (bsc) _h_et->fill(eta-YBEAM, p.Et()/GeV);
      }

      // Categorise as NSD
      const bool nsd = (nminus > 0 && nplus > 0);
      if (nsd) _noe_nsd->fill();

      // Histogram
      for (const Particle& p : particlesByRapidity) {
        if (inel) _h_inel->fill(p.abseta(), p.E()/GeV);
        if (nsd) _h_nsd->fill(p.abseta(), p.E()/GeV);
      }

      // SD selection
      static const double ETAMIN = 3.152;
      static const double ETAMAX = 5.205;
      static const double EMIN   = 5*GeV;
      bool stableParticleEnergyCutMinus = false, stableParticleEnergyCutPlus = false;
      for (const Particle& p : particlesByRapidity) {
        if (p.E() < EMIN) continue;
        if (inRange(p.eta(), -ETAMAX, -ETAMIN)) stableParticleEnergyCutMinus = true;
        if (inRange(p.eta(),  ETAMIN,  ETAMAX)) stableParticleEnergyCutPlus = true;
      }

      // Select SD-enhanced events with the following condition
      const bool sd = (stableParticleEnergyCutPlus != stableParticleEnergyCutMinus); //< bool XOR
      if (sd) {
        _noe_sd->fill();
        for (const Particle& p : particlesByRapidity) {
          if (inRange(p.abseta(), ETAMIN, ETAMAX)) {
            if (stableParticleEnergyCutPlus && p.eta() > 0) _h_sd->fill(p.abseta(), p.E()/GeV);
            if (stableParticleEnergyCutMinus && p.eta() < 0) _h_sd->fill(p.abseta(), p.E()/GeV);
          } else { // CASTOR
            _h_sd->fill(p.abseta(), p.E()/GeV/2);
          }
        }
      }

      // Count how many are NSD and SD
      if (nsd && sd) _noe_nsd_sd->fill();
    }


    void finalize() {
      scale(_h_inel, 0.5/_noe_inel->sumW());
      scale(_h_nsd, 0.5/_noe_nsd->sumW());
      scale(_h_et, 1/_noe_bsc->sumW());
      scale(_h_sd, 1/_noe_sd->sumW());
      MSG_INFO( "Number of events of INEL: " << _noe_inel->effNumEntries() );
      MSG_INFO( "Number of events of NSD: " << _noe_nsd->effNumEntries() );
      MSG_INFO( "Number of events of SD: " << _noe_sd->effNumEntries() );
      MSG_INFO( "Number of events of NSD and SD contribution: " << _noe_nsd_sd->effNumEntries() );
    }


    Histo1DPtr _h_inel, _h_nsd, _h_et, _h_sd;
    CounterPtr _noe_inel, _noe_nsd, _noe_bsc, _noe_sd, _noe_nsd_sd;

  };


  DECLARE_RIVET_PLUGIN(CMS_2018_I1708620);

}