1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  class ALICE_2012_I1181770 : public Analysis {
  9  public:
 10
 11    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2012_I1181770);
 12
 13
 14    void init() {
 15      // Projection setup
 16      declare(ChargedFinalState(), "CFS");
 17
 18      // Book histograms
 19      book(_h_frac_sd_inel, 1, 1, 1);
 20      book(_h_frac_dd_inel, 2, 1, 1);
 21      book(_h_xsec_sd     , 3, 1, 1);
 22      book(_h_xsec_dd     , 4, 1, 1);
 23      book(_h_xsec_inel   , 5, 1, 1);
 24    }
 25
 26
 27    void analyze(const Event& event) {
 28      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 29      if (cfs.size() < 2) vetoEvent; // need at least two particles to calculate gaps
 30
 31      // Fill INEL plots for each event
 32      _h_xsec_inel->fill(sqrtS()/GeV);
 33
 34      // Identify particles with most positive/most negative rapidities
 35      const Particles particlesByRap = cfs.particles(cmpMomByRap);
 36      const Particle pslowest = particlesByRap.front();
 37      const Particle pfastest = particlesByRap.back();
 38
 39      // Find gap sizes
 40      const Particles particlesByEta = cfs.particles(cmpMomByEta); // sorted from minus to plus
 41      const size_t num_particles = particlesByEta.size();
 42      vector<double> gaps;
 43      for (size_t ip = 1; ip < num_particles; ++ip) {
 44        const Particle& p1 = particlesByEta[ip-1];
 45        const Particle& p2 = particlesByEta[ip];
 46        const double gap = p2.eta() - p1.eta();
 47        assert(gap >= 0);
 48        gaps.push_back(gap);
 49      }
 50
 51      // First, last, and largest gaps
 52      const double gapmax = *max_element(gaps.begin(), gaps.end());
 53      const double gapbwd = gaps.front();
 54      const double gapfwd = gaps.back();
 55
 56      // Mx calculation
 57      FourMomentum p4lead;
 58      if (pslowest.pid() == PID::PROTON && pfastest.pid() == PID::PROTON) {
 59        p4lead = (fabs(pslowest.rapidity()) > fabs(pfastest.rapidity())) ? pslowest.momentum() : pfastest.momentum();
 60      } else if (pslowest.pid() == PID::PROTON) {
 61        p4lead = pslowest.momentum();
 62      } else if (pfastest.pid() == PID::PROTON) {
 63        p4lead = pfastest.momentum();
 64      }
 65      const double Mx = sqrt( (sqrtS()-p4lead.E()-p4lead.p3().mod()) * (sqrtS()-p4lead.E()+p4lead.p3().mod()) );
 66
 67      // Fill SD (and escape) if Mx is sufficiently low
 68      if (Mx < 200*GeV) {
 69        _h_xsec_sd->fill(sqrtS()/GeV);
 70        return;
 71      }
 72
 73      // Also remove SD-like events in NSD events
 74      if (fuzzyEquals(gapbwd, gapmax) || fuzzyEquals(gapfwd, gapmax)) vetoEvent;
 75
 76      // Fill DD plots
 77      if (gapmax > 3) _h_xsec_dd->fill(sqrtS()/GeV);
 78    }
 79
 80
 81    void finalize() {
 82
 83      // get the ratio plots: SD/inel, DD/inel
 84      divide(_h_xsec_sd , _h_xsec_inel, _h_frac_sd_inel);
 85      divide(_h_xsec_sd , _h_xsec_inel, _h_frac_dd_inel);
 86
 87      const double scaling = crossSection()/millibarn/sumOfWeights();
 88      scale(_h_xsec_sd,   scaling);
 89      scale(_h_xsec_dd,   scaling);
 90      scale(_h_xsec_inel, scaling);
 91
 92    }
 93
 94  private:
 95
 96    Estimate1DPtr _h_frac_sd_inel;
 97    Estimate1DPtr _h_frac_dd_inel;
 98    Histo1DPtr   _h_xsec_sd;
 99    Histo1DPtr   _h_xsec_dd;
100    Histo1DPtr   _h_xsec_inel;
101
102  };
103
104
105  RIVET_DECLARE_PLUGIN(ALICE_2012_I1181770);
106
107}