1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Tools/Correlators.hh"
  5#include "Rivet/Analyses/AliceCommon.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Multiparticle azimuthal correlations pp, pPb, XeXe and PbPb.
 11  class ALICE_2019_I1723697 : public CumulantAnalysis {
 12  public:
 13
 14    /// Constructor
 15    ALICE_2019_I1723697() :
 16      CumulantAnalysis("ALICE_2019_I1723697") {}
 17
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      // Find system type.
 26      string bOpt = getOption<string>("beam","");
 27      const ParticlePair& beam = beams();
 28      if (beam.first.pid() == PID::PROTON && beam.second.pid() == PID::PROTON)
 29         sysType = pp;
 30      else if (beam.first.pid() == PID::PROTON && beam.second.pid() == PID::LEAD)
 31        sysType = pPb;
 32      else if (beam.first.pid() == PID::XENON && beam.second.pid() == PID::XENON)
 33        sysType = XeXe;
 34      else if (beam.first.pid() == PID::LEAD && beam.second.pid() == PID::LEAD)
 35        sysType = PbPb;
 36      else{
 37        MSG_WARNING("Suspicious beam. You're probably in reentrant mode. Fetching beam type from option.");
 38        if (bOpt == "pp") sysType = pp;
 39        else if (bOpt == "pPb") sysType = pPb;
 40        else if (bOpt == "XeXe") sysType = XeXe;
 41        else if (bOpt == "PbPb") sysType = PbPb;
 42        else MSG_ERROR("Could not decipher beam type. For rivet-merge, set -a ALICE_2019_I1723697:beam=OPT, where opt is pp, pPb, XeXe or PbPb.");
 43        MSG_WARNING("Setting beam type to "+bOpt);
 44      }
 45      // Sanity check
 46      if ((sysType == pp && bOpt != "pp") || (sysType == pPb && bOpt != "pPb") ||
 47        (sysType == XeXe && bOpt != "XeXe") || (sysType == PbPb && bOpt != "PbPb"))
 48          MSG_WARNING("Beam option and registered beam don't match. Are you sure this is intentional?");
 49
 50      // Initialise and register projections
 51      // Declare the trigger projection.
 52      declare<ALICE::V0AndTrigger>(ALICE::V0AndTrigger(),"V0-AND");
 53
 54      // Centrality projection for high multiplicity trigger in pp.
 55      if (sysType == pp)
 56        declareCentrality(ALICE::V0MMultiplicity(),
 57	  "ALICE_2015_PPCentrality", "V0M","V0M");
 58
 59      // The full central charged final state.
 60      const ChargedFinalState& cfs = ChargedFinalState(Cuts::abseta < 0.8 &&
 61        Cuts::pT > 0.2*GeV && Cuts::pT < 5.0*GeV);
 62      declare(cfs, "CFS");
 63
 64      // The positive eta side used for rapidity gap = 1.4.
 65      const ChargedFinalState& cfsp14 = ChargedFinalState(Cuts::eta > 0.5 &&
 66        Cuts::eta < 0.8 && Cuts::pT > 0.2*GeV && Cuts::pT < 5.0*GeV);
 67      declare(cfsp14, "CFSP14");
 68      // ..negative ditto.
 69      const ChargedFinalState& cfsn14 = ChargedFinalState(Cuts::eta < -0.5 &&
 70        Cuts::eta > -0.8 && Cuts::pT > 0.2*GeV && Cuts::pT < 5.0*GeV);
 71      declare(cfsn14, "CFSN14");
 72      // The positive eta side used for rapidity gap = 1.0.
 73      const ChargedFinalState& cfsp10 = ChargedFinalState(Cuts::eta > 0.5 &&
 74        Cuts::eta < 0.8 && Cuts::pT > 0.2*GeV && Cuts::pT < 5.0*GeV);
 75      declare(cfsp10, "CFSP10");
 76      // ..negative ditto.
 77      const ChargedFinalState& cfsn10 = ChargedFinalState(Cuts::eta < -0.5 &&
 78        Cuts::eta > -0.8 && Cuts::pT > 0.2*GeV && Cuts::pT < 5.0*GeV);
 79      declare(cfsn10, "CFSN10");
 80
 81      // Book flow coeff scatters before booking the correlators
 82      // to have access to bin edges.
 83      int n1 = 0;
 84      int n2 = 0;
 85      int n3 = 0;
 86      if (sysType == pPb) n1 = 9, n2 = 1;
 87      else if (sysType == XeXe) n1 = 19, n2 = 1, n3 = 0;
 88      else if (sysType == PbPb) n1 = 30, n2 = 1, n3 = 1;
 89      book(h_v22gap, 1 + n1, 1, 1, true);
 90      book(h_v32gap, 2 + n1, 1, 1, true);
 91      book(h_v42gap, 3 + n1, 1, 1, true);
 92      if (sysType != pp) {
 93        book(h_v24, 4 + n1, 1, 1, true);
 94      }
 95      book(h_v26, 5 + n1 + n2, 1, 1, true);
 96      if (sysType == XeXe || sysType == PbPb) {
 97        book(h_v28, 6 + n1 + n2 + n3, 1, 1, true);
 98      }
 99
100      // Book cumulant scatters.
101      book(h_c22gap, "c22gap", refData(1 + n1, 1, 1).mkScatter());
102      book(h_c32gap, "c32gap", refData(2 + n1, 1, 1).mkScatter());
103      book(h_c42gap, "c42gap", refData(3 + n1, 1, 1).mkScatter());
104      book(h_c24, "c24", refData(4 + n1, 1, 1).mkScatter());
105      book(h_c26, "c26", refData(5 + n1 + n2, 1, 1).mkScatter());
106      if (sysType == XeXe || sysType == PbPb)
107        book(h_c28, "c28", refData(6 + n1 + n2 + n3, 1, 1).mkScatter());
108
109      // Book correlators. First the ungapped ones.
110      ec22 = bookECorrelator<2,2>("ec22", refData(4 + n1, 1, 1));
111      ec24 = bookECorrelator<2,4>("ec24", refData(4 + n1, 1, 1));
112
113      ec622 = bookECorrelator<2,2>("ec622", refData(5 + n1 + n2, 1, 1));
114      ec624 = bookECorrelator<2,4>("ec624", refData(5 + n1 + n2, 1, 1));
115      ec626 = bookECorrelator<2,6>("ec626", refData(5 + n1 + n2, 1, 1));
116
117      // And the ones just valid for XeXe and PbPb.
118      if (sysType == XeXe || sysType == PbPb) {
119	    ec822 = bookECorrelator<2,2>("ec822", refData(6 + n1 + n2 + n3, 1, 1));
120	    ec824 = bookECorrelator<2,4>("ec824", refData(6 + n1 + n2 + n3, 1, 1));
121	    ec826 = bookECorrelator<2,6>("ec826", refData(6 + n1 + n2 + n3, 1, 1));
122	    ec828 = bookECorrelator<2,8>("ec828", refData(6 + n1 + n2 + n3, 1, 1));
123      }
124      // ...and the gapped ones.
125      ec22gap = bookECorrelatorGap<2,2>("ec22gap", refData(1 + n1, 1, 1));
126      ec32gap = bookECorrelatorGap<3,2>("ec32gap", refData(2 + n1, 1, 1));
127      ec42gap = bookECorrelatorGap<4,2>("ec42gap", refData(3 + n1, 1, 1));
128      // Get the max order of booked correlators for the projections.
129      pair<int, int> max = getMaxValues();
130
131      // Declare correlator projections.
132      declare(Correlators(cfs, max.first, max.second), "Correlators");
133      declare(Correlators(cfsp14, max.first, max.second), "CorrelatorsPos14");
134      declare(Correlators(cfsn14, max.first, max.second), "CorrelatorsNeg14");
135      declare(Correlators(cfsp10, max.first, max.second), "CorrelatorsPos10");
136      declare(Correlators(cfsn10, max.first, max.second), "CorrelatorsNeg10");
137    }
138
139
140    /// Perform the per-event analysis
141    void analyze(const Event& event) {
142      // Event trigger.
143      if (!apply<ALICE::V0AndTrigger>(event, "V0-AND")() ) vetoEvent;
144
145      // High multiplicity trigger for pp.
146      if (sysType == pp && apply<CentralityProjection>(event, "V0M")() > 0.1)
147        vetoEvent;
148      // Number of charged particles.
149      const double nch = apply<ChargedFinalState>(event, "CFS").size();
150
151      // The correlators projections.
152      const Correlators& c = apply<Correlators>(event,"Correlators");
153      const Correlators& cp10 =
154        apply<Correlators>(event,"CorrelatorsPos10");
155      const Correlators& cn10 =
156        apply<Correlators>(event,"CorrelatorsNeg10");
157      const Correlators& cp14 =
158        apply<Correlators>(event,"CorrelatorsPos14");
159      const Correlators& cn14 =
160        apply<Correlators>(event,"CorrelatorsNeg14");
161
162      // Fill correlators.
163      ec22->fill(nch, c);
164      ec24->fill(nch, c);
165
166      ec622->fill(nch, c);
167      ec624->fill(nch, c);
168      ec626->fill(nch, c);
169
170      if (sysType == XeXe || sysType == PbPb) {
171        ec822->fill(nch, c);
172        ec824->fill(nch, c);
173        ec826->fill(nch, c);
174        ec828->fill(nch, c);
175      }
176
177      // Fill gapped correlators.
178      ec22gap->fill(nch, cp14, cn14);
179      ec32gap->fill(nch, cp10, cn10);
180      ec42gap->fill(nch, cp10, cn10);
181
182    }
183
184    /// Normalise histograms etc., after the run
185    void finalize() {
186      // Stream correlators to yoda file.
187      // Fill cumulant scatters.
188      cnTwoInt(h_c22gap, ec22gap);
189      cnTwoInt(h_c32gap, ec32gap);
190      cnTwoInt(h_c42gap, ec42gap);
191
192      cnFourInt(h_c24, ec22, ec24);
193      cnSixInt(h_c26, ec622, ec624, ec626);
194      if (sysType == XeXe || sysType == PbPb)
195        cnEightInt(h_c28, ec822, ec824, ec826, ec828);
196
197      // Fill flow scatters.
198      vnTwoInt(h_v22gap, ec22gap);
199      vnTwoInt(h_v32gap, ec32gap);
200      vnTwoInt(h_v42gap, ec42gap);
201
202      if (sysType != pp)
203        vnFourInt(h_v24, ec22, ec24);
204      vnSixInt(h_v26, ec622, ec624, ec626);
205      if (sysType == XeXe || sysType == PbPb)
206        vnEightInt(h_v28, ec822, ec824, ec826, ec828);
207    }
208
209    /// @}
210    // System check enum.
211    enum SysType {pp, pPb, XeXe, PbPb};
212    SysType sysType;
213
214    /// @name Histograms
215    /// @{
216    // Flow coefficients.
217    Scatter2DPtr h_v22gap;
218    Scatter2DPtr h_v32gap;
219    Scatter2DPtr h_v42gap;
220    Scatter2DPtr h_v24;
221    Scatter2DPtr h_v26;
222    Scatter2DPtr h_v28;
223
224    // Cumulants.
225    Scatter2DPtr h_c22gap;
226    Scatter2DPtr h_c32gap;
227    Scatter2DPtr h_c42gap;
228    Scatter2DPtr h_c24;
229    Scatter2DPtr h_c26;
230    Scatter2DPtr h_c28;
231
232    // Correlators.
233    ECorrPtr ec22;
234    ECorrPtr ec24;
235
236    ECorrPtr ec622;
237    ECorrPtr ec624;
238    ECorrPtr ec626;
239
240    ECorrPtr ec822;
241    ECorrPtr ec824;
242    ECorrPtr ec826;
243    ECorrPtr ec828;
244
245    // Gapped correlators.
246    ECorrPtr ec22gap;
247    ECorrPtr ec32gap;
248    ECorrPtr ec42gap;
249
250    /// @}
251
252
253  };
254
255
256  RIVET_DECLARE_PLUGIN(ALICE_2019_I1723697);
257
258
259}