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    double phaseDif(double a1, double a2, const pair<double, double>& edges) {
 26      double dif = a1 - a2;
 27      while (dif < edges.first)
 28        dif += 2*M_PI;
 29      while (dif > edges.second)
 30        dif -= 2*M_PI;
 31      return dif;
 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      // Resize all the analysis object containers to right size.
 69      ratio.resize(refdata.size());
 70      signal.resize(refdata.size());
 71      background.resize(refdata.size());
 72      nsp.resize(refdata.size());
 73      nmp.resize(refdata.size());
 74      for (int i = 0, N = refdata.size(); i < N; ++i) {
 75        const YODA::Estimate1D& tmp = refData(refdata[i]);
 76        // The ratio plots.
 77        book(ratio[i], refdata[i]);
 78        // Signal and mixed background should not be displayed.
 79        book(signal[i], "TMP/" + refdata[i] + "-s", tmp);
 80        book(background[i], "TMP/" + refdata[i] + "-b", tmp);
 81        // Number of signal and mixed pairs for normalization.
 82        book(nsp[i],"TMP/nsp"+std::to_string(i));
 83        book(nmp[i],"TMP/nmp"+std::to_string(i));
 84        // The differing deltaphi histogram edges per pair.
 85        deltaphi.push_back({tmp.xMin(), tmp.xMax()});
 86      }
 87    }
 88
 89
 90    void fillPair(const Particle& p1, const Particle& p2, vector<Histo1DPtr>& histos,
 91      vector<CounterPtr>& sow) {
 92	   if (isSame(p1,p2)) return;
 93          // If the pair is not within eta acceptance, we can continue early.
 94          if (abs(p1.eta() - p2.eta()) > 1.3) return;
 95          // Figure out which pid pair we are looking at.
 96          int iPair = -1;
 97          for (int i = 0, N = pid.size(); i < N; ++i) {
 98            if (pid[i].first == p1.pid() && pid[i].second == p2.pid()) {
 99              iPair = i;
100              break;
101            }
102          }
103          // If the pair is not in the analysis, don't fill anything.
104          if (iPair < 0) return;
105          // Apply min pT cuts, varies for different species.
106          if (p1.pT() < pTcuts[iPair].first || p2.pT() < pTcuts[iPair].second) return;
107          const double dPhi = phaseDif(p1.phi(), p2.phi(), deltaphi[iPair]);
108          histos[iPair]->fill(dPhi);
109          sow[iPair]->fill();
110    }
111
112
113    /// Perform the per-event analysis
114    void analyze(const Event& event) {
115      // Triggering.
116      if (!apply<ALICE::V0AndTrigger>(event, "V0-AND")()) return;
117
118      // The projections for signal and mixed event background.
119      const PrimaryParticles& pp =
120        apply<PrimaryParticles>(event,"APRIM");
121      const EventMixingFinalState& evm =
122        apply<EventMixingFinalState>(event, "EVM");
123
124      // Test if we have enough mixing events available to continue.
125      if (!evm.hasMixingEvents()) return;
126
127      for (const Particle& p1 : pp.particles()) {
128	      // First do the signal histograms.
129        for (const Particle& p2 : pp.particles())
130	        fillPair(p1, p2, signal, nsp);
131	      // Then do the background
132        for (const Particle& p2 : evm.particles())
133	        fillPair(p1, p2, background, nmp);
134      }
135    }
136
137
138    /// Normalise histograms etc., after the run
139    void finalize() {
140      for (int i = 0, N = pid.size(); i < N; ++i) {
141	    // Scaling factor eqns. (2)-(5) in the paper.
142        double sc = nmp[i]->sumW() / nsp[i]->sumW();
143        signal[i]->scaleW(sc);
144        divide(signal[i],background[i],ratio[i]);
145      }
146    }
147
148    /// @}
149
150
151    /// Analysis variables.
152    vector<pair<int, int> > pid;
153    vector<pair<double, double> > pTcuts;
154    vector<pair<double, double> > deltaphi;
155    /// @name Histograms and counters
156    /// @{
157    vector<Histo1DPtr> signal;
158    vector<Histo1DPtr> background;
159    vector<Estimate1DPtr> ratio;
160    vector<CounterPtr> nsp;
161    vector<CounterPtr> nmp;
162
163    /// @}
164
165  };
166
167
168  RIVET_DECLARE_PLUGIN(ALICE_2016_I1507157);
169
170}