Rivet analyses referenceALICE_2012_I930312Particle-yield modification in jet-like azimuthal di-hadron correlations in Pb-Pb collisions at $\sqrt{s_\rm{NN}} = 2.76$ TeVExperiment: ALICE (LHC) Inspire ID: 930312 Status: VALIDATED Authors:
Beam energies: (1380.0, 1380.0); (287040.0, 287040.0) GeV
The yield of charged particles associated with high-pT trigger particles ($8 < p_{\perp} < 15$ GeV/c) is measured with the ALICE detector in Pb-Pb collisions at $\sqrt{s_{NN}} = 2.76$ TeV relative to proton-proton collisions at the same energy. The conditional per-trigger yields are extracted from the narrow jet-like correlation peaks in azimuthal di-hadron correlations. In the 5\% most central collisions, we observe that the yield of associated charged particles with transverse momenta $p_\perp > 3$ GeV/c on the away-side drops to about 60\% of that observed in pp collisions, while on the near-side a moderate enhancement of 20-30\% is found. Source code: ALICE_2012_I930312.cc 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}
|