1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/SingleValueProjection.hh"
  6#include "Rivet/Analyses/AliceCommon.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// ALICE PbPb at 2.76 TeV azimuthal di-hadron correlations
 12  class ALICE_2012_I930312 : public Analysis {
 13  public:
 14
 15    // Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2012_I930312);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Declare centrality projection
 26      declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_CENT_PBPB", "V0M", "V0M");
 27
 28      // Projection for trigger particles: charged, primary particles
 29      // with |eta| < 1.0 and 8 < pT < 15 GeV/c
 30      declare(ALICE::PrimaryParticles(Cuts::abseta < 1.0 && Cuts::abscharge > 0
 31        && Cuts::ptIn(8.*GeV, 15.*GeV)), "APRIMTrig");
 32
 33      // pT bins edges
 34      vector<double> ptBins = { 3., 4., 6., 8., 10. };
 35
 36      // Projections for associated particles: charged, primary particles
 37      // with |eta| < 1.0 and different pT bins
 38      for (int ipt = 0; ipt < PT_BINS; ++ipt) {
 39        Cut cut = Cuts::abseta < 1.0 && Cuts::abscharge > 0 &&
 40          Cuts::ptIn(ptBins[ipt]*GeV, ptBins[ipt+1]*GeV);
 41        declare(ALICE::PrimaryParticles(cut), "APRIMAssoc" + toString(ipt));
 42      }
 43
 44      // Create event strings
 45      vector<string> evString = { "pp", "central", "peripheral" };
 46
 47      // Initialize trigger counters and yield histograms
 48      string title = "Per trigger particle yield";
 49      string xtitle = "$\\Delta\\eta$ (rad)";
 50      string ytitle = "$1 / N_{trig} {\\rm d}N_{assoc} / {\\rm d}\\Delta\\eta$ (rad$^-1$)";
 51      string hYieldName[EVENT_TYPES][PT_BINS];
 52      for (int itype = 0; itype < EVENT_TYPES; ++itype) {
 53        book(_counterTrigger[itype], "counter." + toString(itype));
 54        for (int ipt = 0; ipt < PT_BINS; ++ipt) {
 55          hYieldName[itype][ipt]= "yield." + evString[itype] + ".pt" + toString(ipt);
 56          book(_histYield[itype][ipt], hYieldName[itype][ipt], 36, -0.5*M_PI, 1.5*M_PI);
 57        }
 58      }
 59
 60      // Find out the beam type, also specified from option.
 61      string beamOpt = getOption<string>("beam","NONE");
 62      if (beamOpt != "NONE") {
 63        MSG_WARNING("You are using a specified beam type, instead of using what"
 64                    "is provided by the generator. "
 65                    "Only do this if you are completely sure what you are doing.");
 66        if (beamOpt=="PP") isHI = false;
 67        else if (beamOpt=="HI") isHI = true;
 68        else {
 69          MSG_ERROR("Beam option error. You have specified an unsupported beam.");
 70          return;
 71      	}
 72      }
 73      else {
 74        const ParticlePair& beam = beams();
 75        if (beam.first.pid() == PID::PROTON && beam.second.pid() == PID::PROTON) isHI = false;
 76        else if (beam.first.pid() == PID::LEAD && beam.second.pid() == PID::LEAD)
 77          isHI = true;
 78        else {
 79          MSG_WARNING("Beam unspecified. Assuming you are running rivet-merge.");
 80        }
 81      }
 82
 83      // Initialize IAA and ICP histograms
 84      book(_histIAA[0], 1, 1, 1);
 85      book(_histIAA[1], 2, 1, 1);
 86      book(_histIAA[2], 5, 1, 1);
 87      book(_histIAA[3], 3, 1, 1);
 88      book(_histIAA[4], 4, 1, 1);
 89      book(_histIAA[5], 6, 1, 1);
 90
 91      // Initialize background-subtracted yield histograms
 92      for (int itype = 0; itype < EVENT_TYPES; ++itype) {
 93        for (int ipt = 0; ipt < PT_BINS; ++ipt) {
 94          book(_histYieldNoBkg[itype][ipt], hYieldName[itype][ipt] + ".nobkg", 36, -0.5*M_PI, 1.5*M_PI);
 95        }
 96      }
 97
 98    }
 99
100
101    /// Perform the per-event analysis
102    void analyze(const Event& event) {
103
104      // Trigger particles
105      Particles trigParticles =
106        apply<ALICE::PrimaryParticles>(event,"APRIMTrig").particles();
107
108      // Associated particles
109      Particles assocParticles[PT_BINS];
110      for (int ipt = 0; ipt < PT_BINS; ++ipt) {
111        string pname = "APRIMAssoc" + toString(ipt);
112        assocParticles[ipt] =
113          apply<ALICE::PrimaryParticles>(event,pname).particles();
114      }
115
116      // Check type of event. This may not be a perfect way to check for the
117      // type of event as there might be some weird conditions hidden inside.
118      // For example some HepMC versions check if number of hard collisions
119      // is equal to 0 and assign 'false' in that case, which is usually wrong.
120      // This might be changed in the future
121      int ev_type = 0; // pp
122      if ( isHI ) {
123        // Prepare centrality projection and value
124        const CentralityProjection& centrProj =
125          apply<CentralityProjection>(event, "V0M");
126        double centr = centrProj();
127        // Set the flag for the type of the event
128        if (centr > 0.0 && centr < 5.0)
129          ev_type = 1; // PbPb, central
130        else if (centr > 60.0 && centr < 90.0)
131          ev_type = 2; // PbPb, peripherial
132        else
133          vetoEvent; // PbPb, other, this is not used in the analysis at all
134      }
135
136      // Fill trigger histogram for a proper event type
137      _counterTrigger[ev_type]->fill(trigParticles.size());
138
139      // Loop over trigger particles
140      for (const Particle& trigParticle : trigParticles) {
141        // For each pt bin
142        for (int ipt = 0; ipt < PT_BINS; ++ipt) {
143          // Loop over associated particles
144          for (const Particle& assocParticle : assocParticles[ipt]) {
145            // If associated and trigger particle are not the same particles.
146            if (!isSame(trigParticle, assocParticle)) {
147              // Test trigger particle.
148              if (trigParticle.pt() > assocParticle.pt()) {
149                // Calculate delta phi in range (-0.5*PI, 1.5*PI).
150                double dPhi = deltaPhi(trigParticle, assocParticle, true);
151                if (dPhi < -0.5 * M_PI) dPhi += 2 * M_PI;
152                // Fill yield histogram for calculated delta phi
153                _histYield[ev_type][ipt]->fill(dPhi);
154              }
155            }
156          }
157        }
158      }
159    }
160
161
162    /// Normalise histograms etc., after the run
163    void finalize() {
164
165      // Check for the reentrant finalize
166      bool pp_available = false, PbPb_available = false;
167      for (int itype = 0; itype < EVENT_TYPES; ++itype) {
168        for (int ipt = 0; ipt < PT_BINS; ++ipt) {
169          if (_histYield[itype][ipt]->numEntries() > 0)
170            itype == 0 ? pp_available = true : PbPb_available = true;
171        }
172      }
173      // Skip postprocessing if pp or PbPb histograms are not available
174      if (!(pp_available && PbPb_available))
175        return;
176
177      // Variable for near and away side peak integral calculation
178      double integral[EVENT_TYPES][PT_BINS][2] = { { {0.0} } };
179
180      // Variables for background calculation
181      double bkg = 0.0;
182      double bkgErr[EVENT_TYPES][PT_BINS] = { {0.0} };
183
184      // Variables for integration error calculation
185      double norm[EVENT_TYPES] = {0.0};
186      double numEntries[EVENT_TYPES][PT_BINS][2] = { { {0.0} } };
187      int numBins[EVENT_TYPES][PT_BINS][2] = { { {0} } };
188
189      // For each event type
190      for (int itype = 0; itype < EVENT_TYPES; ++itype) {
191        // Get counter
192        CounterPtr counter = _counterTrigger[itype];
193        // For each pT range
194        for (int ipt = 0; ipt < PT_BINS; ++ipt) {
195
196          // Get yield histograms
197          Histo1DPtr hYield = _histYield[itype][ipt];
198          Histo1DPtr hYieldNoBkg = _histYieldNoBkg[itype][ipt];
199
200          // Check if histograms are fine
201          if (counter->sumW() == 0 || hYield->numEntries() == 0) {
202            MSG_WARNING("There are no entries in one of the histograms");
203            continue;
204          }
205
206          // Scale yield histogram
207          norm[itype] = 1. / counter->sumW();
208          scale(hYield, norm[itype]);
209
210          // Calculate background
211          double sum = 0.0;
212          int nbins = 0;
213          for (size_t ibin = 0; ibin < hYield->numBins(); ++ibin) {
214            double xmid = hYield->bin(ibin).xMid();
215            if (inRange(xmid, -0.5 * M_PI, -0.5 * M_PI + 0.4) ||
216                inRange(xmid, 0.5 * M_PI - 0.4, 0.5 * M_PI + 0.4) ||
217                inRange(xmid, 1.5 * M_PI - 0.4, 1.5 * M_PI)) {
218              sum += hYield->bin(ibin).sumW();
219              nbins += 1;
220            }
221          }
222          if (nbins == 0) {
223            MSG_WARNING("Failed to estimate background!");
224            continue;
225          }
226          bkg = sum / nbins;
227
228          // Calculate background error
229          sum = 0.0;
230          nbins = 0;
231          for (size_t ibin = 0; ibin < hYield->numBins(); ++ibin) {
232            double xmid = hYield->bin(ibin).xMid();
233            if (inRange(xmid, 0.5 * M_PI - 0.4, 0.5 * M_PI + 0.4)) {
234              sum += (hYield->bin(ibin).sumW() - bkg) *
235                     (hYield->bin(ibin).sumW() - bkg);
236              nbins++;
237            }
238          }
239          if (nbins < 2) {
240            MSG_WARNING("Failed to estimate background error!");
241            continue;
242          }
243          bkgErr[itype][ipt] = sqrt(sum / (nbins - 1));
244
245          // Fill histograms with removed background
246          for (const auto& b : hYield->bins()) {
247            hYieldNoBkg->fill(b.xMid(), b.sumW() - bkg);
248          }
249
250          // Integrate near-side yield
251          size_t lowerBin = hYield->indexAt(-0.7 + 0.02);
252          size_t upperBin = hYield->indexAt( 0.7 - 0.02) + 1;
253          nbins = upperBin - lowerBin;
254          numBins[itype][ipt][NEAR] = nbins;
255          integral[itype][ipt][NEAR] =
256            hYield->integralRange(lowerBin, upperBin) - nbins * bkg;
257          numEntries[itype][ipt][NEAR] =
258            hYield->integralRange(lowerBin, upperBin) * counter->sumW();
259
260          // Integrate away-side yield
261          lowerBin = hYield->indexAt(M_PI - 0.7 + 0.02);
262          upperBin = hYield->indexAt(M_PI + 0.7 - 0.02) + 1;
263          nbins = upperBin - lowerBin;
264          numBins[itype][ipt][AWAY] = nbins;
265          integral[itype][ipt][AWAY] =
266            hYield->integralRange(lowerBin, upperBin) - nbins * bkg;
267          numEntries[itype][ipt][AWAY] =
268            hYield->integralRange(lowerBin, upperBin) * counter->sumW();
269
270        }
271      }
272
273      // Variables for IAA/ICP plots
274      double yval[2] = { 0.0, 0.0 };
275      double yerr[2] = { 0.0, 0.0 };
276
277      int types1[3] = {1, 2, 1};
278      int types2[3] = {0, 0, 2};
279
280      // Fill IAA/ICP plots for near and away side peak
281      for (int ihist = 0; ihist < 3; ++ihist) {
282        int type1 = types1[ihist];
283        int type2 = types2[ihist];
284        double norm1 = norm[type1];
285        double norm2 = norm[type2];
286        for (int ipt = 0; ipt < PT_BINS; ++ipt) {
287          double bkgErr1 = bkgErr[type1][ipt];
288          double bkgErr2 = bkgErr[type2][ipt];
289          for (int ina = 0; ina < 2; ++ina) {
290            double integ1 = integral[type1][ipt][ina];
291            double integ2 = integral[type2][ipt][ina];
292            double numEntries1 = numEntries[type1][ipt][ina];
293            double numEntries2 = numEntries[type2][ipt][ina];
294            double numBins1 = numBins[type1][ipt][ina];
295            double numBins2 = numBins[type2][ipt][ina];
296            yval[ina] = integ1 / integ2;
297            yerr[ina] = norm1 * norm1 * numEntries1 +
298              norm2 * norm2 * numEntries2 * integ1 * integ1 / (integ2 * integ2) +
299              numBins1 * numBins1 * bkgErr1 * bkgErr1 +
300              numBins2 * numBins2 * bkgErr2 * bkgErr2 * integ1 * integ1 / (integ2 * integ2);
301            yerr[ina] = sqrt(yerr[ina])/integ2;
302          }
303          _histIAA[ihist]->bin(ipt+1).set(yval[NEAR], yerr[NEAR]);
304          _histIAA[ihist + 3]->bin(ipt+1).set(yval[AWAY], yerr[AWAY]);
305        }
306      }
307
308    }
309
310    /// @}
311
312  private:
313
314    bool isHI;
315    static const int PT_BINS = 4;
316    static const int EVENT_TYPES = 3;
317    static const int NEAR = 0;
318    static const int AWAY = 1;
319
320    /// @name Histograms
321    /// @{
322    Histo1DPtr _histYield[EVENT_TYPES][PT_BINS];
323    Histo1DPtr _histYieldNoBkg[EVENT_TYPES][PT_BINS];
324    CounterPtr _counterTrigger[EVENT_TYPES];
325    Estimate1DPtr _histIAA[6];
326    /// @}
327
328  };
329
330  RIVET_DECLARE_PLUGIN(ALICE_2012_I930312);
331
332}