1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Analyses/AliceCommon.hh"
  4#include "Rivet/Projections/PrimaryParticles.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/EventMixingFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Angular correlations of identified particles in pp at 7 TeV.
 12  ///
 13  /// Also showcasing use of EventMixingFinalState.
 14  class ALICE_2016_I1507157 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2016_I1507157);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// @brief Calculate angular distance between particles,
 25    /// return index of corresponding discrete bin
 26    size_t discPhi(double a1, double a2, const YODA::Axis<double>& axis) const {
 27      double diff = a1 - a2;
 28      while (diff < axis.min(1))               diff += 2*M_PI;
 29      while (diff > axis.max(axis.numBins()))  diff -= 2*M_PI;
 30      size_t idx = axis.index(diff);
 31      return idx > axis.numBins()? 0 : idx;
 32    }
 33
 34
 35    /// Book histograms and initialise projections before the run
 36    void init() {
 37
 38      const double etamax = 0.8;
 39      const double pTmin = 0.2; // GeV
 40      const double pTmax = 2.5; //GeV
 41
 42      // Trigger projection.
 43      declare(ALICE::V0AndTrigger(), "V0-AND");
 44      // Charged tracks used to manage the mixing observable.
 45      ChargedFinalState cfsMult(Cuts::abseta < etamax);
 46      declare(cfsMult, "CFSMult");
 47
 48      // Primary particles.
 49      PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS,
 50	      Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON,
 51	      Rivet::PID::NEUTRON, Rivet::PID::LAMBDA, Rivet::PID::SIGMAMINUS,
 52       	Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0,
 53	      Rivet::PID::OMEGAMINUS},Cuts::abseta < etamax && Cuts::pT > pTmin*GeV && Cuts::pT < pTmax*GeV);
 54      declare(pp,"APRIM");
 55
 56      // The event mixing projection
 57      declare(EventMixingFinalState(cfsMult, pp, 5, 0, 100, 10, defaultWeightIndex()),"EVM");
 58      // The particle pairs.
 59      pid = {{211, -211}, {321, -321}, {2212, -2212}, {3122, -3122}, {211, 211},
 60             {321, 321}, {2212, 2212}, {3122, 3122}, {2212, 3122}, {2212, -3122}};
 61      // The differing pT cuts per pair, in GeV.
 62      pTcuts = {{0.2, 0.2},{0.3, 0.3},{0.5,0.5},{0.6,0.6},{0.2,0.2},
 63	      {0.3,0.3},{0.5,0.5},{0.6,0.6},{0.5,0.6},{0.5,0.6}};
 64      // The associated histograms in the data file.
 65      vector<string> refdata = {"d04-x01-y01","d04-x01-y02","d04-x01-y03",
 66        "d06-x01-y02","d05-x01-y01","d05-x01-y02","d05-x01-y03","d06-x01-y01",
 67        "d01-x01-y02","d02-x01-y02"};
 68      _axes.resize(3);
 69      _axes[0] = YODA::Axis<double>({-1.46, -1.19, -0.97, -0.76, -0.54, -0.33, -0.11, 0.11, 0.33, 0.54,
 70                                      0.76,  0.98,  1.19,  1.41,  1.62,  1.84,  2.06, 2.27, 2.49, 2.71,
 71                                      2.92,  3.14,  3.36,  3.58,  3.79,  4.01,  4.22, 4.44, 4.66, 4.77});
 72      _axes[1] = YODA::Axis<double>({-1.46, -1.19, -0.97, -0.76, -0.54, -0.33, -0.11, 0.11, 0.33, 0.54,
 73                                      0.76, 0.98, 1.19, 1.41, 1.62, 1.84, 2.06, 2.27, 2.49, 2.71, 2.92,
 74                                      3.14, 3.36, 3.58, 3.79, 4.01, 4.22, 4.44, 4.66, 4.77});
 75      _axes[2] = YODA::Axis<double>({-1.14, -0.95, -0.69, -0.41, -0.14, 0.14, 0.41, 0.68, 0.96, 1.23,
 76                                      1.50,  1.77,  2.05,  2.33,  2.59, 2.87, 3.14, 3.42, 3.68, 3.96,
 77                                      4.24, 4.50, 4.78, 4.92});
 78
 79      // Resize all the analysis object containers to right size.
 80      ratio.resize(refdata.size());
 81      signal.resize(refdata.size());
 82      background.resize(refdata.size());
 83      nsp.resize(refdata.size());
 84      nmp.resize(refdata.size());
 85      for (int i = 0, N = refdata.size(); i < N; ++i) {
 86        const auto& tmp = refData<YODA::BinnedEstimate<string>>(refdata[i]);
 87        // The ratio plots.
 88        book(ratio[i], refdata[i]);
 89        // Signal and mixed background should not be displayed.
 90        book(signal[i], "TMP/" + refdata[i] + "-s", tmp);
 91        book(background[i], "TMP/" + refdata[i] + "-b", tmp);
 92        // Number of signal and mixed pairs for normalization.
 93        book(nsp[i],"TMP/nsp"+std::to_string(i));
 94        book(nmp[i],"TMP/nmp"+std::to_string(i));
 95      }
 96    }
 97
 98
 99    void fillPair(const Particle& p1, const Particle& p2, vector<BinnedHistoPtr<string>>& histos, vector<CounterPtr>& sow) {
100	    if (isSame(p1,p2)) return;
101      // If the pair is not within eta acceptance, we can continue early.
102      if (abs(p1.eta() - p2.eta()) > 1.3) return;
103      // Figure out which pid pair we are looking at.
104      int iPair = -1;
105      for (int i = 0, N = pid.size(); i < N; ++i) {
106        if (pid[i].first == p1.pid() && pid[i].second == p2.pid()) {
107          iPair = i;
108          break;
109        }
110      }
111      // If the pair is not in the analysis, don't fill anything.
112      if (iPair < 0) return;
113      // Apply min pT cuts, varies for different species.
114      if (p1.pT() < pTcuts[iPair].first || p2.pT() < pTcuts[iPair].second) return;
115      size_t axis_id = iPair < 3? 0 : 2;
116      if (4 <= iPair && iPair <= 6)  axis_id = 1;
117      size_t idx = discPhi(p1.phi(), p2.phi(), _axes[axis_id]);
118      histos[iPair]->fill(idx? histos[iPair]->bin(idx).xEdge() : "OTHER"s);
119      sow[iPair]->fill();
120    }
121
122
123    /// Perform the per-event analysis
124    void analyze(const Event& event) {
125      // Triggering.
126      if (!apply<ALICE::V0AndTrigger>(event, "V0-AND")()) return;
127
128      // The projections for signal and mixed event background.
129      const PrimaryParticles& pp =
130        apply<PrimaryParticles>(event,"APRIM");
131      const EventMixingFinalState& evm =
132        apply<EventMixingFinalState>(event, "EVM");
133
134      // Test if we have enough mixing events available to continue.
135      if (!evm.hasMixingEvents()) return;
136
137      for (const Particle& p1 : pp.particles()) {
138	      // First do the signal histograms.
139        for (const Particle& p2 : pp.particles())
140	        fillPair(p1, p2, signal, nsp);
141	      // Then do the background
142        for (const Particle& p2 : evm.particles())
143	        fillPair(p1, p2, background, nmp);
144      }
145    }
146
147
148    /// Normalise histograms etc., after the run
149    void finalize() {
150      for (int i = 0, N = pid.size(); i < N; ++i) {
151	    // Scaling factor eqns. (2)-(5) in the paper.
152        double sc = nmp[i]->sumW() / nsp[i]->sumW();
153        signal[i]->scaleW(sc);
154        divide(signal[i],background[i],ratio[i]);
155      }
156    }
157
158    /// @}
159
160
161    /// Analysis variables.
162    vector<pair<int, int> > pid;
163    vector<pair<double, double> > pTcuts;
164    /// @name Histograms and counters
165    /// @{
166    vector<BinnedHistoPtr<string>> signal;
167    vector<BinnedHistoPtr<string>> background;
168    vector<BinnedEstimatePtr<string>> ratio;
169    vector<CounterPtr> nsp;
170    vector<CounterPtr> nmp;
171    vector<YODA::Axis<double>> _axes;
172
173    /// @}
174
175  };
176
177
178  RIVET_DECLARE_PLUGIN(ALICE_2016_I1507157);
179
180}