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