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