1/// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// Forward energy flow at 13 TeV with CMS
  9  ///
 10  /// @note Rivet 3 conversion by A Buckley
 11  class CMS_2018_I1708620 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2018_I1708620);
 16
 17
 18    /// Initialise
 19    void init() {
 20
 21      declare(FinalState(), "FS");
 22      const ChargedFinalState cfsBSCplus(Cuts::eta > 3.9 && Cuts::eta < 4.4);
 23      declare(cfsBSCplus, "cfsBSCplus");
 24      const ChargedFinalState cfsBSCminus(Cuts::eta > -4.4 && Cuts::eta < -3.9);
 25      declare(cfsBSCminus, "cfsBSCminus");
 26
 27      book(_noe_inel, "TMP/Ninel");
 28      book(_noe_nsd, "TMP/Nnsd");
 29      book(_noe_bsc, "TMP/Nbsc");
 30      book(_noe_sd, "TMP/Nsd");
 31      book(_noe_nsd_sd, "TMP/Nnsdsd");
 32
 33      book(_h_inel, 1, 1, 1);
 34      book(_h_nsd , 2, 1, 1);
 35      book(_h_et  , 3, 1, 1);
 36      book(_h_sd  , 4, 1, 1);
 37    }
 38
 39
 40    void analyze(const Event& event) {
 41
 42      const ChargedFinalState& cfsBSCplus = apply<ChargedFinalState>(event, "cfsBSCplus");
 43      const ChargedFinalState& cfsBSCminus = apply<ChargedFinalState>(event, "cfsBSCminus");
 44      const bool bscplus = !cfsBSCplus.empty();
 45      const bool bscminus = !cfsBSCminus.empty();
 46
 47      // Find final-state particles
 48      const FinalState& fs = apply<FinalState>(event, "FS");
 49      // const Particles particlesByRapidity = fs.particlesByPt();
 50      // sortBy(particlesByRapidity, cmpMomByRap);
 51      const Particles particlesByRapidity = fs.particles(cmpMomByRap);
 52      const size_t num_particles = particlesByRapidity.size();
 53
 54      // Find gaps, and choose the middle one as the event gap centre
 55      vector<double> gaps, midpoints;
 56      for (size_t ip = 1; ip < num_particles; ++ip) {
 57        const Particle& p1 = particlesByRapidity[ip-1];
 58        const Particle& p2 = particlesByRapidity[ip];
 59        const double gap = p2.rapidity() - p1.rapidity();
 60        const double mid = (p2.rapidity() + p1.rapidity()) / 2;
 61        gaps.push_back(gap);
 62        midpoints.push_back(mid);
 63      }
 64      const size_t imid = std::distance(gaps.begin(), max_element(gaps.begin(), gaps.end()));
 65      const double gapcenter = midpoints[imid];
 66
 67      // Assign particles to sides and compute Mx, My
 68      FourMomentum v4mx, v4my;
 69      for (const Particle& p : particlesByRapidity) {
 70        (p.rapidity() > gapcenter ? v4mx : v4my) += p.momentum();
 71      }
 72      const double mx = v4mx.mass();
 73      const double my = v4my.mass();
 74
 75      // Compute xi variables
 76      const double xix = sqr(mx/(sqrtS()/GeV));
 77      const double xiy = sqr(my/(sqrtS()/GeV));
 78      const double xi_sd = max(xix, xiy);
 79
 80      // Compute if inelastic, and other variables
 81      const bool inel = (xi_sd > 1e-6);
 82      if (inel) _noe_inel->fill();
 83      const bool bsc = (bscplus && bscminus);
 84      if (bsc) _noe_bsc->fill(); ///< @todo Not re-entry safe: FIX
 85
 86      // Count/histogram backward and forward Et
 87      static const double YBEAM = 9.54;
 88      int nplus  = 0, nminus = 0;
 89      for (const Particle& p : particlesByRapidity) {
 90        const double eta = p.eta();
 91        if (!p.isVisible()) continue;
 92        if (inRange(eta, 2.866, 5.205)) nplus += 1;
 93        if (inRange(eta, -5.205, -2.866)) nminus += 1;
 94        if (bsc) _h_et->fill(eta-YBEAM, p.Et()/GeV);
 95      }
 96
 97      // Categorise as NSD
 98      const bool nsd = (nminus > 0 && nplus > 0);
 99      if (nsd) _noe_nsd->fill();
100
101      // Histogram
102      for (const Particle& p : particlesByRapidity) {
103        if (inel) _h_inel->fill(p.abseta(), p.E()/GeV);
104        if (nsd) _h_nsd->fill(p.abseta(), p.E()/GeV);
105      }
106
107      // SD selection
108      static const double ETAMIN = 3.152;
109      static const double ETAMAX = 5.205;
110      static const double EMIN   = 5*GeV;
111      bool stableParticleEnergyCutMinus = false, stableParticleEnergyCutPlus = false;
112      for (const Particle& p : particlesByRapidity) {
113        if (p.E() < EMIN) continue;
114        if (inRange(p.eta(), -ETAMAX, -ETAMIN)) stableParticleEnergyCutMinus = true;
115        if (inRange(p.eta(),  ETAMIN,  ETAMAX)) stableParticleEnergyCutPlus = true;
116      }
117
118      // Select SD-enhanced events with the following condition
119      const bool sd = (stableParticleEnergyCutPlus != stableParticleEnergyCutMinus); //< bool XOR
120      if (sd) {
121        _noe_sd->fill();
122        for (const Particle& p : particlesByRapidity) {
123          if (inRange(p.abseta(), ETAMIN, ETAMAX)) {
124            if (stableParticleEnergyCutPlus && p.eta() > 0) _h_sd->fill(p.abseta(), p.E()/GeV);
125            if (stableParticleEnergyCutMinus && p.eta() < 0) _h_sd->fill(p.abseta(), p.E()/GeV);
126          } else { // CASTOR
127            _h_sd->fill(p.abseta(), p.E()/GeV/2);
128          }
129        }
130      }
131
132      // Count how many are NSD and SD
133      if (nsd && sd) _noe_nsd_sd->fill();
134    }
135
136
137    void finalize() {
138      scale(_h_inel, 0.5/_noe_inel->sumW());
139      scale(_h_nsd, 0.5/_noe_nsd->sumW());
140      scale(_h_et, 1/_noe_bsc->sumW());
141      scale(_h_sd, 1/_noe_sd->sumW());
142      MSG_INFO( "Number of events of INEL: " << _noe_inel->effNumEntries() );
143      MSG_INFO( "Number of events of NSD: " << _noe_nsd->effNumEntries() );
144      MSG_INFO( "Number of events of SD: " << _noe_sd->effNumEntries() );
145      MSG_INFO( "Number of events of NSD and SD contribution: " << _noe_nsd_sd->effNumEntries() );
146    }
147
148
149    Histo1DPtr _h_inel, _h_nsd, _h_et, _h_sd;
150    CounterPtr _noe_inel, _noe_nsd, _noe_bsc, _noe_sd, _noe_nsd_sd;
151
152  };
153
154
155  RIVET_DECLARE_PLUGIN(CMS_2018_I1708620);
156
157}