96bec6bf1cdbc1cc28c87320a14994a335acc35c
[u/mrichter/AliRoot.git] / PWGCF / EBYE / MeanPtFluctuations / AliAnalysisTaskPtFluc.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TH2F.h"
5 #include "TList.h"
6 #include "TRandom3.h"
7 #include "iostream"
8
9 #include "AliAnalysisTaskSE.h"
10 #include "AliAnalysisManager.h"
11
12 #include "AliESDEvent.h"
13 #include "AliESDInputHandler.h"
14 #include "AliESDtrackCuts.h"
15
16 #include "AliMCEvent.h"
17 #include "AliStack.h"
18 #include "AliMCEventHandler.h"
19
20 #include "AliAnalysisTaskPtFluc.h"
21
22
23 using namespace std;
24
25 // Analysis of Pt Fluctuations (pp)
26 // Author: Stefan Heckel
27 // Version of pp task:   11.4, 11.06.2012
28
29
30 ClassImp(AliAnalysisTaskPtFluc)
31
32 //________________________________________________________________________
33 AliAnalysisTaskPtFluc::AliAnalysisTaskPtFluc(const char *name)
34   :AliAnalysisTaskSE(name),
35   fESD(0),
36   fMCev(0),
37   fRandom3(0),
38   fOutputList(0),
39   fPtSpec(0),
40   fPtSpec2(0),
41   fMult(0),
42   fEta(0),
43   fEtaPhiPlus(0),
44   fEtaPhiMinus(0),
45   fVtxZ(0),
46   fVtxZCut(0),
47   fVtxZCont(0),
48   fVtxZCutDiff(0),
49   fVtxZPileup(0),
50   fVtxZTrackCuts(0),
51   fVtxZDiff1(0),
52   fVtxZDiff2(0),
53   fVtxZDiff3(0),
54   fVtxZDiff1b(0),
55   fVtxZDiff2b(0),
56   fVtxZDiff3b(0),
57   fEventMeanPt(0),
58   fEventMeanPtSq(0),
59   fEventMeanPtMult(0),
60   fMultEventMeanPt(0),
61   fMultEventMeanPtSq(0),
62   fTwoPartCorrEv(0),
63   fTwoPartCorrEvSq(0),
64   fTwoPartCorrEvSample(0),
65   fTwoPartCorrEvSampleSq(0),
66   fESDTrackCuts(0),
67   fMaxVertexZ(0),
68   fMaxVertexZDiff1(0),
69   fNContributors(0),
70   fMC(0),
71   fMCType(0)
72 {
73   DefineOutput(1, TList::Class());
74 }
75
76 //________________________________________________________________________
77 AliAnalysisTaskPtFluc::~AliAnalysisTaskPtFluc()
78 {
79   if (fRandom3) delete fRandom3; fRandom3 = 0;
80
81   if (fOutputList) {
82     fOutputList->Clear();
83     delete fOutputList;
84   }
85   fOutputList = 0;
86
87   if (fESDTrackCuts) delete fESDTrackCuts; fESDTrackCuts = 0;
88
89 }
90
91 //________________________________________________________________________
92 void AliAnalysisTaskPtFluc::UserCreateOutputObjects()
93 {
94   // Create histograms
95   // Called once
96
97   OpenFile(1, "RECREATE");
98   fOutputList = new TList();
99   fOutputList->SetOwner();
100
101
102   fPtSpec = new TH1F("fPtSpec","Pt spectrum",100,0,2.5);
103   fPtSpec2 = new TH1F("fPtSpec2","Pt spectrum 2 - MC ESD",100,0,2.5);
104   fMult = new TH1F("fMult","Multiplicity distribution",120,-0.5,119.5);
105
106   fEta = new TH1F("fEta","Eta distribution",80,-2,2);
107   fEtaPhiPlus = new TH1F("fEtaPhiPlus","Phi distribution for positive eta",62,0,6.2);
108   fEtaPhiMinus = new TH1F("fEtaPhiMinus","Phi distribution for negative eta",62,0,6.2);
109
110   fVtxZ = new TH1F("fVtxZ","Vertex Z distribution before cuts",100,-20,20);
111   fVtxZCut = new TH1F("fVtxZCut","Vertex Z distribution after vtxZ cut",110,-11,11);
112   fVtxZCont = new TH1F("fVtxZCont","Vertex Z distribution after nCont cut",110,-11,11);
113   fVtxZCutDiff = new TH1F("fVtxZCutDiff","Vertex Z distribution after cut on vtx Z Diff",110,-11,11);
114   fVtxZPileup = new TH1F("fVtxZPileup","Vertex Z distribution after pileup rejection",110,-11,11);
115   fVtxZTrackCuts = new TH1F("fVtxZTrackCuts","Vertex Z distribution after track cuts",110,-11,11);
116
117   fVtxZDiff1 = new TH1F("fVtxZDiff1","Difference 1 between vertex Z distributions",100,-5,5);
118   fVtxZDiff2 = new TH1F("fVtxZDiff2","Difference 2 between vertex Z distributions",100,-5,5);
119   fVtxZDiff3 = new TH1F("fVtxZDiff3","Difference 3 between vertex Z distributions",100,-5,5);
120   fVtxZDiff1b = new TH1F("fVtxZDiff1b","Difference 1 between vertex Z distributions after all cuts",100,-5,5);
121   fVtxZDiff2b = new TH1F("fVtxZDiff2b","Difference 2 between vertex Z distributions after all cuts",100,-5,5);
122   fVtxZDiff3b = new TH1F("fVtxZDiff3b","Difference 3 between vertex Z distributions after all cuts",100,-5,5);
123
124
125   fEventMeanPt = new TH1F("fEventMeanPt","Mean-Pt distribution",50,0,2.5);
126   fEventMeanPtSq = new TH1F("fEventMeanPtSq","Mean-Pt squared distribution",100,0,5);
127   fEventMeanPtMult = new TH2F("fEventMeanPtMult","Mean-Pt for single events vs. multiplicity",20,0,100,200,0.,2.);
128
129   fMultEventMeanPt = new TH1F("fMultEventMeanPt","Mean-Pt event by event vs. Multiplicity",100,-0.5,99.5);
130   fMultEventMeanPtSq = new TH1F("fMultEventMeanPtSq","Mean-Pt event by event squared vs. Multiplicity",100,-0.5,99.5);
131
132
133   fTwoPartCorrEv = new TH1F("fTwoPartCorrEv","Two-particle correlator vs. Multiplicity",100,-0.5,99.5);
134   fTwoPartCorrEvSq = new TH1F("fTwoPartCorrEvSq","Two-particle correlator squared vs. Multiplicity",100,-0.5,99.5);
135   fTwoPartCorrEvSample = new TH1F("fTwoPartCorrEvSample","Two-particle correlator vs. Multiplicity",100,-0.5,99.5);
136   fTwoPartCorrEvSampleSq = new TH1F("fTwoPartCorrEvSampleSq","Two-particle correlator squared vs. Multiplicity",100,-0.5,99.5);
137
138
139   // Add histograms to the output list
140   fOutputList->Add(fPtSpec);
141   fOutputList->Add(fPtSpec2);
142   fOutputList->Add(fMult);
143   fOutputList->Add(fEta);
144   fOutputList->Add(fEtaPhiPlus);
145   fOutputList->Add(fEtaPhiMinus);
146   fOutputList->Add(fVtxZ);
147   fOutputList->Add(fVtxZCut);
148   fOutputList->Add(fVtxZCont);
149   fOutputList->Add(fVtxZCutDiff);
150   fOutputList->Add(fVtxZPileup);
151   fOutputList->Add(fVtxZTrackCuts);
152   fOutputList->Add(fVtxZDiff1);
153   fOutputList->Add(fVtxZDiff2);
154   fOutputList->Add(fVtxZDiff3);
155   fOutputList->Add(fVtxZDiff1b);
156   fOutputList->Add(fVtxZDiff2b);
157   fOutputList->Add(fVtxZDiff3b);
158   fOutputList->Add(fEventMeanPt);
159   fOutputList->Add(fEventMeanPtSq);
160   fOutputList->Add(fEventMeanPtMult);
161   fOutputList->Add(fMultEventMeanPt);
162   fOutputList->Add(fMultEventMeanPtSq);
163   fOutputList->Add(fTwoPartCorrEv);
164   fOutputList->Add(fTwoPartCorrEvSq);
165   fOutputList->Add(fTwoPartCorrEvSample);
166   fOutputList->Add(fTwoPartCorrEvSampleSq);
167
168   // Post output data (if histograms are not used later, PostData is at least called here)
169   PostData(1, fOutputList);
170
171 }
172
173 //________________________________________________________________________
174 void AliAnalysisTaskPtFluc::UserExec(Option_t *)
175 {
176   // Main loop
177   // Called for each event
178
179
180   // --- ESD event handler ---
181   if (!fMC || fMCType < 3) {
182
183 //     Printf("ESD handler");
184
185     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
186
187     if (!esdH) {
188       Printf("ERROR: Could not get ESDInputHandler");
189     }
190     else {
191       fESD = esdH->GetEvent();
192     }
193
194     if (!fESD) {
195       Printf("ERROR: fESD not available");
196       return;
197     }
198
199     if(!fESDTrackCuts) Printf("ERROR: No esd track cut");
200
201   } // --- End ESD event handler ---
202
203
204   // --- MC event handler ---
205
206   AliStack *stack = NULL;
207
208   if (fMC) {
209
210 //     Printf("MC handler");
211
212     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
213
214     if (!mcH) {
215       Printf("ERROR: Could not get MCInputHandler");
216       return;
217     }
218     else {
219      fMCev = mcH->MCEvent();
220     }
221
222     if (!fMCev) {
223       Printf("ERROR: fMCev not available");
224       return;
225     }
226
227    stack = fMCev->Stack();
228
229   }
230
231   // --- End MC event handler ---
232
233
234   // --- Vertex cuts ---
235
236   Double_t vtxZGlobal=0., vtxZSPD=0., vtxZTPC=0.;
237   Double_t vtxNContGlobal=0.;
238 //   Double_t  vtxNContSPD=0., vtxNContTPC=0.;
239   Double_t vtxZ=0., vtxNCont=0.;
240   Double_t vtxZdiff1=0., vtxZdiff2=0., vtxZdiff3=0.;
241
242  if (!fMC || fMCType < 3) {
243
244   // Global vertex
245   const AliESDVertex* vtxESD = fESD->GetPrimaryVertexTracks();
246   vtxZGlobal = vtxESD->GetZv();
247   vtxNContGlobal = vtxESD->GetNContributors();
248
249   // SPD vertex
250   const AliESDVertex* vtxESDSPD = fESD->GetPrimaryVertexSPD();
251   vtxZSPD = vtxESDSPD->GetZv();
252 //   vtxNContSPD = vtxESDSPD->GetNContributors();
253
254   // TPC vertex
255   const AliESDVertex* vtxESDTPC = fESD->GetPrimaryVertexTPC();
256   vtxZTPC = vtxESDTPC->GetZv();
257 //   vtxNContTPC = vtxESDTPC->GetNContributors();
258
259   vtxZ = vtxZGlobal;
260   vtxNCont = vtxNContGlobal;
261
262   fVtxZ->Fill(vtxZ); // VtxZ before cuts
263
264   // Differences between vertices
265   vtxZdiff1 = vtxZTPC - vtxZGlobal;
266   vtxZdiff2 = vtxZTPC - vtxZSPD;
267   vtxZdiff3 = vtxZGlobal - vtxZSPD;
268   fVtxZDiff1->Fill(vtxZdiff1);
269   fVtxZDiff2->Fill(vtxZdiff2);
270   fVtxZDiff3->Fill(vtxZdiff3);
271
272   // Event cut on the z-position of the vertex
273   if(vtxZ > fMaxVertexZ || vtxZ < (-1.*fMaxVertexZ)) {
274 //     Printf("VertexZ out of range, Zv = %f",vtxZ);
275     return;
276   }
277
278   fVtxZCut->Fill(vtxZ); // VtxZ after cut on vtxZ
279
280   // Event cut on the number of contributors
281   if(vtxNCont < fNContributors) {
282 //     Printf("Vertex has no contributors");
283     return;
284   }
285
286   fVtxZCont->Fill(vtxZ); // VtxZ after cut on nContributors
287
288   // Event cut on the difference of the z-positions of the vertices (TPC - global)
289   if (fMaxVertexZDiff1 > 0.) {
290     if(vtxZdiff1 > fMaxVertexZDiff1 || vtxZdiff1 < (-1.*fMaxVertexZDiff1)) {
291 //       Printf("VertexZ Diff (TPC - global) out of range, vtxZdiff1 = %f",vtxZ);
292       return;
293     }
294   }
295
296   fVtxZCutDiff->Fill(vtxZ); // VtxZ after cut on vtxZDiff
297
298   } // --- End vertex cuts ---
299
300
301   // --- Pileup rejection ---
302   if (!fMC || fMCType < 3) {
303
304     if (fESD->IsPileupFromSPD(3,0.8,3.,2.,5.)) {
305 //       Printf("IsPileupFromSPD");
306       return;
307     }
308
309     fVtxZPileup->Fill(vtxZ); // VtxZ after pileup rejection
310
311   } // --- End pileup rejection ---
312
313
314   Int_t nrTracks = 0;           // count for number of tracks which pass the track cuts
315 //   Int_t nrTracks2 = 0;               // count for number of tracks which pass the track cuts - for MC ESD vs. MC truth
316   Int_t nrESDTracks = 0;        // number of tracks in the ESD file
317   Int_t nrMCTracks = 0;         // number of tracks in the MC file
318   Double_t sumPt = 0.;          // sum of track Pts
319   Double_t twoPartCorrEv = 0.;  // Two-particle correlator of one event for multiplicity bin analysis
320   Double_t twoPartCorrEvSample = 0.; // Two-part. corr. of one event for whole sample analysis
321
322   Double_t tracks[120] = {0.};  // array of track Pts, needed for the two-particle correlator
323
324   Double_t trackPt=0., trackEta=0., trackPhi=0.;
325   Double_t trackPt2=0.; // for MC ESD
326   Double_t eventMeanPt=0., eventMeanPtSq=0., evMptMult=0.;
327   Double_t twoPartCorrPair=0., twoPartCorrEvSq=0.;
328   Double_t twoPartCorrPairSample=0., twoPartCorrEvSampleSq=0.;
329 //   Double_t nrPairs=0.;
330
331   Double_t *nbins = 0x0; // Mean pT values for multiplicity bin analysis
332   Double_t evMptSample=0.;  // Mean pT value for whole sample analysis
333
334   Double_t random=0., trackPt40=0., ptRatio=0.;
335   Int_t trackPtBin=0;
336
337   fRandom3 = new TRandom3();
338
339   if (fMCType == 2 || fMCType == 4) { // MC, Type = mod. MC truth
340    fRandom3->SetSeed(0);
341   }
342
343
344   Double_t energy = 0.; // beam energy ---> has to be set by the user (in GeV) for MCType >= 3 (w/o ESD) <---
345
346   if (!fMC || fMCType < 3) {
347     energy = fESD->GetBeamEnergy();
348   }
349
350
351 // --- Mean pT values ---
352
353 // !!!!! Have to be calculated in a first run - for each sample, which should be analysed, separately! !!!!!
354
355
356 // -- Mean pT values for the multiplicity bins (starting with N = 0, therefore the first value has to be 0.000) --
357
358 // MC 900 GeV 10e13 (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
359 Double_t nbinsMC900[100] = {0.000, 0.452, 0.468, 0.482, 0.497, 0.511, 0.524, 0.535, 0.544, 0.550, 0.557, 0.563, 0.567, 0.571, 0.576, 0.580, 0.583, 0.586, 0.590, 0.592, 0.595, 0.599, 0.600, 0.602, 0.604, 0.611, 0.604, 0.618, 0.614, 0.612, 0.632, 0.623, 0.617, 0.612, 0.640, 0.632, 0.649, 0.609, 0.658, 0.615, 0.677, 0.614, 0.646, 0.000, 0.000, 0.000, 0.000, 0.000, 0.672, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
360
361 // MC 2.76 TeV 11b10a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
362 Double_t nbinsMC276[100] = {0.000, 0.451, 0.467, 0.485, 0.505, 0.523, 0.540, 0.551, 0.563, 0.569, 0.577, 0.583, 0.588, 0.594, 0.598, 0.603, 0.608, 0.611, 0.615, 0.617, 0.622, 0.622, 0.626, 0.630, 0.633, 0.632, 0.637, 0.633, 0.641, 0.643, 0.650, 0.649, 0.655, 0.650, 0.669, 0.664, 0.653, 0.671, 0.668, 0.683, 0.693, 0.662, 0.765, 0.643, 0.677, 0.660, 0.000, 0.689, 0.708, 0.653, 0.635, 0.000, 0.719, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
363
364 // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
365 Double_t nbinsMC7[100] = {0.000, 0.447, 0.465, 0.485, 0.508, 0.530, 0.548, 0.562, 0.573, 0.581, 0.589, 0.595, 0.601, 0.607, 0.612, 0.617, 0.621, 0.625, 0.629, 0.633, 0.637, 0.640, 0.643, 0.646, 0.649, 0.652, 0.654, 0.657, 0.659, 0.662, 0.664, 0.666, 0.668, 0.670, 0.672, 0.674, 0.677, 0.677, 0.680, 0.682, 0.684, 0.685, 0.687, 0.688, 0.693, 0.692, 0.692, 0.693, 0.696, 0.698, 0.696, 0.700, 0.702, 0.697, 0.701, 0.708, 0.712, 0.714, 0.705, 0.717, 0.702, 0.701, 0.716, 0.716, 0.761, 0.706, 0.750, 0.715, 0.731, 0.728, 0.711, 0.733, 0.000, 0.714, 0.704, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.652, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
366 // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> values for MC truth <---
367 // Double_t nbinsMC7[100] = {0.000, 0.441, 0.455, 0.473, 0.493, 0.515, 0.534, 0.550, 0.561, 0.570, 0.578, 0.585, 0.591, 0.597, 0.602, 0.607, 0.612, 0.616, 0.621, 0.624, 0.628, 0.631, 0.635, 0.638, 0.641, 0.644, 0.646, 0.649, 0.652, 0.654, 0.656, 0.659, 0.661, 0.663, 0.665, 0.666, 0.669, 0.670, 0.672, 0.674, 0.676, 0.678, 0.679, 0.683, 0.683, 0.685, 0.686, 0.686, 0.692, 0.693, 0.690, 0.691, 0.690, 0.692, 0.701, 0.701, 0.700, 0.704, 0.701, 0.704, 0.709, 0.698, 0.707, 0.711, 0.692, 0.739, 0.689, 0.715, 0.691, 0.720, 0.790, 0.670, 0.000, 0.717, 0.833, 0.716, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.619, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
368
369
370 // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
371 // Double_t nbinsMC7[100] = {0.000, 0.483, 0.490, 0.502, 0.515, 0.527, 0.535, 0.541, 0.545, 0.548, 0.552, 0.555, 0.558, 0.561, 0.563, 0.566, 0.568, 0.571, 0.573, 0.574, 0.576, 0.578, 0.579, 0.581, 0.582, 0.583, 0.584, 0.585, 0.585, 0.586, 0.586, 0.587, 0.588, 0.589, 0.588, 0.588, 0.588, 0.588, 0.587, 0.589, 0.588, 0.588, 0.587, 0.587, 0.586, 0.583, 0.588, 0.587, 0.589, 0.583, 0.586, 0.582, 0.584, 0.583, 0.585, 0.582, 0.587, 0.576, 0.580, 0.588, 0.566, 0.581, 0.594, 0.585, 0.531, 0.000, 0.580, 0.499, 0.534, 0.534, 0.580, 0.432, 0.674, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.655, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
372 // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> values for MC truth <---
373 // Double_t nbinsMC7[100] = {0.000, 0.479, 0.483, 0.493, 0.506, 0.518, 0.527, 0.533, 0.537, 0.540, 0.543, 0.546, 0.549, 0.551, 0.554, 0.556, 0.558, 0.560, 0.562, 0.564, 0.566, 0.568, 0.569, 0.570, 0.571, 0.572, 0.573, 0.574, 0.575, 0.576, 0.576, 0.576, 0.577, 0.577, 0.577, 0.579, 0.577, 0.577, 0.577, 0.576, 0.577, 0.577, 0.577, 0.575, 0.575, 0.575, 0.573, 0.575, 0.571, 0.575, 0.568, 0.569, 0.570, 0.576, 0.566, 0.570, 0.568, 0.575, 0.565, 0.562, 0.564, 0.545, 0.563, 0.527, 0.551, 0.540, 0.571, 0.561, 0.562, 0.553, 0.000, 0.586, 0.000, 0.603, 0.480, 0.489, 0.566, 0.513, 0.715, 0.000, 0.000, 0.000, 0.441, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
374
375
376 // Data 900 GeV 10c.pass3 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
377 Double_t nbinsData900[100] = {0.000, 0.488, 0.486, 0.489, 0.495, 0.504, 0.513, 0.521, 0.528, 0.535, 0.541, 0.545, 0.549, 0.554, 0.558, 0.561, 0.565, 0.568, 0.570, 0.573, 0.576, 0.577, 0.580, 0.583, 0.582, 0.587, 0.587, 0.589, 0.591, 0.594, 0.601, 0.599, 0.598, 0.601, 0.608, 0.609, 0.600, 0.606, 0.611, 0.605, 0.612, 0.603, 0.576, 0.606, 0.617, 0.609, 0.572, 0.607, 0.546, 0.000, 0.604, 0.000, 0.689, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
378
379 // Data 2.76 TeV 11a.pass2 w/oSDD (Pt-Range: 0.15 - 2) ---> TPC standalone <---
380 Double_t nbinsData276[100] = {0.000, 0.492, 0.491, 0.496, 0.503, 0.513, 0.523, 0.533, 0.541, 0.549, 0.555, 0.560, 0.565, 0.570, 0.574, 0.578, 0.582, 0.585, 0.589, 0.591, 0.594, 0.597, 0.599, 0.601, 0.604, 0.606, 0.608, 0.610, 0.612, 0.613, 0.615, 0.618, 0.619, 0.621, 0.622, 0.624, 0.625, 0.626, 0.626, 0.628, 0.632, 0.633, 0.633, 0.635, 0.635, 0.638, 0.634, 0.641, 0.642, 0.644, 0.638, 0.646, 0.644, 0.647, 0.644, 0.646, 0.640, 0.650, 0.650, 0.655, 0.666, 0.650, 0.659, 0.644, 0.646, 0.657, 0.629, 0.626, 0.602, 0.652, 0.629, 0.000, 0.659, 0.580, 0.000, 0.000, 0.000, 0.000, 0.000, 0.706, 0.000, 0.596, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
381
382 // Data 7 TeV 10d.pass2 + 10e.pass2 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
383 Double_t nbinsData7[100] = {0.000, 0.492, 0.492, 0.497, 0.505, 0.516, 0.527, 0.538, 0.547, 0.555, 0.562, 0.568, 0.573, 0.578, 0.583, 0.587, 0.591, 0.595, 0.598, 0.601, 0.604, 0.607, 0.610, 0.612, 0.615, 0.617, 0.619, 0.621, 0.623, 0.625, 0.627, 0.629, 0.631, 0.632, 0.634, 0.636, 0.637, 0.639, 0.640, 0.642, 0.643, 0.645, 0.646, 0.647, 0.648, 0.650, 0.650, 0.652, 0.653, 0.654, 0.655, 0.657, 0.657, 0.659, 0.659, 0.659, 0.660, 0.663, 0.663, 0.662, 0.663, 0.665, 0.668, 0.665, 0.667, 0.668, 0.672, 0.671, 0.669, 0.672, 0.672, 0.671, 0.675, 0.666, 0.683, 0.670, 0.673, 0.676, 0.684, 0.672, 0.673, 0.675, 0.689, 0.678, 0.678, 0.720, 0.675, 0.690, 0.664, 0.622, 0.000, 0.698, 0.678, 0.686, 0.751, 0.657, 0.000, 0.663, 0.000, 0.677};
384
385 // -- End mean pT values for the multiplicity bins --
386
387
388 // -- Selection of MC/Data and energy; whole sample values --
389
390 if (fMC) { // - MC -
391
392   if (energy < 1000.) {
393
394 //     Printf(" -- MC, 900 GeV -- ");
395
396     evMptSample = 0.5281; // MC 900 GeV 10e13 (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
397
398     nbins = nbinsMC900;
399   }
400   else if (energy > 1000. && energy < 3000.) {
401
402 //     Printf(" -- MC, 2.76 TeV -- ");
403
404     evMptSample = 0.5574; // MC 2.76 TeV 11b10a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
405
406     nbins = nbinsMC276;
407   }
408   else {
409
410 //     Printf(" -- MC, 7 TeV -- ");
411
412     evMptSample = 0.5809; // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
413 //     evMptSample = 0.5732; // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> value for MC truth <---
414
415 //     evMptSample = 0.5481; // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
416 //     evMptSample = 0.5409; // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> value for MC truth <---
417
418     nbins = nbinsMC7;
419   }
420
421 } // - End MC -
422 else { // - Data -
423
424   if (energy < 1000.) {
425
426 //     Printf(" -- Data, 900 GeV -- ");
427
428     evMptSample = 0.5263; // Data 900 GeV 10c.pass3 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
429
430     nbins = nbinsData900;
431   }
432   else if (energy > 1000. && energy < 3000.) {
433
434 //     Printf(" -- Data, 2.76 TeV -- ");
435
436     evMptSample = 0.5548; // Data 2.76 TeV 11a.pass2 w/oSDD (Pt-Range: 0.15 - 2) ---> TPC standalone <---
437
438     nbins = nbinsData276;
439   }
440   else {
441
442 //     Printf(" -- Data, 7 TeV -- ");
443
444     evMptSample = 0.5738; // Data 7 TeV 10d.pass2 + 10e.pass2 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
445
446     nbins = nbinsData7;
447   }
448
449 } // - End data -
450
451 // -- End selection of MC/Data and energy; whole sample values --
452
453 // --- End mean pT values ---
454
455
456 // --- Ratio of pT distributions MC ESD / MC truth ---
457
458 // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2, array-range: 0 - 2.1 GeV/c)
459 Double_t nbinsPtRatio[84] = {0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.823, 0.908, 0.908, 0.901, 0.903, 0.910, 0.918, 0.931, 0.941, 0.947, 0.951, 0.953, 0.956, 0.959, 0.959, 0.960, 0.961, 0.962, 0.963, 0.963, 0.963, 0.964, 0.964, 0.964, 0.963, 0.963, 0.963, 0.964, 0.964, 0.962, 0.963, 0.964, 0.963, 0.963, 0.962, 0.962, 0.961, 0.963, 0.964, 0.964, 0.965, 0.965, 0.965, 0.964, 0.964, 0.964, 0.964, 0.964, 0.963, 0.964, 0.961, 0.957, 0.958, 0.959, 0.959, 0.958, 0.958, 0.957, 0.959, 0.956, 0.956, 0.954, 0.954, 0.956, 0.955, 0.954, 0.954, 0.954, 0.955, 0.954, 0.951, 0.952, 0.951, 0.949, 0.000, 0.000, 0.000, 0.000};
460
461 // --- End ratio of pT distributions MC ESD / MC truth ---
462
463
464   if (!fMC || fMCType < 3) {
465     nrESDTracks = fESD->GetNumberOfTracks();
466   }
467
468   if (fMC) {
469     nrMCTracks = stack->GetNtrack();
470   }
471
472 //   Printf("\n\n nrESDTracks: %i   nrMCTracks : %i \n",nrESDTracks,nrMCTracks);
473
474
475 if (fMC) { // - MC -
476
477  if (fMCType == 0) { // MC, Type = ESD
478
479
480   // --- Loop over all tracks of one ESD event ---
481   for (Int_t iTracks = 0; iTracks < nrESDTracks; iTracks++) {
482
483     AliESDtrack* track = fESD->GetTrack(iTracks);
484
485     if (!track) {
486       Printf("ERROR: Could not receive track %d\n", iTracks);
487       continue;
488     }
489
490       if(!fESDTrackCuts->AcceptTrack(track))continue;
491
492         trackPt = track->Pt();
493         fPtSpec->Fill(trackPt);
494         tracks[nrTracks] = trackPt;
495
496         trackEta = track->Eta();
497         fEta->Fill(trackEta);
498
499         trackPhi = track->Phi();
500
501         if (trackEta > 0.) {
502           fEtaPhiPlus->Fill(trackPhi);
503         }
504         else if (trackEta < 0.) {
505           fEtaPhiMinus->Fill(trackPhi);
506         }
507
508         sumPt = sumPt + trackPt;
509         nrTracks++;
510
511 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);
512
513   } // --- End track loop ---
514
515
516  } // End MC, Type = ESD
517  else if (fMCType > 0) { // MC, Type = MC truth
518
519
520   // --- Loop over all tracks of one MC event - MC truth / stack ---
521   for (Int_t iTracks = 0; iTracks < nrMCTracks; iTracks++) {
522
523     TParticle *track = stack->Particle(iTracks);
524
525     if (!track) {
526       Printf("ERROR: Could not receive track %d\n", iTracks);
527       continue;
528     }
529
530         // only charged particles
531         TParticlePDG* pdg = track->GetPDG();
532         if(!pdg) continue;
533
534         Double_t charge = pdg->Charge()/3.;
535         if ( TMath::Abs(charge) < 0.001 ) continue;
536
537         // only physical primary
538         Bool_t prim = stack->IsPhysicalPrimary(iTracks);
539         if(!prim) continue;
540
541         // get pt, eta and phi and cut in eta and pt
542         trackPt = track->Pt();
543         trackEta = track->Eta();
544         trackPhi = track->Phi();
545
546 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f  ",trackPt,trackEta,trackPhi);
547
548         if (trackEta > 0.8) continue;
549         if (trackEta < -0.8) continue;
550
551         if (trackPt < 0.15) continue;
552         if (trackPt > 2.) continue;
553
554 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f  ",trackPt,trackEta,trackPhi);
555 //      Printf("-> Track accepted!");
556
557
558         // -- Reject tracks according to ratio of pT distribution MC ESD / MC truth --
559         if (fMCType == 2 || fMCType == 4) { // MC, Type = mod. MC truth
560
561           trackPt40 = trackPt*40.;
562           trackPtBin = floor(trackPt40);
563
564           ptRatio = nbinsPtRatio[trackPtBin];
565
566 //        Printf("trackPt: %f   trackPt40: %f   trackPtBin: %i   ptRatio: %.3f   ",trackPt,trackPt40,trackPtBin,ptRatio);
567
568           random = fRandom3->Rndm();
569 //        Printf("Random: %f  \n",random);
570
571           if (random > ptRatio) continue;
572 //        Printf("accepted random: %f \n",random);
573
574         } // -- End reject tracks according to ratio of pT distribution MC ESD / MC truth --
575
576
577         fPtSpec->Fill(trackPt);
578         tracks[nrTracks] = trackPt;
579
580         fEta->Fill(trackEta);
581
582         if (trackEta > 0.) {
583           fEtaPhiPlus->Fill(trackPhi);
584         }
585         else if (trackEta < 0.) {
586           fEtaPhiMinus->Fill(trackPhi);
587         }
588
589         sumPt = sumPt + trackPt;
590         nrTracks++;
591
592 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);
593
594   } // --- End stack track loop ---
595
596 //   Printf("nrTracks MC truth: %i \n",nrTracks);
597
598
599   // MC, Type = MC truth with corresponding ESD
600   if (fMCType < 3) {
601
602     // --- Loop over all tracks of one MC event - MC ESD ---
603     for (Int_t jTracks = 0; jTracks < nrESDTracks; jTracks++) {
604       AliESDtrack* track = fESD->GetTrack(jTracks);
605       if (!track) {
606         Printf("ERROR: Could not receive track %d\n", jTracks);
607         continue;
608       }
609
610         if(!fESDTrackCuts->AcceptTrack(track))continue;
611
612           trackPt2 = track->Pt();
613           fPtSpec2->Fill(trackPt2);
614
615 //        nrTracks2++;
616
617 //        Printf("Track Pt ESD = %.3f\n\n",trackPt2);
618
619     } // --- End ESD track loop ---
620
621 //     Printf("nrTracks ESD: %i \n\n",nrTracks2);
622
623   } // End MC, Tpye = MC truth with corresponding ESD
624
625  } // End MC, Type = MC truth
626
627
628 } // - End MC -
629 else { // - Data -
630
631
632   // --- Loop over all tracks of one ESD event ---
633   for (Int_t iTracks = 0; iTracks < nrESDTracks; iTracks++) {
634
635     AliESDtrack* track = fESD->GetTrack(iTracks);
636
637     if (!track) {
638       Printf("ERROR: Could not receive track %d\n", iTracks);
639       continue;
640     }
641
642       if(!fESDTrackCuts->AcceptTrack(track))continue;
643
644         trackPt = track->Pt();
645         fPtSpec->Fill(trackPt);
646         tracks[nrTracks] = trackPt;
647
648         trackEta = track->Eta();
649         fEta->Fill(trackEta);
650
651         trackPhi = track->Phi();
652
653         if (trackEta > 0.) {
654           fEtaPhiPlus->Fill(trackPhi);
655         }
656         else if (trackEta < 0.) {
657           fEtaPhiMinus->Fill(trackPhi);
658         }
659
660         sumPt = sumPt + trackPt;
661         nrTracks++;
662
663 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);
664
665   } // --- End track loop ---
666
667 } // - End data -
668
669
670   // --- Calculation of various values and filling of histograms (for remaining events with N_acc > 0) ---
671
672   if(nrTracks != 0) {
673
674     if (!fMC || fMCType < 3) {
675
676       // VtxZ after all track cuts
677       fVtxZTrackCuts->Fill(vtxZ);
678
679       // Differences between vertices after all cuts
680       fVtxZDiff1b->Fill(vtxZdiff1);
681       fVtxZDiff2b->Fill(vtxZdiff2);
682       fVtxZDiff3b->Fill(vtxZdiff3);
683
684     }
685
686     // Multiplicity distribution
687     fMult->Fill(nrTracks);
688
689 //     // Number of pairs in event
690 //     nrPairs = 0.5 * nrTracks * (nrTracks-1);
691
692     // Calculation of mean Pt and mean Pt Squared
693     eventMeanPt = sumPt / nrTracks;
694     eventMeanPtSq = eventMeanPt * eventMeanPt;
695
696     // Mean-Pt and Mean-Pt squared
697     fEventMeanPt->Fill(eventMeanPt);
698     fEventMeanPtSq->Fill(eventMeanPtSq);
699     fEventMeanPtMult->Fill(nrTracks,eventMeanPt);
700
701     // Mean-Pt and Mean-Pt squared depending on multiplicity
702     fMultEventMeanPt->Fill(nrTracks,eventMeanPt);
703     fMultEventMeanPtSq->Fill(nrTracks,eventMeanPtSq);
704
705 //     Printf("nrTracks: %i   sumPt: %.8f   meanPt: %.8f   meanPtSq: %.8f\n",nrTracks,sumPt,eventMeanPt,eventMeanPtSq);
706
707
708     // --- Two-particle correlator ---
709
710     evMptMult = nbins[nrTracks];
711 //     Printf("nrTracks = %3i   evMptMult = %.3f   evMptSample = %.3f \n\n",nrTracks,evMptMult,evMptSample);
712
713     for (int i=0; i<nrTracks; i++) {
714
715       for (int j=i+1; j<nrTracks; j++) {
716
717         twoPartCorrPair = (tracks[i] - evMptMult) * (tracks[j] - evMptMult); // Multiplicity bins
718         twoPartCorrEv = twoPartCorrEv + twoPartCorrPair;
719
720         twoPartCorrPairSample = (tracks[i] - evMptSample) * (tracks[j] - evMptSample); // Whole sample
721         twoPartCorrEvSample = twoPartCorrEvSample + twoPartCorrPairSample;
722
723       }
724     }
725
726     // Multiplicity bins
727     twoPartCorrEvSq = twoPartCorrEv * twoPartCorrEv;
728     fTwoPartCorrEv->Fill(nrTracks,twoPartCorrEv);
729     fTwoPartCorrEvSq->Fill(nrTracks,twoPartCorrEvSq);
730
731     // Whole sample
732     twoPartCorrEvSampleSq = twoPartCorrEvSample * twoPartCorrEvSample;
733     fTwoPartCorrEvSample->Fill(nrTracks,twoPartCorrEvSample);
734     fTwoPartCorrEvSampleSq->Fill(nrTracks,twoPartCorrEvSampleSq);
735
736     // --- End two-particle correlator ---
737
738
739   } // --- End calculation of various values and filling of histograms ---
740
741
742   // Post output data
743   PostData(1, fOutputList);
744
745   // Delet pointer
746   if (fRandom3) delete fRandom3; fRandom3 = 0;
747
748 }
749
750 void AliAnalysisTaskPtFluc::Terminate(Option_t *)
751 {
752   // Called once at the end of the query
753
754   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
755   if (!fOutputList) {
756     Printf("ERROR: fOutputList not available\n");
757     return;
758   }
759
760 }