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