AliPhysics  a76316e (a76316e)
AliAnalysisTaskFlowModesMC.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // AliAnalysisTaskFlowModesMC - ALICE Flow framework
17 //
18 // ALICE analysis task for universal study of flow.
19 // Note: So far implemented only for AOD analysis!
20 // Author: Vojtech Pacik
21 // Modified by: Naghmeh Mohammadi to include non-linear terms, Nikhef, 2017
22 // Generic framework by: You Zhou
23 
24 #include <vector>
25 
26 #include <TDatabasePDG.h>
27 #include <TPDGCode.h>
28 
29 #include "TFile.h"
30 #include "TChain.h"
31 #include "TH1D.h"
32 #include "TH2D.h"
33 #include "TH3D.h"
34 #include "TProfile.h"
35 #include "TProfile2D.h"
36 #include "TProfile3D.h"
37 #include "TList.h"
38 #include "TComplex.h"
39 #include "TRandom3.h"
40 
41 #include "AliAnalysisTask.h"
42 #include "AliAnalysisManager.h"
43 #include "AliAODInputHandler.h"
44 #include "AliAnalysisUtils.h"
45 #include "AliMultSelection.h"
46 #include "AliPIDResponse.h"
47 #include "AliPIDCombined.h"
48 #include "AliAnalysisTaskFlowModesMC.h"
49 #include "AliLog.h"
50 #include "AliAODEvent.h"
51 #include "AliESDEvent.h"
52 #include "AliMCEvent.h"
53 #include "AliAODTrack.h"
54 #include "AliAODv0.h"
55 #include "AliVTrack.h"
56 #include "AliESDpid.h"
57 #include "AliFlowBayesianPID.h"
58 #include "AliCollisionGeometry.h"
59 #include "AliGenEventHeader.h"
60 #include "AliGenAmpt.h"
61 
62 
63 class AliAnalysisTaskFlowModesMC;
64 
65 ClassImp(AliAnalysisTaskFlowModesMC); // classimp: necessary for root
66 
67 Int_t AliAnalysisTaskFlowModesMC::fHarmonics[] = {2,3,4,5,6};
68 Int_t AliAnalysisTaskFlowModesMC::fMixedHarmonics[] = {422,633,523,6222};// 422: v4{psi2}, 523: v5{psi2, psi3}, 633: v6{psi3} and 6222: v6{psi2}
69 
70 Double_t AliAnalysisTaskFlowModesMC::fEtaGap[] = {0.};//,0.4,0.8};
71 
72 
73 
74 AliAnalysisTaskFlowModesMC::AliAnalysisTaskFlowModesMC() : AliAnalysisTaskSE(),
75 fEventAOD(0x0),
76 fEventMC(0x0),
77 fPIDResponse(0x0),
78 fPIDCombined(0x0),
79 fFlowNUAWeightsFile(0x0),
80 fFlowNUEWeightsFile(0x0),
81 fInit(kFALSE),
82 fIndexSampling(0),
83 fIndexCentrality(-1),
84 fIndexImpactParameter(-1),
85 fEventCounter(0),
86 fNumEventsAnalyse(50),
87 fRunNumber(-1),
88 fExtraPileUp(kFALSE),
89 fPDGMassPion(TDatabasePDG::Instance()->GetParticle(211)->Mass()),
90 fPDGMassKaon(TDatabasePDG::Instance()->GetParticle(321)->Mass()),
91 fPDGMassProton(TDatabasePDG::Instance()->GetParticle(2212)->Mass()),
92 
93 // FlowPart containers
94 fVectorCharged(0x0),
95 fVectorPion(0x0),
96 fVectorKaon(0x0),
97 fVectorProton(0x0),
98 
99 // analysis selection
100 fRunMode(kFull),
101 fAnalType(kMC),
102 fSampling(kFALSE),
103 fCutImpactParameter(kTRUE),
104 fFillQA(kTRUE),
105 fProcessCharged(kFALSE),
106 fProcessPID(kFALSE),
107 fBayesianResponse(NULL),
108 
109 // flow related
110 fCutFlowRFPsPtMin(0),
111 fCutFlowRFPsPtMax(0),
112 fFlowPOIsPtMin(0),
113 fFlowPOIsPtMax(20.),
114 fFlowCentMin(0),
115 fFlowCentMax(100),
116 fFlowCentNumBins(1000),
117 fImpactParameterMin(0),
118 fImpactParameterMax(50),
119 fCutFlowDoFourCorrelations(kTRUE),
120 fDoOnlyMixedCorrelations(kFALSE),
121 fFlowFillWeights(kFALSE),
122 fFlowUseNUAWeights(kFALSE),
123 fFlowUseNUEWeights(kFALSE),
124 fFlowNUAWeightsPath(),
125 fFlowNUEWeightsPath(),
126 fPositivelyChargedRef(kFALSE),
127 fNegativelyChargedRef(kFALSE),
128 fPositivelyChargedPOI(kFALSE),
129 fNegativelyChargedPOI(kFALSE),
130 
131 
132 
133 // events selection
134 fPVtxCutZ(0.),
135 fColSystem(kPbPb),
136 fMultEstimator(),
137 fTrigger(0),
138 fFullCentralityRange(kTRUE),
139 // charged tracks selection
140 fCutChargedEtaMax(0),
141 fCutChargedPtMax(0),
142 fCutChargedPtMin(0),
143 fCutChargedDCAzMax(0),
144 fCutChargedDCAxyMax(0),
145 fCutChargedTrackFilterBit(0),
146 fCutChargedNumTPCclsMin(0),
147 fMaxChi2perTPCcls(0),
148 
149 // PID tracks selection
150 fESDpid(),
151 fCutPIDUseAntiProtonOnly(kFALSE),
152 fPIDnsigma(kTRUE),
153 fPIDnsigmaCombination(2),
154 fPIDbayesian(kFALSE),
155 fCutPIDnSigmaPionMax(3),
156 fCutPIDnSigmaKaonMax(3),
157 fCutPIDnSigmaProtonMax(3),
158 fCutPIDnSigmaTPCRejectElectron(3),
159 fCutPIDnSigmaCombinedNoTOFrejection(kFALSE),
160 fCurrCentr(0.0),
161 fParticleProbability(0.9),
162 
163 // output lists
164 fQAEvents(0x0),
165 fQACharged(0x0),
166 fQAPID(0x0),
167 fFlowWeights(0x0),
168 fFlowRefs(0x0),
169 fFlowCharged(0x0),
170 fFlowPID(0x0),
171 
172 // flow histograms & profiles
173 fh3BeforeNUAWeightsRefs(0x0),
174 fh3BeforeNUAWeightsCharged(0x0),
175 fh3BeforeNUAWeightsPion(0x0),
176 fh3BeforeNUAWeightsKaon(0x0),
177 fh3BeforeNUAWeightsProton(0x0),
178 
179 fh3AfterNUAWeightsRefs(0x0),
180 fh3AfterNUAWeightsCharged(0x0),
181 fh3AfterNUAWeightsPion(0x0),
182 fh3AfterNUAWeightsKaon(0x0),
183 fh3AfterNUAWeightsProton(0x0),
184 
185 fhBeforeNUEWeightsRefs(0x0),
186 fhBeforeNUEWeightsCharged(0x0),
187 fhBeforeNUEWeightsPion(0x0),
188 fhBeforeNUEWeightsKaon(0x0),
189 fhBeforeNUEWeightsProton(0x0),
190 
191 fhAfterNUEWeightsRefs(0x0),
192 fhAfterNUEWeightsCharged(0x0),
193 fhAfterNUEWeightsPion(0x0),
194 fhAfterNUEWeightsKaon(0x0),
195 fhAfterNUEWeightsProton(0x0),
196 
197 fh3NUAWeightRefsPlus(0x0),
198 fh3NUAWeightChargedPlus(0x0),
199 fh3NUAWeightPionPlus(0x0),
200 fh3NUAWeightKaonPlus(0x0),
201 fh3NUAWeightProtonPlus(0x0),
202 
203 fh3NUAWeightRefsMinus(0x0),
204 fh3NUAWeightChargedMinus(0x0),
205 fh3NUAWeightPionMinus(0x0),
206 fh3NUAWeightKaonMinus(0x0),
207 fh3NUAWeightProtonMinus(0x0),
208 
209 // event histograms
210 fhEventSampling(0x0),
211 fhEventCentrality(0x0),
212 fh2EventCentralityNumSelCharged(0x0),
213 fhEventCounter(0x0),
214 
215 // charged histogram
216 fh2RefsMult(0x0),
217 fh2RefsPt(0x0),
218 fh2RefsEta(0x0),
219 fh2RefsPhi(0x0),
220 fhChargedCounter(0x0),
221 
222 // PID histogram
223 fh2PIDPionMult(0x0),
224 fh2PIDPionPt(0x0),
225 fh2PIDPionPhi(0x0),
226 fh2PIDPionEta(0x0),
227 fhPIDPionCharge(0x0),
228 fh2PIDKaonMult(0x0),
229 fh2PIDKaonPt(0x0),
230 fh2PIDKaonPhi(0x0),
231 fh2PIDKaonEta(0x0),
232 fhPIDKaonCharge(0x0),
233 fh2PIDProtonMult(0x0),
234 fh2PIDProtonPt(0x0),
235 fh2PIDProtonPhi(0x0),
236 fh2PIDProtonEta(0x0),
237 fhPIDProtonCharge(0x0),
238 fh2PIDPionTPCdEdx(0x0),
239 fh2PIDPionTOFbeta(0x0),
240 fh2PIDKaonTPCdEdx(0x0),
241 fh2PIDKaonTOFbeta(0x0),
242 fh2PIDProtonTPCdEdx(0x0),
243 fh2PIDProtonTOFbeta(0x0),
244 fh2PIDPionTPCnSigmaPion(0x0),
245 fh2PIDPionTOFnSigmaPion(0x0),
246 fh2PIDPionTPCnSigmaKaon(0x0),
247 fh2PIDPionTOFnSigmaKaon(0x0),
248 fh2PIDPionTPCnSigmaProton(0x0),
249 fh2PIDPionTOFnSigmaProton(0x0),
250 fh2PIDKaonTPCnSigmaPion(0x0),
251 fh2PIDKaonTOFnSigmaPion(0x0),
252 fh2PIDKaonTPCnSigmaKaon(0x0),
253 fh2PIDKaonTOFnSigmaKaon(0x0),
254 fh2PIDKaonTPCnSigmaProton(0x0),
255 fh2PIDKaonTOFnSigmaProton(0x0),
256 fh2PIDProtonTPCnSigmaPion(0x0),
257 fh2PIDProtonTOFnSigmaPion(0x0),
258 fh2PIDProtonTPCnSigmaKaon(0x0),
259 fh2PIDProtonTOFnSigmaKaon(0x0),
260 fh2PIDProtonTPCnSigmaProton(0x0),
261 fh2PIDProtonTOFnSigmaProton(0x0)
262 
263 {
264  SetPriors(); //init arrays
265  // New PID procedure (Bayesian Combined PID)
266  // allocating here is necessary because we don't
267  // stream this member
268  fBayesianResponse = new AliFlowBayesianPID();
269  fBayesianResponse->SetNewTrackParam();
270  // default constructor, don't allocate memory here!
271  // this is used by root for IO purposes, it needs to remain empty
272 }
273 //_____________________________________________________________________________
274 AliAnalysisTaskFlowModesMC::AliAnalysisTaskFlowModesMC(const char* name) : AliAnalysisTaskSE(name),
275 fEventAOD(0x0),
276 fEventMC(0x0),
277 fPIDResponse(0x0),
278 fPIDCombined(0x0),
279 fFlowNUAWeightsFile(0x0),
280 fFlowNUEWeightsFile(0x0),
281 fInit(kFALSE),
282 fIndexSampling(0),
283 fIndexCentrality(-1),
284 fIndexImpactParameter(-1),
285 fEventCounter(0),
286 fNumEventsAnalyse(50),
287 fRunNumber(-1),
288 fExtraPileUp(kFALSE),
289 fPDGMassPion(TDatabasePDG::Instance()->GetParticle(211)->Mass()),
290 fPDGMassKaon(TDatabasePDG::Instance()->GetParticle(321)->Mass()),
291 fPDGMassProton(TDatabasePDG::Instance()->GetParticle(2212)->Mass()),
292 
293 // FlowPart containers
294 fVectorCharged(0x0),
295 fVectorPion(0x0),
296 fVectorKaon(0x0),
297 fVectorProton(0x0),
298 
299 // analysis selection
300 fRunMode(kFull),
301 fAnalType(kMC),
302 fSampling(kFALSE),
303 fCutImpactParameter(kTRUE),
304 fFillQA(kTRUE),
305 fProcessCharged(kFALSE),
306 fProcessPID(kFALSE),
307 fBayesianResponse(NULL),
308 
309 // flow related
310 fCutFlowRFPsPtMin(0),
311 fCutFlowRFPsPtMax(0),
312 fFlowPOIsPtMin(0),
313 fFlowPOIsPtMax(20.),
314 fFlowCentMin(0),
315 fFlowCentMax(100),
316 fFlowCentNumBins(1000),
317 fImpactParameterMin(0),
318 fImpactParameterMax(50),
319 fCutFlowDoFourCorrelations(kTRUE),
320 fDoOnlyMixedCorrelations(kFALSE),
321 fFlowFillWeights(kFALSE),
322 fFlowUseNUAWeights(kFALSE),
323 fFlowUseNUEWeights(kFALSE),
324 fFlowNUAWeightsPath(),
325 fFlowNUEWeightsPath(),
326 fPositivelyChargedRef(kFALSE),
327 fNegativelyChargedRef(kFALSE),
328 fPositivelyChargedPOI(kFALSE),
329 fNegativelyChargedPOI(kFALSE),
330 
331 // events selection
332 fPVtxCutZ(0.),
333 fColSystem(kPbPb),
334 fMultEstimator(),
335 fTrigger(0),
336 fFullCentralityRange(kTRUE),
337 // charged tracks selection
338 fCutChargedEtaMax(0),
339 fCutChargedPtMax(0),
340 fCutChargedPtMin(0),
341 fCutChargedDCAzMax(0),
342 fCutChargedDCAxyMax(0),
343 fCutChargedTrackFilterBit(0),
344 fCutChargedNumTPCclsMin(0),
345 fMaxChi2perTPCcls(0),
346 // PID tracks selection
347 fESDpid(),
348 fCutPIDUseAntiProtonOnly(kFALSE),
349 fPIDnsigma(kTRUE),
350 fPIDnsigmaCombination(2.),
351 fPIDbayesian(kFALSE),
352 fCutPIDnSigmaPionMax(3),
353 fCutPIDnSigmaKaonMax(3),
354 fCutPIDnSigmaProtonMax(3),
355 fCutPIDnSigmaTPCRejectElectron(3),
356 fCutPIDnSigmaCombinedNoTOFrejection(kFALSE),
357 fCurrCentr(0.0),
358 fParticleProbability(0.9),
359 
360 
361 // output lists
362 fQAEvents(0x0),
363 fQACharged(0x0),
364 fQAPID(0x0),
365 fFlowWeights(0x0),
366 fFlowRefs(0x0),
367 fFlowCharged(0x0),
368 fFlowPID(0x0),
369 
370 // flow histograms & profiles
371 fh3BeforeNUAWeightsRefs(0x0),
372 fh3BeforeNUAWeightsCharged(0x0),
373 fh3BeforeNUAWeightsPion(0x0),
374 fh3BeforeNUAWeightsKaon(0x0),
375 fh3BeforeNUAWeightsProton(0x0),
376 
377 fh3AfterNUAWeightsRefs(0x0),
378 fh3AfterNUAWeightsCharged(0x0),
379 fh3AfterNUAWeightsPion(0x0),
380 fh3AfterNUAWeightsKaon(0x0),
381 fh3AfterNUAWeightsProton(0x0),
382 
383 fhBeforeNUEWeightsRefs(0x0),
384 fhBeforeNUEWeightsCharged(0x0),
385 fhBeforeNUEWeightsPion(0x0),
386 fhBeforeNUEWeightsKaon(0x0),
387 fhBeforeNUEWeightsProton(0x0),
388 
389 fhAfterNUEWeightsRefs(0x0),
390 fhAfterNUEWeightsCharged(0x0),
391 fhAfterNUEWeightsPion(0x0),
392 fhAfterNUEWeightsKaon(0x0),
393 fhAfterNUEWeightsProton(0x0),
394 
395 fh3NUAWeightRefsPlus(0x0),
396 fh3NUAWeightChargedPlus(0x0),
397 fh3NUAWeightPionPlus(0x0),
398 fh3NUAWeightKaonPlus(0x0),
399 fh3NUAWeightProtonPlus(0x0),
400 
401 fh3NUAWeightRefsMinus(0x0),
402 fh3NUAWeightChargedMinus(0x0),
403 fh3NUAWeightPionMinus(0x0),
404 fh3NUAWeightKaonMinus(0x0),
405 fh3NUAWeightProtonMinus(0x0),
406 
407 // event histograms
408 fhEventSampling(0x0),
409 fhEventCentrality(0x0),
410 fh2EventCentralityNumSelCharged(0x0),
411 fhEventCounter(0x0),
412 
413 // charged histogram
414 fh2RefsMult(0x0),
415 fh2RefsPt(0x0),
416 fh2RefsEta(0x0),
417 fh2RefsPhi(0x0),
418 fhChargedCounter(0x0),
419 
420 // PID histogram
421 fh2PIDPionMult(0x0),
422 fh2PIDPionPt(0x0),
423 fh2PIDPionPhi(0x0),
424 fh2PIDPionEta(0x0),
425 fhPIDPionCharge(0x0),
426 fh2PIDKaonMult(0x0),
427 fh2PIDKaonPt(0x0),
428 fh2PIDKaonPhi(0x0),
429 fh2PIDKaonEta(0x0),
430 fhPIDKaonCharge(0x0),
431 fh2PIDProtonMult(0x0),
432 fh2PIDProtonPt(0x0),
433 fh2PIDProtonPhi(0x0),
434 fh2PIDProtonEta(0x0),
435 fhPIDProtonCharge(0x0),
436 fh2PIDPionTPCdEdx(0x0),
437 fh2PIDPionTOFbeta(0x0),
438 fh2PIDKaonTPCdEdx(0x0),
439 fh2PIDKaonTOFbeta(0x0),
440 fh2PIDProtonTPCdEdx(0x0),
441 fh2PIDProtonTOFbeta(0x0),
442 fh2PIDPionTPCnSigmaPion(0x0),
443 fh2PIDPionTOFnSigmaPion(0x0),
444 fh2PIDPionTPCnSigmaKaon(0x0),
445 fh2PIDPionTOFnSigmaKaon(0x0),
446 fh2PIDPionTPCnSigmaProton(0x0),
447 fh2PIDPionTOFnSigmaProton(0x0),
448 fh2PIDKaonTPCnSigmaPion(0x0),
449 fh2PIDKaonTOFnSigmaPion(0x0),
450 fh2PIDKaonTPCnSigmaKaon(0x0),
451 fh2PIDKaonTOFnSigmaKaon(0x0),
452 fh2PIDKaonTPCnSigmaProton(0x0),
453 fh2PIDKaonTOFnSigmaProton(0x0),
454 fh2PIDProtonTPCnSigmaPion(0x0),
455 fh2PIDProtonTOFnSigmaPion(0x0),
456 fh2PIDProtonTPCnSigmaKaon(0x0),
457 fh2PIDProtonTOFnSigmaKaon(0x0),
458 fh2PIDProtonTPCnSigmaProton(0x0),
459 fh2PIDProtonTOFnSigmaProton(0x0)
460 {
461 
462  SetPriors(); //init arrays
463  // New PID procedure (Bayesian Combined PID)
464  fBayesianResponse = new AliFlowBayesianPID();
465  fBayesianResponse->SetNewTrackParam();
466  for(Short_t iCent(0); iCent<fNumCentralities; iCent++){
467  fhNUEWeightRefsPlus[iCent] = 0x0;
468  fhNUEWeightRefsMinus[iCent] = 0x0;
469  fhNUEWeightChargedPlus[iCent] = 0x0;
470  fhNUEWeightChargedMinus[iCent] = 0x0;
471  fhNUEWeightPionPlus[iCent] = 0x0;
472  fhNUEWeightPionMinus[iCent] = 0x0;
473  fhNUEWeightKaonPlus[iCent] = 0x0;
474  fhNUEWeightKaonMinus[iCent] = 0x0;
475  fhNUEWeightProtonPlus[iCent] = 0x0;
476  fhNUEWeightProtonMinus[iCent] = 0x0;
477  }
478 
479  // Flow vectors
480  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
481  {
482  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
483  {
484  fFlowVecQpos[iHarm][iPower] = TComplex(0,0,kFALSE);
485  fFlowVecQneg[iHarm][iPower] = TComplex(0,0,kFALSE);
486 
487  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
488  {
489  fFlowVecPpos[iHarm][iPower][iPt] = TComplex(0,0,kFALSE);
490  fFlowVecPneg[iHarm][iPower][iPt] = TComplex(0,0,kFALSE);
491  fFlowVecS[iHarm][iPower][iPt] = TComplex(0,0,kFALSE);
492  }//endfor(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
493  }//endfor(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
494  }//endfor(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
495 
496  // Flow profiles & histograms
497  if(fDoOnlyMixedCorrelations){
498  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics; iMixedHarm++)
499  {
500  for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
501  {
502  for(Short_t iSample(0); iSample < fNumSamples; iSample++)
503  {
504  fpMixedRefsCor4[iSample][iGap][iMixedHarm] = 0x0;
505  if(iMixedHarm==3)fpMixedRefsCor6[iSample][iGap] = 0x0;
506 
507  fpMixedChargedCor3Pos[iSample][iGap][iMixedHarm] = 0x0;
508  fpMixedChargedCor3Neg[iSample][iGap][iMixedHarm] = 0x0;
509 
510  fpMixedPionCor3Pos[iSample][iGap][iMixedHarm] = 0x0;
511  fpMixedPionCor3Neg[iSample][iGap][iMixedHarm] = 0x0;
512  fpMixedKaonCor3Pos[iSample][iGap][iMixedHarm] = 0x0;
513  fpMixedKaonCor3Neg[iSample][iGap][iMixedHarm] = 0x0;
514  fpMixedProtonCor3Pos[iSample][iGap][iMixedHarm] = 0x0;
515  fpMixedProtonCor3Neg[iSample][iGap][iMixedHarm] = 0x0;
516 
517  if(iMixedHarm==3){
518  fpMixedChargedCor4Pos[iSample][iGap] = 0x0;
519  fpMixedChargedCor4Neg[iSample][iGap] = 0x0;
520 
521  fpMixedPionCor4Pos[iSample][iGap] = 0x0;
522  fpMixedPionCor4Neg[iSample][iGap] = 0x0;
523  fpMixedKaonCor4Pos[iSample][iGap] = 0x0;
524  fpMixedKaonCor4Neg[iSample][iGap] = 0x0;
525  fpMixedProtonCor4Pos[iSample][iGap] = 0x0;
526  fpMixedProtonCor4Neg[iSample][iGap] = 0x0;
527  }
528  }
529 
530  }//endfor(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
531  }//endfor(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics; iMixedHarm++)
532  }//endif(fDoOnlyMixedCorrelations)
533 
534  if(!fDoOnlyMixedCorrelations){
535  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++){
536  //fpRefsCor4[iHarm] = 0x0;
537  //fp2ChargedCor4[iHarm] = 0x0;
538  //fp2PionCor4[iHarm] = 0x0;
539  //fp2KaonCor4[iHarm] = 0x0;
540  //fp2ProtonCor4[iHarm] = 0x0;
541  for(Short_t iGap(0); iGap < fNumEtaGap; iGap++){
542  // mean Qx,Qy
543  fpMeanQxRefsPos[iGap][iHarm] = 0x0;
544  fpMeanQxRefsNeg[iGap][iHarm] = 0x0;
545  fpMeanQyRefsPos[iGap][iHarm] = 0x0;
546  fpMeanQyRefsNeg[iGap][iHarm] = 0x0;
547  for(Short_t iSample(0); iSample < fNumSamples; iSample++)
548  {
549 
550  fpRefsCor2[iSample][iGap][iHarm] = 0x0;
551  fp2ChargedCor2Pos[iSample][iGap][iHarm] = 0x0;
552  fp2ChargedCor2Neg[iSample][iGap][iHarm] = 0x0;
553  fp2PionCor2Pos[iSample][iGap][iHarm] = 0x0;
554  fp2PionCor2Neg[iSample][iGap][iHarm] = 0x0;
555  fp2KaonCor2Pos[iSample][iGap][iHarm] = 0x0;
556  fp2KaonCor2Neg[iSample][iGap][iHarm] = 0x0;
557  fp2ProtonCor2Pos[iSample][iGap][iHarm] = 0x0;
558  fp2ProtonCor2Neg[iSample][iGap][iHarm] = 0x0;
559  }
560  }//endfor(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
561  }//endfor(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
562  }//endif(!fDoOnlyMixedCorrelations)
563 
564  // QA histograms
565  for(Short_t iQA(0); iQA < fiNumIndexQA; iQA++)
566  {
567  // Event histograms
568  fhQAEventsPVz[iQA] = 0x0;
569  fhQAEventsNumContrPV[iQA] = 0x0;
570  fhQAEventsNumSPDContrPV[iQA] = 0x0;
571  fhQAEventsDistPVSPD[iQA] = 0x0;
572  fhQAEventsSPDresol[iQA] = 0x0;
573  fhQAEventsPileUp[iQA] = 0x0;
574  fhQAEventsCentralityOutliers[iQA] = 0x0;
575  fhEventsMultTOFFilterbit32[iQA] = 0x0;
576  // charged
577  fhQAChargedMult[iQA] = 0x0;
578  fhQAChargedPt[iQA] = 0x0;
579  fhQAChargedEta[iQA] = 0x0;
580  fhQAChargedPhi[iQA] = 0x0;
581  fhQAChargedCharge[iQA] = 0x0;
582  fhQAChargedFilterBit[iQA] = 0x0;
583  fhQAChargedNumTPCcls[iQA] = 0x0;
584  fhQAChargedDCAxy[iQA] = 0x0;
585  fhQAChargedDCAz[iQA] = 0x0;
586 
587  // PID
588  fhQAPIDTPCstatus[iQA] = 0x0;
589  fhQAPIDTOFstatus[iQA] = 0x0;
590  fhQAPIDTPCdEdx[iQA] = 0x0;
591  fhQAPIDTOFbeta[iQA] = 0x0;
592 
593  fh3PIDPionTPCTOFnSigmaPion[iQA] = 0x0;
594  fh3PIDPionTPCTOFnSigmaKaon[iQA] = 0x0;
595  fh3PIDPionTPCTOFnSigmaProton[iQA] = 0x0;
596 
597  fh3PIDKaonTPCTOFnSigmaPion[iQA] = 0x0;
598  fh3PIDKaonTPCTOFnSigmaKaon[iQA] = 0x0;
599  fh3PIDKaonTPCTOFnSigmaProton[iQA] = 0x0;
600 
601  fh3PIDProtonTPCTOFnSigmaPion[iQA] = 0x0;
602  fh3PIDProtonTPCTOFnSigmaKaon[iQA] = 0x0;
603  fh3PIDProtonTPCTOFnSigmaProton[iQA] = 0x0;
604 
605 
606  }//endfor(Short_t iQA(0); iQA < fiNumIndexQA; iQA++)
607 
608  // defining input/output
609  DefineInput(0, TChain::Class());
610  DefineOutput(1, TList::Class());
611  DefineOutput(2, TList::Class());
612  DefineOutput(3, TList::Class());
613  DefineOutput(4, TList::Class());
614  DefineOutput(5, TList::Class());
615  DefineOutput(6, TList::Class());
616  DefineOutput(7, TList::Class());
617 }
618 //_____________________________________________________________________________
619 AliAnalysisTaskFlowModesMC::~AliAnalysisTaskFlowModesMC()
620 {
621  // destructor
622  // if(fPIDCombined)
623  // {
624  // delete fPIDCombined;
625  // }
626 
627  // deleting FlowPart vectors (containers)
628  if(fVectorCharged) delete fVectorCharged;
629  if(fVectorPion) delete fVectorPion;
630  if(fVectorKaon) delete fVectorKaon;
631  if(fVectorProton) delete fVectorProton;
632 
633  // deleting output lists
634  if(fFlowWeights) delete fFlowWeights;
635  if(fFlowRefs) delete fFlowRefs;
636  if(fFlowCharged) delete fFlowCharged;
637  if(fFlowPID) delete fFlowPID;
638 
639  if(fQAEvents) delete fQAEvents;
640  if(fQACharged) delete fQACharged;
641  if(fQAPID) delete fQAPID;
642 }
643 //_____________________________________________________________________________
644 void AliAnalysisTaskFlowModesMC::UserCreateOutputObjects()
645 {
646  // create output objects
647  // this function is called ONCE at the start of your analysis (RUNTIME)
648  // *************************************************************
649 
650  // list all parameters used in this analysis
651  ListParameters();
652 
653  // task initialization
654  fInit = InitializeTask();
655  if(!fInit) return;
656 
657  // creating output lists
658  fFlowRefs = new TList();
659  fFlowRefs->SetOwner(kTRUE);
660  fFlowRefs->SetName("fFlowRefs");
661  fFlowCharged = new TList();
662  fFlowCharged->SetOwner(kTRUE);
663  fFlowCharged->SetName("fFlowCharged");
664  fFlowPID = new TList();
665  fFlowPID->SetOwner(kTRUE);
666  fFlowPID->SetName("fFlowPID");
667  fFlowWeights = new TList();
668  fFlowWeights->SetOwner(kTRUE);
669  fFlowWeights->SetName("fFlowWeights");
670 
671  fQAEvents = new TList();
672  fQAEvents->SetOwner(kTRUE);
673  fQACharged = new TList();
674  fQACharged->SetOwner(kTRUE);
675  fQAPID = new TList();
676  fQAPID->SetOwner(kTRUE);
677 
678  // creating particle vectors
679  fVectorCharged = new std::vector<FlowPart>;
680  fVectorPion = new std::vector<FlowPart>;
681  fVectorKaon = new std::vector<FlowPart>;
682  fVectorProton = new std::vector<FlowPart>;
683 
684  fVectorCharged->reserve(300);
685  if(fProcessPID) { fVectorPion->reserve(200); fVectorKaon->reserve(100); fVectorProton->reserve(100); }
686 
687  // creating histograms
688  // event histogram
689 
690  fhEventSampling = new TH2D("fhEventSampling","Event sampling; centrality/multiplicity/impactParameter; sample index", fFlowCentNumBins,0,100, fNumSamples,0,fNumSamples);
691  fQAEvents->Add(fhEventSampling);
692  fhEventCentrality = new TH1D("fhEventCentrality",Form("Event centrality (%s); centrality/multiplicity/impactParameter",fMultEstimator.Data()), fFlowCentNumBins,0,100);
693  fQAEvents->Add(fhEventCentrality);
694  fh2EventCentralityNumSelCharged = new TH2D("fh2EventCentralityNumSelCharged",Form("Event centrality (%s) or impact parameter (b) vs. N^{sel}_{ch}; N^{sel}_{ch}; centrality/multiplicity/impactParameter",fMultEstimator.Data()), 3000,0,3000, fFlowCentNumBins,0,100);
695  fQAEvents->Add(fh2EventCentralityNumSelCharged);
696 
697  const Short_t iEventCounterBins = 11;//10 it was
698  TString sEventCounterLabel[iEventCounterBins] = {"Input","Physics selection OK","Centr. Est. Consis. OK","PV OK","SPD Vtx OK","Pileup MV OK","Out-of-bunch Pileup OK","Vtx Consis. OK","PV #it{z} OK","ESD TPC Mult. Diff. OK","Selected"};
699  fhEventCounter = new TH1D("fhEventCounter","Event Counter",iEventCounterBins,0,iEventCounterBins);
700  for(Short_t i(0); i < iEventCounterBins; i++) fhEventCounter->GetXaxis()->SetBinLabel(i+1, sEventCounterLabel[i].Data() );
701  fQAEvents->Add(fhEventCounter);
702 
703  // flow histograms & profiles
704  // weights
705  if(fAnalType == kAOD && (fFlowFillWeights || fRunMode == kFillWeights))
706  {
707  fh3BeforeNUAWeightsRefs = new TH3D("fh3BeforeNUAWeightsRefs","Weights: Refs; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5, 20, -10,10);
708  fh3BeforeNUAWeightsRefs->Sumw2();
709  fFlowWeights->Add(fh3BeforeNUAWeightsRefs);
710  fh3BeforeNUAWeightsCharged = new TH3D("fh3BeforeNUAWeightsCharged","Weights: Charged; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
711  fh3BeforeNUAWeightsCharged->Sumw2();
712  fFlowWeights->Add(fh3BeforeNUAWeightsCharged);
713  fh3BeforeNUAWeightsPion = new TH3D("fh3BeforeNUAWeightsPion","Weights: #pi; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
714  fh3BeforeNUAWeightsPion->Sumw2();
715  fFlowWeights->Add(fh3BeforeNUAWeightsPion);
716  fh3BeforeNUAWeightsKaon = new TH3D("fh3BeforeNUAWeightsKaon","Weights: K; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
717  fh3BeforeNUAWeightsKaon->Sumw2();
718  fFlowWeights->Add(fh3BeforeNUAWeightsKaon);
719  fh3BeforeNUAWeightsProton = new TH3D("fh3BeforeNUAWeightsProton","Weights: p; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
720  fh3BeforeNUAWeightsProton->Sumw2();
721  fFlowWeights->Add(fh3BeforeNUAWeightsProton);
722 
723  fhBeforeNUEWeightsRefs = new TH1D("fhBeforeNUEWeightsRefs","Weights: Refs; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
724  fhBeforeNUEWeightsRefs->Sumw2();
725  fFlowWeights->Add(fhBeforeNUEWeightsRefs);
726  fhBeforeNUEWeightsCharged = new TH1D("fhBeforeNUEWeightsCharged","Weights: Charged; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
727  fhBeforeNUEWeightsCharged->Sumw2();
728  fFlowWeights->Add(fhBeforeNUEWeightsCharged);
729  fhBeforeNUEWeightsPion = new TH1D("fhBeforeNUEWeightsPion","Weights: #pi; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
730  fhBeforeNUEWeightsPion->Sumw2();
731  fFlowWeights->Add(fhBeforeNUEWeightsPion);
732  fhBeforeNUEWeightsKaon = new TH1D("fhBeforeNUEWeightsKaon","Weights: K; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
733  fhBeforeNUEWeightsKaon->Sumw2();
734  fFlowWeights->Add(fhBeforeNUEWeightsKaon);
735  fhBeforeNUEWeightsProton = new TH1D("fhBeforeNUEWeightsProton","Weights: p; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
736  fhBeforeNUEWeightsProton->Sumw2();
737  fFlowWeights->Add(fhBeforeNUEWeightsProton);
738  }//if(fFlowFillWeights || fRunMode == kFillWeights)
739 
740  if(fAnalType == kAOD && fFlowUseNUAWeights)
741  {
742  fh3AfterNUAWeightsRefs = new TH3D("fh3AfterNUAWeightsRefs","Weights: Refs; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5, 20, -10,10);
743  fh3AfterNUAWeightsRefs->Sumw2();
744  fFlowWeights->Add(fh3AfterNUAWeightsRefs);
745  fh3AfterNUAWeightsCharged = new TH3D("fh3AfterNUAWeightsCharged","Weights: Charged; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
746  fh3AfterNUAWeightsCharged->Sumw2();
747  fFlowWeights->Add(fh3AfterNUAWeightsCharged);
748  fh3AfterNUAWeightsPion = new TH3D("fh3AfterNUAWeightsPion","Weights: #pi; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
749  fh3AfterNUAWeightsPion->Sumw2();
750  fFlowWeights->Add(fh3AfterNUAWeightsPion);
751  fh3AfterNUAWeightsKaon = new TH3D("fh3AfterNUAWeightsKaon","Weights: K; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
752  fh3AfterNUAWeightsKaon->Sumw2();
753  fFlowWeights->Add(fh3AfterNUAWeightsKaon);
754  fh3AfterNUAWeightsProton = new TH3D("fh3AfterNUAWeightsProton","Weights: p; #varphi; #eta; Prim. vtx_{z} (cm)", 100,0,TMath::TwoPi(), 151,-1.5,1.5,20, -10,10);
755  fh3AfterNUAWeightsProton->Sumw2();
756  fFlowWeights->Add(fh3AfterNUAWeightsProton);
757  }//if(fFlowUseNUAWeights)
758  if(fAnalType == kAOD && fFlowUseNUEWeights){
759  fhAfterNUEWeightsRefs = new TH1D("fhAfterNUEWeightsRefs","Weights: Refs; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
760  fhAfterNUEWeightsRefs->Sumw2();
761  fFlowWeights->Add(fhAfterNUEWeightsRefs);
762  fhAfterNUEWeightsCharged = new TH1D("fhAfterNUEWeightsCharged","Weights: Charged; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
763  fhAfterNUEWeightsCharged->Sumw2();
764  fFlowWeights->Add(fhAfterNUEWeightsCharged);
765  fhAfterNUEWeightsPion = new TH1D("fhAfterNUEWeightsPion","Weights: #pi; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
766  fhAfterNUEWeightsPion->Sumw2();
767  fFlowWeights->Add(fhAfterNUEWeightsPion);
768  fhAfterNUEWeightsKaon = new TH1D("fhAfterNUEWeightsKaon","Weights: K; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
769  fhAfterNUEWeightsKaon->Sumw2();
770  fFlowWeights->Add(fhAfterNUEWeightsKaon);
771  fhAfterNUEWeightsProton = new TH1D("fhAfterNUEWeightsProton","Weights: p; #it{p}_{T} (GeV/#it{c})",fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
772  fhAfterNUEWeightsProton->Sumw2();
773  fFlowWeights->Add(fhAfterNUEWeightsProton);
774  }//if(fFlowUseNUEWeights)
775 
776 
777  //mixed correlations
778  //
779  if(fDoOnlyMixedCorrelations){
780  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics; iMixedHarm++)
781  {
782  for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)//For now only for nonoverlapping subevents...
783  {
784  for(Short_t iSample(0); iSample < fNumSamples; iSample++)
785  {
786  //reference flow for mixed harmonics 422,633,523
787  if(iMixedHarm!=3){
788  fpMixedRefsCor4[iSample][iGap][iMixedHarm] = new TProfile(Form("fpRefs_<4>_MixedHarm%d_gap%02.2g_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Ref: <<4>> | Gap %g | v%d | sample %d ; centrality/multiplicity;",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax);
789  fpMixedRefsCor4[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
790  fFlowRefs->Add(fpMixedRefsCor4[iSample][iGap][iMixedHarm]);
791  }
792  if(iMixedHarm==3){
793  //reference flow for mixed harmonics 6222
794  fpMixedRefsCor6[iSample][iGap] = new TProfile(Form("fpRefs_<6>_MixedHarm%d_gap%02.2g_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Ref: <<6>> | Gap %g | v%d | sample %d ; centrality/multiplicity;",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax);
795  fpMixedRefsCor6[iSample][iGap]->Sumw2(kTRUE);
796  fFlowRefs->Add(fpMixedRefsCor6[iSample][iGap]);
797  }
798 
799  if(fProcessCharged){
800  if(iMixedHarm!=3){
801  fpMixedChargedCor3Pos[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Charged_<3>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Charged: <<3'>> | Gap %g | v%d | POIs pos | sample %d ; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
802  fpMixedChargedCor3Pos[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
803  fFlowCharged->Add(fpMixedChargedCor3Pos[iSample][iGap][iMixedHarm]);
804 
805  fpMixedChargedCor3Neg[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Charged_<3>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Charged: <<3'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
806  fpMixedChargedCor3Neg[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
807  fFlowCharged->Add(fpMixedChargedCor3Neg[iSample][iGap][iMixedHarm]);
808  }
809  if(iMixedHarm==3){
810  fpMixedChargedCor4Pos[iSample][iGap] = new TProfile2D(Form("fp2Charged_<4>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Charged: <<4'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
811  fpMixedChargedCor4Pos[iSample][iGap]->Sumw2(kTRUE);
812  fFlowCharged->Add(fpMixedChargedCor4Pos[iSample][iGap]);
813 
814  fpMixedChargedCor4Neg[iSample][iGap] = new TProfile2D(Form("fp2Charged_<4>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Charged: <<4'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
815  fpMixedChargedCor4Neg[iSample][iGap]->Sumw2(kTRUE);
816  fFlowCharged->Add(fpMixedChargedCor4Neg[iSample][iGap]);
817  }
818  }//if(fProcessCharged)
819  if(fProcessPID){
820  if(iMixedHarm!=3){
821  fpMixedPionCor3Pos[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Pion_<3>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Pion: <<3'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
822  fpMixedPionCor3Pos[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
823  fFlowPID->Add(fpMixedPionCor3Pos[iSample][iGap][iMixedHarm]);
824 
825  fpMixedPionCor3Neg[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Pion_<3>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Pion: <<3'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
826  fpMixedPionCor3Neg[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
827  fFlowPID->Add(fpMixedPionCor3Neg[iSample][iGap][iMixedHarm]);
828 
829  fpMixedKaonCor3Pos[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Kaon_<3>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Pion: <<3'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
830  fpMixedKaonCor3Pos[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
831  fFlowPID->Add(fpMixedKaonCor3Pos[iSample][iGap][iMixedHarm]);
832 
833  fpMixedKaonCor3Neg[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Kaon_<3>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Kaon: <<3'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
834  fpMixedKaonCor3Neg[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
835  fFlowPID->Add(fpMixedKaonCor3Neg[iSample][iGap][iMixedHarm]);
836 
837  fpMixedProtonCor3Pos[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Proton_<3>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Proton: <<3'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
838  fpMixedProtonCor3Pos[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
839  fFlowPID->Add(fpMixedProtonCor3Pos[iSample][iGap][iMixedHarm]);
840 
841  fpMixedProtonCor3Neg[iSample][iGap][iMixedHarm] = new TProfile2D(Form("fp2Proton_<3>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Proton: <<3'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
842  fpMixedProtonCor3Neg[iSample][iGap][iMixedHarm]->Sumw2(kTRUE);
843  fFlowPID->Add(fpMixedProtonCor3Neg[iSample][iGap][iMixedHarm]);
844  }
845  if(iMixedHarm==3){
846  fpMixedPionCor4Pos[iSample][iGap] = new TProfile2D(Form("fp2Pion_<4>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Pion: <<4'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
847  fpMixedPionCor4Pos[iSample][iGap]->Sumw2(kTRUE);
848  fFlowPID->Add(fpMixedPionCor4Pos[iSample][iGap]);
849 
850  fpMixedPionCor4Neg[iSample][iGap] = new TProfile2D(Form("fp2Pion_<4>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Pion: <<4'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
851  fpMixedPionCor4Neg[iSample][iGap]->Sumw2(kTRUE);
852  fFlowPID->Add(fpMixedPionCor4Neg[iSample][iGap]);
853 
854  fpMixedKaonCor4Pos[iSample][iGap] = new TProfile2D(Form("fp2Kaon_<4>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Kaon: <<4'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
855  fpMixedKaonCor4Pos[iSample][iGap]->Sumw2(kTRUE);
856  fFlowPID->Add(fpMixedKaonCor4Pos[iSample][iGap]);
857 
858  fpMixedKaonCor4Neg[iSample][iGap] = new TProfile2D(Form("fp2Kaon_<4>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Kaon: <<4'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
859  fpMixedKaonCor4Neg[iSample][iGap]->Sumw2(kTRUE);
860  fFlowPID->Add(fpMixedKaonCor4Neg[iSample][iGap]);
861 
862  fpMixedProtonCor4Pos[iSample][iGap] = new TProfile2D(Form("fp2Proton_<4>_MixedHarm%d_gap%02.2g_Pos_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Proton: <<4'>> | Gap %g | v%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
863  fpMixedProtonCor4Pos[iSample][iGap]->Sumw2(kTRUE);
864  fFlowPID->Add(fpMixedProtonCor4Pos[iSample][iGap]);
865 
866  fpMixedProtonCor4Neg[iSample][iGap] = new TProfile2D(Form("fp2Proton_<4>_MixedHarm%d_gap%02.2g_Neg_sample%d",fMixedHarmonics[iMixedHarm],10*fEtaGap[iGap],iSample),Form("Proton: <<4'>> | Gap %g | v%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fMixedHarmonics[iMixedHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
867  fpMixedProtonCor4Neg[iSample][iGap]->Sumw2(kTRUE);
868  fFlowPID->Add(fpMixedProtonCor4Neg[iSample][iGap]);
869  }
870  }//if(fProcessPID)
871  }//for(Short_t iSample(0); iSample < fNumSamples; iSample++)
872  }//for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
873  }//for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics; iMixedHarm++)
874  }//if(fDoOnlyMixedCorrelations)
875 
876  if(!fDoOnlyMixedCorrelations){
877  // correlations
878  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
879  {
880  for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
881  {
882  for(Short_t iSample(0); iSample < fNumSamples; iSample++)
883  {
884  fpRefsCor2[iSample][iGap][iHarm] = new TProfile(Form("fpRefs_<2>_harm%d_gap%02.2g_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("Ref: <<2>> | Gap %g | n=%d | sample %d; centrality/multiplicity;",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax);
885  fpRefsCor2[iSample][iGap][iHarm]->Sumw2(kTRUE);
886  fFlowRefs->Add(fpRefsCor2[iSample][iGap][iHarm]);
887 
888  if(fProcessCharged)
889  {
890  fp2ChargedCor2Pos[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Charged_<2>_harm%d_gap%02.2g_Pos_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("Charged: <<2'>> | Gap %g | n=%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
891  fp2ChargedCor2Pos[iSample][iGap][iHarm]->Sumw2(kTRUE);
892  fFlowCharged->Add(fp2ChargedCor2Pos[iSample][iGap][iHarm]);
893 
894 
895  fp2ChargedCor2Neg[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Charged_<2>_harm%d_gap%02.2g_Neg_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("Charged: <<2'>> | Gap %g | n=%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
896  fp2ChargedCor2Neg[iSample][iGap][iHarm]->Sumw2(kTRUE);
897  fFlowCharged->Add(fp2ChargedCor2Neg[iSample][iGap][iHarm]);
898 
899  }//if(fProcessCharged)
900 
901  if(fProcessPID)
902  {
903  fp2PionCor2Pos[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Pion_<2>_harm%d_gap%02.2g_Pos_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("PID #pi: <<2'>> | Gap %g | n=%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
904  fp2PionCor2Pos[iSample][iGap][iHarm]->Sumw2(kTRUE);
905  fFlowPID->Add(fp2PionCor2Pos[iSample][iGap][iHarm]);
906 
907  fp2KaonCor2Pos[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Kaon_<2>_harm%d_gap%02.2g_Pos_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("PID K: <<2'>> | Gap %g | n=%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
908  fp2KaonCor2Pos[iSample][iGap][iHarm]->Sumw2(kTRUE);
909  fFlowPID->Add(fp2KaonCor2Pos[iSample][iGap][iHarm]);
910 
911  fp2ProtonCor2Pos[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Proton_<2>_harm%d_gap%02.2g_Pos_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("PID p: <<2'>> | Gap %g | n=%d | POIs pos | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
912  fp2ProtonCor2Pos[iSample][iGap][iHarm]->Sumw2(kTRUE);
913  fFlowPID->Add(fp2ProtonCor2Pos[iSample][iGap][iHarm]);
914 
915  fp2PionCor2Neg[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Pion_<2>_harm%d_gap%02.2g_Neg_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("PID #pi: <<2'>> | Gap %g | n=%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
916  fp2PionCor2Neg[iSample][iGap][iHarm]->Sumw2(kTRUE);
917  fFlowPID->Add(fp2PionCor2Neg[iSample][iGap][iHarm]);
918 
919  fp2KaonCor2Neg[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Kaon_<2>_harm%d_gap%02.2g_Neg_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("PID K: <<2'>> | Gap %g | n=%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
920  fp2KaonCor2Neg[iSample][iGap][iHarm]->Sumw2(kTRUE);
921  fFlowPID->Add(fp2KaonCor2Neg[iSample][iGap][iHarm]);
922 
923  fp2ProtonCor2Neg[iSample][iGap][iHarm] = new TProfile2D(Form("fp2Proton_<2>_harm%d_gap%02.2g_Neg_sample%d",fHarmonics[iHarm],10*fEtaGap[iGap],iSample),Form("PID p: <<2'>> | Gap %g | n=%d | POIs neg | sample %d; centrality/multiplicity; #it{p}_{T} (GeV/c)",fEtaGap[iGap],fHarmonics[iHarm],iSample), fFlowCentNumBins,fFlowCentMin,fFlowCentMax, fFlowPOIsPtNumBins,fFlowPOIsPtMin,fFlowPOIsPtMax);
924  fp2ProtonCor2Neg[iSample][iGap][iHarm]->Sumw2(kTRUE);
925  fFlowPID->Add(fp2ProtonCor2Neg[iSample][iGap][iHarm]);
926  }//endif(fProcessPID)
927  }//for(Short_t iSample(0); iSample < fNumSamples; iSample++)
928  }//for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
929  }//for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
930  }//endif(!fDoOnlyMixedCorrelations)
931 
932  // charged (tracks) histograms
933  fh2RefsMult = new TH2D("fh2RefsMult","RFPs: Centrality: Multiplicity; centrality; multiplicity",fFlowCentNumBins,0,100,1000,0,1000);
934  fQACharged->Add(fh2RefsMult);
935  fh2RefsPt = new TH2D("fh2RefsPt","RFPs: #it{p}_{T}; centrality; #it{p}_{T} (GeV/#it{c})",fFlowCentNumBins,0,100,300,0,30);
936  fQACharged->Add(fh2RefsPt);
937  fh2RefsEta = new TH2D("fh2RefsEta","RFPs: #eta; centrality; #eta",fFlowCentNumBins,0,100, 151,-1.5,1.5);
938  fQACharged->Add(fh2RefsEta);
939  fh2RefsPhi = new TH2D("fh2RefsPhi","RFPs: #varphi; centrality; #varphi",fFlowCentNumBins,0,100,100,0,TMath::TwoPi());
940  fQACharged->Add(fh2RefsPhi);
941 
942  if(fProcessCharged)
943  {
944  TString sChargedCounterLabel[] = {"Input","FB","#TPC-Cls","TPC-Chi2","DCA-z","DCA-xy","Eta","Selected"};
945  const Short_t iNBinsChargedCounter = sizeof(sChargedCounterLabel)/sizeof(sChargedCounterLabel[0]);
946  fhChargedCounter = new TH1D("fhChargedCounter","Charged tracks: Counter",iNBinsChargedCounter,0,iNBinsChargedCounter);
947  for(Short_t i(0); i < iNBinsChargedCounter; i++) fhChargedCounter->GetXaxis()->SetBinLabel(i+1, sChargedCounterLabel[i].Data() );
948  fQACharged->Add(fhChargedCounter);
949  } // endif {fProcessCharged}
950 
951  // PID tracks histograms
952  if(fProcessPID)
953  {
954  fh2PIDPionMult = new TH2D("fh2PIDPionMult","PID: #pi: Multiplicity; centrality; multiplicity", fFlowCentNumBins,0,100,200,0,200);
955  fQAPID->Add(fh2PIDPionMult);
956  fh2PIDKaonMult = new TH2D("fh2PIDKaonMult","PID: K: Multiplicity; centrality; multiplicity", fFlowCentNumBins,0,100,100,0,100);
957  fQAPID->Add(fh2PIDKaonMult);
958  fh2PIDProtonMult = new TH2D("fh2PIDProtonMult","PID: p: Multiplicity; centrality; multiplicity", fFlowCentNumBins,0,100,100,0,100);
959  fQAPID->Add(fh2PIDProtonMult);
960 
961  if(fFillQA)
962  {
963  fh2PIDPionPt = new TH2D("fh2PIDPionPt","PID: #pi: centrality vs. #it{p}_{T}; centrality; #it{p}_{T}", fFlowCentNumBins,0,100,150,0.,30.);
964  fQAPID->Add(fh2PIDPionPt);
965  fh2PIDPionPhi = new TH2D("fh2PIDPionPhi","PID: #pi: centrality vs. #varphi; centrality; #varphi", fFlowCentNumBins,0,100,100,0,TMath::TwoPi());
966  fQAPID->Add(fh2PIDPionPhi);
967  fh2PIDPionEta = new TH2D("fh2PIDPionEta","PID: #pi: centrality vs. #eta; centrality; #eta", fFlowCentNumBins,0,100,151,-1.5,1.5);
968  fQAPID->Add(fh2PIDPionEta);
969  fhPIDPionCharge = new TH1D("fhPIDPionCharge","PID: #pi: charge; charge", 3,-1.5,1.5);
970  fQAPID->Add(fhPIDPionCharge);
971  fh2PIDKaonPt = new TH2D("fh2PIDKaonPt","PID: K: centrality vs. #it{p}_{T}; centrality; #it{p}_{T}", fFlowCentNumBins,0,100,150,0.,30.);
972  fQAPID->Add(fh2PIDKaonPt);
973  fh2PIDKaonPhi = new TH2D("fh2PIDKaonPhi","PID: K: centrality vs. #varphi; centrality; #varphi", fFlowCentNumBins,0,100,100,0,TMath::TwoPi());
974  fQAPID->Add(fh2PIDKaonPhi);
975  fh2PIDKaonEta = new TH2D("fh2PIDKaonEta","PID: K: centrality vs. #eta; centrality; #eta", fFlowCentNumBins,0,100,151,-1.5,1.5);
976  fQAPID->Add(fh2PIDKaonEta);
977  fhPIDKaonCharge = new TH1D("fhPIDKaonCharge","PID: K: charge; charge", 3,-1.5,1.5);
978  fQAPID->Add(fhPIDKaonCharge);
979  fh2PIDProtonPt = new TH2D("fh2PIDProtonPt","PID: p: centrality vs. #it{p}_{T}; centrality; #it{p}_{T}", fFlowCentNumBins,0,100,150,0.,30.);
980  fQAPID->Add(fh2PIDProtonPt);
981  fh2PIDProtonPhi = new TH2D("fh2PIDProtonPhi","PID: p: centrality vs. #varphi; centrality; #varphi", fFlowCentNumBins,0,100,100,0,TMath::TwoPi());
982  fQAPID->Add(fh2PIDProtonPhi);
983  fh2PIDProtonEta = new TH2D("fh2PIDProtonEta","PID: p: centrality vs. #eta; centrality; #eta", fFlowCentNumBins,0,100,151,-1.5,1.5);
984  fQAPID->Add(fh2PIDProtonEta);
985  fhPIDProtonCharge = new TH1D("fhPIDProtonCharge","PID: p: charge; charge", 3,-1.5,1.5);
986  fQAPID->Add(fhPIDProtonCharge);
987  fh2PIDPionTPCdEdx = new TH2D("fh2PIDPionTPCdEdx","PID: #pi: TPC dE/dx; #it{p} (GeV/#it{c}); TPC dE/dx", 200,0,20, 131,-10,1000);
988  fQAPID->Add(fh2PIDPionTPCdEdx);
989  fh2PIDPionTOFbeta = new TH2D("fh2PIDPionTOFbeta","PID: #pi: TOF #beta; #it{p} (GeV/#it{c});TOF #beta", 200,0,20, 101,-0.1,1.5);
990  fQAPID->Add(fh2PIDPionTOFbeta);
991  fh2PIDKaonTPCdEdx = new TH2D("fh2PIDKaonTPCdEdx","PID: K: TPC dE/dx; #it{p} (GeV/#it{c}); TPC dE/dx", 200,0,20, 131,-10,1000);
992  fQAPID->Add(fh2PIDKaonTPCdEdx);
993  fh2PIDKaonTOFbeta = new TH2D("fh2PIDKaonTOFbeta","PID: K: TOF #beta; #it{p} (GeV/#it{c});TOF #beta", 200,0,20, 101,-0.1,1.5);
994  fQAPID->Add(fh2PIDKaonTOFbeta);
995  fh2PIDProtonTPCdEdx = new TH2D("fh2PIDProtonTPCdEdx","PID: p: TPC dE/dx; #it{p} (GeV/#it{c}); TPC dE/dx", 200,0,20, 131,-10,1000);
996  fQAPID->Add(fh2PIDProtonTPCdEdx);
997  fh2PIDProtonTOFbeta = new TH2D("fh2PIDProtonTOFbeta","PID: p: TOF #beta; #it{p} (GeV/#it{c});TOF #beta", 200,0,20, 101,-0.1,1.5);
998  fQAPID->Add(fh2PIDProtonTOFbeta);
999  fh2PIDPionTPCnSigmaPion = new TH2D("fh2PIDPionTPCnSigmaPion","PID: #pi: TPC n#sigma (#pi hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1000  fQAPID->Add(fh2PIDPionTPCnSigmaPion);
1001  fh2PIDPionTOFnSigmaPion = new TH2D("fh2PIDPionTOFnSigmaPion","PID: #pi: TOF n#sigma (#pi hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1002  fQAPID->Add(fh2PIDPionTOFnSigmaPion);
1003  fh2PIDPionTPCnSigmaKaon = new TH2D("fh2PIDPionTPCnSigmaKaon","PID: #pi: TPC n#sigma (K hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1004  fQAPID->Add(fh2PIDPionTPCnSigmaKaon);
1005  fh2PIDPionTOFnSigmaKaon = new TH2D("fh2PIDPionTOFnSigmaKaon","PID: #pi: TOF n#sigma (K hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1006  fQAPID->Add(fh2PIDPionTOFnSigmaKaon);
1007  fh2PIDPionTPCnSigmaProton = new TH2D("fh2PIDPionTPCnSigmaProton","PID: #pi: TPC n#sigma (p hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1008  fQAPID->Add(fh2PIDPionTPCnSigmaProton);
1009  fh2PIDPionTOFnSigmaProton = new TH2D("fh2PIDPionTOFnSigmaProton","PID: #pi: TOF n#sigma (p hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1010  fQAPID->Add(fh2PIDPionTOFnSigmaProton);
1011  fh2PIDKaonTPCnSigmaPion = new TH2D("fh2PIDKaonTPCnSigmaPion","PID: K: TPC n#sigma (#pi hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1012  fQAPID->Add(fh2PIDKaonTPCnSigmaPion);
1013  fh2PIDKaonTOFnSigmaPion = new TH2D("fh2PIDKaonTOFnSigmaPion","PID: K: TOF n#sigma (#pi hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1014  fQAPID->Add(fh2PIDKaonTOFnSigmaPion);
1015  fh2PIDKaonTPCnSigmaKaon = new TH2D("fh2PIDKaonTPCnSigmaKaon","PID: K: TPC n#sigma (K hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1016  fQAPID->Add(fh2PIDKaonTPCnSigmaKaon);
1017  fh2PIDKaonTOFnSigmaKaon = new TH2D("fh2PIDKaonTOFnSigmaKaon","PID: K: TOF n#sigma (K hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1018  fQAPID->Add(fh2PIDKaonTOFnSigmaKaon);
1019  fh2PIDKaonTPCnSigmaProton = new TH2D("fh2PIDKaonTPCnSigmaProton","PID: K: TPC n#sigma (p hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1020  fQAPID->Add(fh2PIDKaonTPCnSigmaProton);
1021  fh2PIDKaonTOFnSigmaProton = new TH2D("fh2PIDKaonTOFnSigmaProton","PID: K: TOF n#sigma (p hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1022  fQAPID->Add(fh2PIDKaonTOFnSigmaProton);
1023  fh2PIDProtonTPCnSigmaPion = new TH2D("fh2PIDProtonTPCnSigmaPion","PID: p: TPC n#sigma (#pi hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1024  fQAPID->Add(fh2PIDProtonTPCnSigmaPion);
1025  fh2PIDProtonTOFnSigmaPion = new TH2D("fh2PIDProtonTOFnSigmaPion","PID: p: TOF n#sigma (#pi hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1026  fQAPID->Add(fh2PIDProtonTOFnSigmaPion);
1027  fh2PIDProtonTPCnSigmaKaon = new TH2D("fh2PIDProtonTPCnSigmaKaon","PID: p: TPC n#sigma (K hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1028  fQAPID->Add(fh2PIDProtonTPCnSigmaKaon);
1029  fh2PIDProtonTOFnSigmaKaon = new TH2D("fh2PIDProtonTOFnSigmaKaon","PID: p: TOF n#sigma (K hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1030  fQAPID->Add(fh2PIDProtonTOFnSigmaKaon);
1031  fh2PIDProtonTPCnSigmaProton = new TH2D("fh2PIDProtonTPCnSigmaProton","PID: p: TPC n#sigma (p hyp.); #it{p}_{T} (GeV/#it{c}); TPC n#sigma", 200,0,20, 42,-11,10);
1032  fQAPID->Add(fh2PIDProtonTPCnSigmaProton);
1033  fh2PIDProtonTOFnSigmaProton = new TH2D("fh2PIDProtonTOFnSigmaProton","PID: p: TOF n#sigma (p hyp.); #it{p}_{T} (GeV/#it{c}); TOF n#sigma", 200,0,20, 42,-11,10);
1034  fQAPID->Add(fh2PIDProtonTOFnSigmaProton);
1035  }
1036 
1037  } //endif {fProcessPID}
1038 
1039  const Short_t iNBinsPIDstatus = 4;
1040  TString sPIDstatus[iNBinsPIDstatus] = {"kDetNoSignal","kDetPidOk","kDetMismatch","kDetNoParams"};
1041  const Short_t iNFilterMapBinBins = 32;
1042 
1043  // QA histograms
1044  if(fFillQA)
1045  {
1046  TString sQAindex[fiNumIndexQA] = {"Before", "After"};
1047  for(Short_t iQA(0); iQA < fiNumIndexQA; iQA++)
1048  {
1049  // EVENTs QA histograms
1050  fhQAEventsPVz[iQA] = new TH1D(Form("fhQAEventsPVz_%s",sQAindex[iQA].Data()), "QA Events: PV-#it{z}", 101,-50,50);
1051  fQAEvents->Add(fhQAEventsPVz[iQA]);
1052  fhQAEventsNumContrPV[iQA] = new TH1D(Form("fhQAEventsNumContrPV_%s",sQAindex[iQA].Data()), "QA Events: Number of contributors to AOD PV", 20,0,20);
1053  fQAEvents->Add(fhQAEventsNumContrPV[iQA]);
1054  fhQAEventsNumSPDContrPV[iQA] = new TH1D(Form("fhQAEventsNumSPDContrPV_%s",sQAindex[iQA].Data()), "QA Events: SPD contributors to PV", 20,0,20);
1055  fQAEvents->Add(fhQAEventsNumSPDContrPV[iQA]);
1056  fhQAEventsDistPVSPD[iQA] = new TH1D(Form("fhQAEventsDistPVSPD_%s",sQAindex[iQA].Data()), "QA Events: PV SPD vertex", 50,0,5);
1057  fQAEvents->Add(fhQAEventsDistPVSPD[iQA]);
1058  fhQAEventsSPDresol[iQA] = new TH1D(Form("fhQAEventsSPDresol_%s",sQAindex[iQA].Data()), "QA Events: SPD resolution", 150,0,15);
1059  fQAEvents->Add(fhQAEventsSPDresol[iQA]);
1060 
1061  fhQAEventsCentralityOutliers[iQA] = new TH2D(Form("fhQAEventsCentralityOutliers_%s",sQAindex[iQA].Data()),"QA Events: Centrality distribution; centrality percentile (V0M);centrality percentile (CL1)",100,0,100,100,0,100);
1062  fQAEvents->Add(fhQAEventsCentralityOutliers[iQA]);
1063 
1064  fhQAEventsPileUp[iQA] = new TH2D(Form("fhQAEventsPileUp_%s",sQAindex[iQA].Data()),"QA Events: TPC vs. ESD multiplicity; TPC multiplicity; ESD multiplicity",500,0,6000,500,0,6000);
1065  fQAEvents->Add(fhQAEventsPileUp[iQA]);
1066  fhEventsMultTOFFilterbit32[iQA] = new TH2D(Form("fhEventsMultTOFFilterbit32_%s",sQAindex[iQA].Data()),"filterbit32 vs. TOF multiplicity; multiplicity(fb32);multiplicity(fb32+TOF)", 4000,0,4000,2000,0,2000);
1067  fQAEvents->Add(fhEventsMultTOFFilterbit32[iQA]);
1068  // Charged tracks QA
1069  if(fProcessCharged)
1070  {
1071  fhQAChargedMult[iQA] = new TH1D(Form("fhQAChargedMult_%s",sQAindex[iQA].Data()),"QA Charged: Number of Charged in selected events; #it{N}^{Charged}", 1500,0,1500);
1072  fQACharged->Add(fhQAChargedMult[iQA]);
1073  fhQAChargedCharge[iQA] = new TH1D(Form("fhQAChargedCharge_%s",sQAindex[iQA].Data()),"QA Charged: Track charge; charge;", 3,-1.5,1.5);
1074  fQACharged->Add(fhQAChargedCharge[iQA]);
1075  fhQAChargedPt[iQA] = new TH1D(Form("fhQAChargedPt_%s",sQAindex[iQA].Data()),"QA Charged: Track #it{p}_{T}; #it{p}_{T} (GeV/#it{c})", 300,0.,30.);
1076  fQACharged->Add(fhQAChargedPt[iQA]);
1077  fhQAChargedEta[iQA] = new TH1D(Form("fhQAChargedEta_%s",sQAindex[iQA].Data()),"QA Charged: Track #it{#eta}; #it{#eta}", 151,-1.5,1.5);
1078  fQACharged->Add(fhQAChargedEta[iQA]);
1079  fhQAChargedPhi[iQA] = new TH1D(Form("fhQAChargedPhi_%s",sQAindex[iQA].Data()),"QA Charged: Track #it{#varphi}; #it{#varphi}", 100,0.,TMath::TwoPi());
1080  fQACharged->Add(fhQAChargedPhi[iQA]);
1081  fhQAChargedFilterBit[iQA] = new TH1D(Form("fhQAChargedFilterBit_%s",sQAindex[iQA].Data()), "QA Charged: Filter bit",iNFilterMapBinBins,0,iNFilterMapBinBins);
1082  for(Int_t j = 0x0; j < iNFilterMapBinBins; j++) fhQAChargedFilterBit[iQA]->GetXaxis()->SetBinLabel(j+1, Form("%g",TMath::Power(2,j)));
1083  fQACharged->Add(fhQAChargedFilterBit[iQA]);
1084  fhQAChargedNumTPCcls[iQA] = new TH1D(Form("fhQAChargedNumTPCcls_%s",sQAindex[iQA].Data()),"QA Charged: Track number of TPC clusters; #it{N}^{TPC clusters}", 160,0,160);
1085  fQACharged->Add(fhQAChargedNumTPCcls[iQA]);
1086  fhQAChargedDCAxy[iQA] = new TH1D(Form("fhQAChargedDCAxy_%s",sQAindex[iQA].Data()),"QA Charged: Track DCA-xy; DCA_{#it{xy}} (cm)", 100,0.,10);
1087  fQACharged->Add(fhQAChargedDCAxy[iQA]);
1088  fhQAChargedDCAz[iQA] = new TH1D(Form("fhQAChargedDCAz_%s",sQAindex[iQA].Data()),"QA Charged: Track DCA-z; DCA_{#it{z}} (cm)", 200,-10.,10.);
1089  fQACharged->Add(fhQAChargedDCAz[iQA]);
1090  } // endif {fProcessCharged}
1091 
1092  // PID tracks QA
1093  if(fProcessPID)
1094  {
1095  fhQAPIDTPCstatus[iQA] = new TH1D(Form("fhQAPIDTPCstatus_%s",sQAindex[iQA].Data()),"QA PID: PID status: TPC;", iNBinsPIDstatus,0,iNBinsPIDstatus);
1096  fQAPID->Add(fhQAPIDTPCstatus[iQA]);
1097  fhQAPIDTPCdEdx[iQA] = new TH2D(Form("fhQAPIDTPCdEdx_%s",sQAindex[iQA].Data()),"QA PID: TPC PID information; #it{p} (GeV/#it{c}); TPC dEdx (au)", 100,0,10, 131,-10,1000);
1098  fQAPID->Add(fhQAPIDTPCdEdx[iQA]);
1099  fhQAPIDTOFstatus[iQA] = new TH1D(Form("fhQAPIDTOFstatus_%s",sQAindex[iQA].Data()),"QA PID: PID status: TOF;", iNBinsPIDstatus,0,iNBinsPIDstatus);
1100  fQAPID->Add(fhQAPIDTOFstatus[iQA]);
1101  fhQAPIDTOFbeta[iQA] = new TH2D(Form("fhQAPIDTOFbeta_%s",sQAindex[iQA].Data()),"QA PID: TOF #beta information; #it{p} (GeV/#it{c}); TOF #beta", 100,0,10, 101,-0.1,1.5);
1102  fQAPID->Add(fhQAPIDTOFbeta[iQA]);
1103 
1104  fh3PIDPionTPCTOFnSigmaPion[iQA] = new TH3D(Form("fh3PIDPionTPCTOFnSigmaPion_%s",sQAindex[iQA].Data()),"QA PID: #pi: TPC-TOF n#sigma (#pi hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1105  fQAPID->Add(fh3PIDPionTPCTOFnSigmaPion[iQA]);
1106  fh3PIDPionTPCTOFnSigmaKaon[iQA] = new TH3D(Form("fh3PIDPionTPCTOFnSigmaKaon_%s",sQAindex[iQA].Data()),"QA PID: #pi: TPC-TOF n#sigma (K hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1107  fQAPID->Add(fh3PIDPionTPCTOFnSigmaKaon[iQA]);
1108  fh3PIDPionTPCTOFnSigmaProton[iQA] = new TH3D(Form("fh3PIDPionTPCTOFnSigmaProton_%s",sQAindex[iQA].Data()),"QA PID: #pi: TPC-TOF n#sigma (p hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1109  fQAPID->Add(fh3PIDPionTPCTOFnSigmaProton[iQA]);
1110 
1111  fh3PIDKaonTPCTOFnSigmaPion[iQA] = new TH3D(Form("fh3PIDKaonTPCTOFnSigmaPion_%s",sQAindex[iQA].Data()),"QA PID: K: TPC-TOF n#sigma (#pi hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1112  fQAPID->Add(fh3PIDKaonTPCTOFnSigmaPion[iQA]);
1113  fh3PIDKaonTPCTOFnSigmaKaon[iQA] = new TH3D(Form("fh3PIDKaonTPCTOFnSigmaKaon_%s",sQAindex[iQA].Data()),"QA PID: K: TPC-TOF n#sigma (K hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1114  fQAPID->Add(fh3PIDKaonTPCTOFnSigmaKaon[iQA]);
1115  fh3PIDKaonTPCTOFnSigmaProton[iQA] = new TH3D(Form("fh3PIDKaonTPCTOFnSigmaProton_%s",sQAindex[iQA].Data()),"QA PID: K: TPC-TOF n#sigma (p hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1116  fQAPID->Add(fh3PIDKaonTPCTOFnSigmaProton[iQA]);
1117 
1118  fh3PIDProtonTPCTOFnSigmaPion[iQA] = new TH3D(Form("fh3PIDProtonTPCTOFnSigmaPion_%s",sQAindex[iQA].Data()),"QA PID: p: TPC-TOF n#sigma (#pi hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1119  fQAPID->Add(fh3PIDProtonTPCTOFnSigmaPion[iQA]);
1120  fh3PIDProtonTPCTOFnSigmaKaon[iQA] = new TH3D(Form("fh3PIDProtonTPCTOFnSigmaKaon_%s",sQAindex[iQA].Data()),"QA PID: p: TPC-TOF n#sigma (K hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1121  fQAPID->Add(fh3PIDProtonTPCTOFnSigmaKaon[iQA]);
1122  fh3PIDProtonTPCTOFnSigmaProton[iQA] = new TH3D(Form("fh3PIDProtonTPCTOFnSigmaProton_%s",sQAindex[iQA].Data()),"QA PID: p: TPC-TOF n#sigma (p hyp.); n#sigma^{TPC}; n#sigma^{TOF}; #it{p}_{T}", 22,-11,10, 22,-11,10, 100,0,10);
1123  fQAPID->Add(fh3PIDProtonTPCTOFnSigmaProton[iQA]);
1124 
1125 
1126 
1127  for(Int_t j = 0x0; j < iNBinsPIDstatus; j++)
1128  {
1129  fhQAPIDTOFstatus[iQA]->GetXaxis()->SetBinLabel(j+1, sPIDstatus[j].Data() );
1130  fhQAPIDTPCstatus[iQA]->GetXaxis()->SetBinLabel(j+1, sPIDstatus[j].Data() );
1131  }
1132  } // endif {fProcessPID}
1133  }
1134  }
1135 
1136  // posting data (mandatory)
1137  PostData(1, fFlowRefs);
1138  PostData(2, fFlowCharged);
1139  PostData(3, fFlowPID);
1140  PostData(4, fQAEvents);
1141  PostData(5, fQACharged);
1142  PostData(6, fQAPID);
1143  PostData(7, fFlowWeights);
1144 
1145  return;
1146 }//AliAnalysisTaskFlowModesMC::UserCreateOutputObjects()
1147 //_____________________________________________________________________________
1148 void AliAnalysisTaskFlowModesMC::ListParameters()
1149 {
1150  // lists all task parameters
1151  // *************************************************************
1152  printf("\n======= List of parameters ========================================\n");
1153  printf(" -------- Analysis task ---------------------------------------\n");
1154  printf(" fRunMode: (RunMode) %d\n", fRunMode);
1155  printf(" fAnalType: (AnalType) %d\n", fAnalType);
1156  printf(" fFillQA: (Bool_t) %s\n", fFillQA ? "kTRUE" : "kFALSE");
1157  printf(" fProcessCharged: (Bool_t) %s\n", fProcessCharged ? "kTRUE" : "kFALSE");
1158  printf(" fProcessPID: (Bool_t) %s\n", fProcessPID ? "kTRUE" : "kFALSE");
1159  printf(" -------- Flow related ----------------------------------------\n");
1160  printf(" fCutFlowDoFourCorrelations: (Bool_t) %s\n", fCutFlowDoFourCorrelations ? "kTRUE" : "kFALSE");
1161  printf(" fCutFlowRFPsPtMin: (Float_t) %g (GeV/c)\n", fCutFlowRFPsPtMin);
1162  printf(" fCutFlowRFPsPtMax: (Float_t) %g (GeV/c)\n", fCutFlowRFPsPtMax);
1163  printf(" fFlowPOIsPtMin: (Float_t) %g (GeV/c)\n", fFlowPOIsPtMin);
1164  printf(" fFlowPOIsPtMax: (Float_t) %g (GeV/c)\n", fFlowPOIsPtMax);
1165  printf(" fFlowCentNumBins: (Int_t) %d (GeV/c)\n", fFlowCentNumBins);
1166  printf(" fFlowCentMin: (Int_t) %d (GeV/c)\n", fFlowCentMin);
1167  printf(" fFlowCentMax: (Int_t) %d (GeV/c)\n", fFlowCentMax);
1168  printf(" fFlowUseNUAWeights: (Bool_t) %s\n", fFlowUseNUAWeights ? "kTRUE" : "kFALSE");
1169  printf(" fFlowUseNUEWeights: (Bool_t) %s\n", fFlowUseNUEWeights ? "kTRUE" : "kFALSE");
1170  printf(" fPositivelyChargedRef: (Bool_t) %s\n", fPositivelyChargedRef ? "kTRUE" : "kFALSE");
1171  printf(" fNegativelyChargedRef: (Bool_t) %s\n", fNegativelyChargedRef ? "kTRUE" : "kFALSE");
1172  printf(" fPositivelyChargedPOI: (Bool_t) %s\n", fPositivelyChargedPOI ? "kTRUE" : "kFALSE");
1173  printf(" fNegativelyChargedPOI: (Bool_t) %s\n", fNegativelyChargedPOI ? "kTRUE" : "kFALSE");
1174  printf(" fFlowNUAWeightsPath: (TString) '%s' \n", fFlowNUAWeightsPath.Data());
1175  printf(" fFlowNUEWeightsPath: (TString) '%s' \n", fFlowNUEWeightsPath.Data());
1176  printf(" -------- Events ----------------------------------------------\n");
1177  printf(" fTrigger: (Short_t) %d\n", fTrigger);
1178  printf(" fMultEstimator: (TString) '%s'\n", fMultEstimator.Data());
1179  printf(" fPVtxCutZ: (Double_t) %g (cm)\n", fPVtxCutZ);
1180  printf(" fFullCentralityRange: (Bool_t) runs over %s centrality range \n", fFullCentralityRange ? "0-100%" : "50-100%");
1181  printf(" -------- Charge tracks ---------------------------------------\n");
1182  printf(" fCutChargedTrackFilterBit: (UInt) %d\n", fCutChargedTrackFilterBit);
1183  printf(" fCutChargedNumTPCclsMin: (UShort_t) %d\n", fCutChargedNumTPCclsMin);
1184  printf(" fCutChargedEtaMax: (Float_t) %g\n", fCutChargedEtaMax);
1185  printf(" fCutChargedPtMin: (Float_t) %g (GeV/c)\n", fCutChargedPtMin);
1186  printf(" fCutChargedPtMax: (Float_t) %g (GeV/c)\n", fCutChargedPtMax);
1187  printf(" fCutChargedDCAzMax: (Float_t) %g (cm)\n", fCutChargedDCAzMax);
1188  printf(" fCutChargedDCAxyMax: (Float_t) %g (cm)\n", fCutChargedDCAxyMax);
1189  printf(" -------- PID (pi,K,p) tracks ---------------------------------\n");
1190  printf(" fCutPIDUseAntiProtonOnly: (Bool_t) %s\n", fCutPIDUseAntiProtonOnly ? "kTRUE" : "kFALSE");
1191  printf(" fCutPIDnSigmaCombinedNoTOFrejection: (Bool_t) %s\n", fCutPIDnSigmaCombinedNoTOFrejection ? "kTRUE" : "kFALSE");
1192  printf(" fCutPIDnSigmaPionMax: (Double_t) %g\n", fCutPIDnSigmaPionMax);
1193  printf(" fCutPIDnSigmaKaonMax: (Double_t) %g\n", fCutPIDnSigmaKaonMax);
1194  printf(" fCutPIDnSigmaProtonMax: (Double_t) %g\n", fCutPIDnSigmaProtonMax);
1195  printf("=====================================================================\n\n");
1196 
1197  return;
1198 }
1199 //_____________________________________________________________________________
1200 Bool_t AliAnalysisTaskFlowModesMC::InitializeTask()
1201 {
1202  // called once on beginning of task (within UserCreateOutputObjects method)
1203  // check if task parameters are specified and valid
1204  // returns kTRUE if succesfull
1205  // *************************************************************
1206 
1207  printf("====== InitializeTask AliAnalysisTaskFlowModesMC =========================\n");
1208  if(fAnalType != kESD && fAnalType != kAOD && fAnalType != kMC && fAnalType != kMCAOD)
1209  {
1210  ::Error("InitializeTask","Analysis type not specified! Terminating!");
1211  return kFALSE;
1212  }
1213 
1214  if(fAnalType == kESD)
1215  {
1216  ::Error("InitializeTask","Analysis type: ESD not implemented! Terminating!");
1217  return kFALSE;
1218  }
1219 
1220  if(fAnalType == kMCAOD)
1221  {
1222  ::Error("InitializeTask","Analysis type: MCAOD not implemented! Terminating!");
1223  return kFALSE;
1224  }
1225 
1226  if(fColSystem != kPP && fColSystem != kPbPb)
1227  {
1228  ::Error("InitializeTask","Collisional system not specified! Terminating!");
1229  return kFALSE;
1230  }
1231  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1232  AliInputEventHandler* inputHandler = (AliInputEventHandler*)mgr->GetInputEventHandler();
1233  // checking PID response
1234  if(fAnalType == kAOD){ //Add later kMCAOD
1235  fPIDResponse = inputHandler->GetPIDResponse();
1236  if(!fPIDResponse)
1237  {
1238  ::Error("InitializeTask","AliPIDResponse object not found! Terminating!");
1239  return kFALSE;
1240  }
1241  fPIDCombined = new AliPIDCombined();
1242  if(!fPIDCombined)
1243  {
1244  ::Error("InitializeTask","AliPIDCombined object not found! Terminating!");
1245  return kFALSE;
1246  }
1247  fPIDCombined->SetDefaultTPCPriors();
1248  fPIDCombined->SetSelectedSpecies(5); // all particle species
1249  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF); // setting TPC + TOF mask
1250  // checking cut setting
1251  ::Info("InitializeTask","Checking task parameters setting conflicts (ranges, etc)");
1252 
1253  // upper-case for multiplicity estimator
1254  fMultEstimator.ToUpper();
1255 
1256  // checking for weights source file
1257  if(fFlowUseNUAWeights && !fFlowNUAWeightsPath.EqualTo(""))
1258  {
1259  fFlowNUAWeightsFile = TFile::Open(Form("alien:///%s",fFlowNUAWeightsPath.Data()));
1260  if(!fFlowNUAWeightsFile)
1261  {
1262  ::Error("InitializeTask","NUA flow weights file not found");
1263  return kFALSE;
1264  }
1265  }
1266  if(fFlowUseNUEWeights && !fFlowNUEWeightsPath.EqualTo(""))
1267  {
1268  fFlowNUEWeightsFile = TFile::Open(Form("alien:///%s",fFlowNUEWeightsPath.Data()));
1269  if(!fFlowNUEWeightsFile)
1270  {
1271  ::Error("InitializeTask","NUE flow weights file not found");
1272  return kFALSE;
1273  }
1274  }
1275  }
1276 
1277 
1278 
1279 
1280  if(fCutFlowRFPsPtMin > 0. && fCutFlowRFPsPtMax > 0. && fCutFlowRFPsPtMin > fCutFlowRFPsPtMax)
1281  {
1282  ::Error("InitializeTask","Cut: RFPs Pt range wrong!");
1283  return kFALSE;
1284  }
1285 
1286  ::Info("InitializeTask","Initialization succesfull!");
1287  printf("======================================================================\n\n");
1288  return kTRUE;
1289 }
1290 //_____________________________________________________________________________
1291 void AliAnalysisTaskFlowModesMC::UserExec(Option_t *)
1292 {
1293  // main method called for each event (event loop)
1294  // *************************************************************
1295 
1296  if(!fInit) return; // check if initialization succesfull
1297 
1298  // local event counter check: if running in test mode, it runs until the 50 events are succesfully processed
1299  if(fRunMode == kTest && fEventCounter >= fNumEventsAnalyse) return;
1300  // event selection
1301  if(fAnalType == kAOD) fEventAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1302 
1303  if(fAnalType == kMC){fEventMC = dynamic_cast<AliMCEvent*>(MCEvent());}
1304 
1305  if(!EventSelection()) return;
1306 
1307  // processing of selected event
1308  if(!ProcessEvent()) return;
1309 
1310  // posting data (mandatory)
1311  PostData(1, fFlowRefs);
1312  PostData(2, fFlowCharged);
1313  PostData(3, fFlowPID);
1314  PostData(4, fQAEvents);
1315  PostData(5, fQACharged);
1316  PostData(6, fQAPID);
1317  PostData(7, fFlowWeights);
1318 
1319  return;
1320 }
1321 //_____________________________________________________________________________
1322 Bool_t AliAnalysisTaskFlowModesMC::EventSelection()
1323 {
1324  // main (envelope) method for event selection
1325  // Specific event selection methods are called from here
1326  // returns kTRUE if event pass all selection criteria
1327  // *************************************************************
1328 
1329  Bool_t eventSelected = kFALSE;
1330 
1331  if(fAnalType == kAOD && !fEventAOD) return kFALSE;
1332  if(fAnalType == kMC && !fEventMC) return kFALSE;
1333  // Fill event QA BEFORE cuts
1334  if(fAnalType == kAOD && fFillQA) FillEventsQA(0);
1335  if(fAnalType == kMC){
1336  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1337  AliInputEventHandler* inputHandler = (AliInputEventHandler*) mgr->GetInputEventHandler();
1338  eventSelected = kTRUE;
1339  //if(fFillQA) FillEventsQA(1);
1340  }
1341  if(fAnalType == kAOD){
1342  // event selection for small systems pp in Run2
1343  if(fColSystem == kPP) eventSelected = IsEventSelected_pp();
1344  // event selection for PbPb in Run2
1345  if(fColSystem == kPbPb) eventSelected = IsEventSelected_PbPb();
1346  }
1347  return eventSelected;
1348 }
1349 //_____________________________________________________________________________
1350 Bool_t AliAnalysisTaskFlowModesMC::IsEventSelected_PbPb()
1351 {
1352  // Event selection for PbPb collisions recorded in Run 2 year 2015
1353  // PbPb (LHC15o)
1354  // return kTRUE if event passes all criteria, kFALSE otherwise
1355  // *************************************************************
1356  fhEventCounter->Fill("Input",1);
1357 
1358  // Physics selection (trigger)
1359  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1360  AliInputEventHandler* inputHandler = (AliInputEventHandler*) mgr->GetInputEventHandler();
1361  UInt_t fSelectMask = inputHandler->IsEventSelected();
1362 
1363  Bool_t isTriggerSelected = kFALSE;
1364  switch(fTrigger) // check for high multiplicity trigger
1365  {
1366  case 0:
1367  isTriggerSelected = fSelectMask& AliVEvent::kINT7;
1368  break;
1369 
1370  case 1:
1371  isTriggerSelected = fSelectMask& AliVEvent::kHighMultV0;
1372  break;
1373 
1374  case 2:
1375  isTriggerSelected = fSelectMask& AliVEvent::kHighMultSPD;
1376  break;
1377 
1378  default: isTriggerSelected = kFALSE;
1379  }
1380  if(!isTriggerSelected){ return kFALSE;}
1381 
1382  // events passing physics selection
1383  fhEventCounter->Fill("Physics selection OK",1);
1384 
1385  // get centrality from AliMultSelection
1386 
1387  AliMultSelection* fMultSelection = 0x0;
1388  fMultSelection = (AliMultSelection*) fEventAOD->FindListObject("MultSelection");
1389  Double_t centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1390  Double_t centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
1391 
1392  // cut on consistency between centrality estimators: VOM vs CL1
1393  if(fabs(centrV0M-centrCL1)>7.5) return kFALSE;
1394  fhEventCounter->Fill("Centr. Est. Consis. OK",1);
1395 
1396 
1397  // primary vertex selection
1398  const AliAODVertex* vtx = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertex());
1399  if(!vtx || vtx->GetNContributors() < 1){ return kFALSE;}
1400  fhEventCounter->Fill("PV OK",1);
1401 
1402  // SPD vertex selection
1403  const AliAODVertex* vtxSPD = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertexSPD());
1404 
1405 
1406  if (((AliAODHeader*)fEventAOD->GetHeader())->GetRefMultiplicityComb08() < 0) return kFALSE;
1407 
1408  if (fEventAOD->IsIncompleteDAQ()) return kFALSE;
1409 
1410  if (vtx->GetNContributors() < 2 || vtxSPD->GetNContributors()<1) return kFALSE;
1411  double cov[6]={0}; double covSPD[6]={0};
1412  vtx->GetCovarianceMatrix(cov);
1413  vtxSPD->GetCovarianceMatrix(covSPD);
1414 
1415  //Special selections for SPD vertex
1416  double zRes = TMath::Sqrt(covSPD[5]);
1417  double fMaxResol=0.25;
1418  if ( vtxSPD->IsFromVertexerZ() && (zRes>fMaxResol)) return kFALSE;
1419  fhEventCounter->Fill("SPD Vtx OK",1);
1420 
1421 
1422  // check for multi-vertexer pile-up
1423  /*
1424  const int kMinPileUpContrib = 5;
1425  const double kMaxPileUpChi2 = 5.0;
1426  const double kMinWDist = 15;
1427  const AliAODVertex* vtxPileUp = 0;
1428  int nPileUp = 0;
1429 
1430  nPileUp=fEventAOD->GetNumberOfPileupVerticesTracks();
1431 
1432  if (nPileUp) {
1433  if (vtx == vtxSPD) return kTRUE; // there are pile-up vertices but no primary
1434  Int_t bcPrim = vtPrm->GetBC();
1435  for (int iPileUp=0;iPileUp<nPileUp;iPileUp++) {
1436  vtxPileUp = (const AliAODVertex*)fEventAOD->GetPileupVertexTracks(iPileUp);
1437  if (vtxPileUp->GetNContributors() < kMinPileUpContrib) continue;
1438  if (vtxPileUp->GetChi2perNDF() > kMaxPileUpChi2) continue;
1439  int bcPlp = vtxPileUp->GetBC(); ///newly added
1440  if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other Bunch crossing (BC)
1441  double wDst = GetWDist(vtx,vtxPileUp);
1442  if (wDst<kMinWDist) continue;
1443  return kTRUE; // pile-up: well separated vertices
1444  }
1445  }
1446  */
1447  // fhEventCounter->Fill("Pileup MV OK",1);
1449  AliAnalysisUtils utils;
1450  utils.SetMinPlpContribMV(5);
1451  utils.SetMaxPlpChi2MV(5);
1452  utils.SetMinWDistMV(15);
1453  utils.SetCheckPlpFromDifferentBCMV(kTRUE);
1454 
1455  Bool_t isPileupFromMV = utils.IsPileUpMV(fEventAOD);
1456 
1457  if(isPileupFromMV) return kFALSE;
1458  fhEventCounter->Fill("Pileup MV OK",1);
1459 
1460  //Bool_t fRejectOutOfBunchPileUp = kFALSE;
1461  //if(fRejectOutOfBunchPileUp) // out-of-bunch rejection (provided by Christian)
1462  //{
1463  //out-of-bunch
1464  // if (utils.IsOutOfBunchPileUp(fEventAOD)){ return kFALSE; }
1465  //}
1466  fhEventCounter->Fill("Out-of-bunch Pileup OK",1);
1468  // check vertex consistency
1469  double dz = vtx->GetZ() - vtxSPD->GetZ();
1470  double errTot = TMath::Sqrt(cov[5]+covSPD[5]);
1471  double err = TMath::Sqrt(cov[5]);
1472  double nsigTot = dz/errTot;
1473  double nsig = dz/err;
1474  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsig)>20) return kFALSE;
1475  fhEventCounter->Fill("Vtx Consis. OK",1);
1476 
1477 
1478  const Double_t aodVtxZ = vtx->GetZ();
1479  if( TMath::Abs(aodVtxZ) > fPVtxCutZ ){return kFALSE;}
1480  fhEventCounter->Fill("PV #it{z} OK",1);
1481 
1482  // cut on # ESD tracks vs # TPC only tracks
1483  const Int_t nTracks = fEventAOD->GetNumberOfTracks();
1484  Int_t multEsd = ((AliAODHeader*)fEventAOD->GetHeader())->GetNumberOfESDTracks();
1485  Int_t multTrk = 0;
1486  //Int_t multTrkBefC = 0;
1487  Int_t multTrkTOF = 0;
1488  Int_t multTPC = 0;
1489  for (Int_t it = 0; it < nTracks; it++) {
1490  AliAODTrack* AODTrk = (AliAODTrack*)fEventAOD->GetTrack(it);
1491  if (!AODTrk){ delete AODTrk; continue; }
1492  if (AODTrk->TestFilterBit(128)) {multTPC++;}
1493  if (AODTrk->TestFilterBit(32)){
1494  multTrk++;
1495  if ( TMath::Abs(AODTrk->GetTOFsignalDz()) <= 10 && AODTrk->GetTOFsignal() >= 12000 && AODTrk->GetTOFsignal() <= 25000) multTrkTOF++;
1496  }
1497  } // end of for (Int_t it = 0; it < nTracks; it++)
1498  Double_t multTPCn = multTPC;
1499  Double_t multEsdn = multEsd;
1500 
1501  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
1502 
1503  if (multESDTPCDif > 15000.) return kFALSE;
1504 
1505  fhEventCounter->Fill("ESD TPC Mult. Diff. OK",1);
1506 
1507  Double_t multTrkn = multTrk;
1508  Double_t multTrkTOFn = multTrkTOF;
1509 
1510  if(fExtraPileUp && multTrkTOFn< (-32+ 0.32*multTrkn+0.000037*multTrkn*multTrkn)) return kFALSE;
1511  if(fExtraPileUp && multTrkTOFn> (13+0.46*multTrkn+0.000018*multTrkn*multTrkn)) return kFALSE;
1512 
1513  fhEventCounter->Fill("TOF fb32 Mult. correlation OK",1);
1514 
1515  fhEventCounter->Fill("Selected",1);
1516 
1517  // Fill event QA AFTER cuts
1518  if(fFillQA) FillEventsQA(1);
1519 
1520  return kTRUE;
1521 
1522 }
1523 //_____________________________________________________________________________
1524 Bool_t AliAnalysisTaskFlowModesMC::IsEventSelected_pp()
1525 {
1526  // Event selection for small system collision recorder in Run 2 year 2016
1527  // pp (LHC16kl), pPb (LHC16rqts)
1528  // return kTRUE if event passes all criteria, kFALSE otherwise
1529  // *************************************************************
1530 
1531  fhEventCounter->Fill("Input",1);
1532 
1533  // Physics selection (trigger)
1534  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1535  AliInputEventHandler* inputHandler = (AliInputEventHandler*) mgr->GetInputEventHandler();
1536  UInt_t fSelectMask = inputHandler->IsEventSelected();
1537 
1538  Bool_t isTriggerSelected = kFALSE;
1539  switch(fTrigger) // check for high multiplicity trigger
1540  {
1541  case 0:
1542  isTriggerSelected = fSelectMask& AliVEvent::kINT7;
1543  break;
1544 
1545  case 1:
1546  isTriggerSelected = fSelectMask& AliVEvent::kHighMultV0;
1547  break;
1548 
1549  case 2:
1550  isTriggerSelected = fSelectMask& AliVEvent::kHighMultSPD;
1551  break;
1552 
1553  default: isTriggerSelected = kFALSE;
1554  }
1555 
1556  if(!isTriggerSelected)
1557  return kFALSE;
1558 
1559  // events passing physics selection
1560  fhEventCounter->Fill("Physics selection OK",1);
1561 
1562  // primary vertex selection
1563  const AliAODVertex* vtx = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertex());
1564  if(!vtx || vtx->GetNContributors() < 1)
1565  return kFALSE;
1566  fhEventCounter->Fill("PV OK",1);
1567 
1568  // SPD vertex selection
1569  const AliAODVertex* vtxSPD = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertexSPD());
1570 
1571  Double_t dMaxResol = 0.25; // suggested from DPG
1572  Double_t cov[6] = {0};
1573  vtxSPD->GetCovarianceMatrix(cov);
1574  Double_t zRes = TMath::Sqrt(cov[5]);
1575  if ( vtxSPD->IsFromVertexerZ() && (zRes > dMaxResol)) return kFALSE;
1576  fhEventCounter->Fill("SPD Vtx OK",1);
1577 
1578  // PileUp rejection included in Physics selection
1579  // but with values for high mult pp (> 5 contrib) => for low ones: do manually (> 3 contrib)
1580 
1581  /*
1582  if(fTrigger == 0 && fAOD->IsPileupFromSPD(3,0.8) )
1583  {
1584  return kFALSE;
1585  }
1586  */
1587 
1588  //fhEventCounter->Fill("Pileup SPD OK",1);
1589 
1590  // pileup rejection from multivertexer
1591 
1592  AliAnalysisUtils utils;
1593  utils.SetMinPlpContribMV(5);
1594  utils.SetMaxPlpChi2MV(5);
1595  utils.SetMinWDistMV(15);
1596  utils.SetCheckPlpFromDifferentBCMV(kFALSE);
1597  Bool_t isPileupFromMV = utils.IsPileUpMV(fEventAOD);
1598  utils.SetCheckPlpFromDifferentBCMV(kTRUE);
1599 
1600  if(isPileupFromMV) return kFALSE;
1601  fhEventCounter->Fill("Pileup MV OK",1);
1602 
1603  //Bool_t fRejectOutOfBunchPileUp = kTRUE; //you have to add this to the task and set from header
1604  //if(fRejectOutOfBunchPileUp) // out-of-bunch rejection (provided by Christian)
1605  //{
1606  //out-of-bunch
1607  // if (utils.IsOutOfBunchPileUp(fEventAOD)){ return kFALSE; }
1608  //}
1609  fhEventCounter->Fill("Out-of-bunch Pileup OK",1);
1610 
1611  // if (utils.IsSPDClusterVsTrackletBG(fEventAOD))
1612  // {
1613  // return kFALSE;
1614  // }
1615  //
1616  // fhEventCounter->Fill("SPDClTrBG OK",1);
1617  //
1618  // // SPD pileup
1619  // if (utils.IsPileUpSPD(fEventAOD))
1620  // {
1621  // return kFALSE;
1622  // }
1623  //
1624  // fhEventCounter->Fill("SPDPU OK",1);
1625  // }
1626 
1627  //fhEventCounter->Fill("Utils OK",1);
1628 
1629  // cutting on PV z-distance
1630  const Double_t aodVtxZ = vtx->GetZ();
1631  if( TMath::Abs(aodVtxZ) > fPVtxCutZ )
1632  {
1633  return kFALSE;
1634  }
1635  fhEventCounter->Fill("PV #it{z} OK",1);
1636 
1637  fhEventCounter->Fill("Selected",1);
1638  return kTRUE;
1639 }
1640 //_____________________________________________________________________________
1641 void AliAnalysisTaskFlowModesMC::FillEventsQA(const Short_t iQAindex)
1642 {
1643  // Filling various QA plots related with event selection
1644  // *************************************************************
1645 
1646  const AliAODVertex* aodVtx = fEventAOD->GetPrimaryVertex();
1647  const Double_t dVtxZ = aodVtx->GetZ();
1648  const Int_t iNumContr = aodVtx->GetNContributors();
1649  const AliAODVertex* spdVtx = fEventAOD->GetPrimaryVertexSPD();
1650  const Int_t iNumContrSPD = spdVtx->GetNContributors();
1651  const Double_t spdVtxZ = spdVtx->GetZ();
1652 
1653  fhQAEventsPVz[iQAindex]->Fill(dVtxZ);
1654  fhQAEventsNumContrPV[iQAindex]->Fill(iNumContr);
1655  fhQAEventsNumSPDContrPV[iQAindex]->Fill(iNumContrSPD);
1656  fhQAEventsDistPVSPD[iQAindex]->Fill(TMath::Abs(dVtxZ - spdVtxZ));
1657 
1658  // SPD vertexer resolution
1659  Double_t cov[6] = {0};
1660  spdVtx->GetCovarianceMatrix(cov);
1661  Double_t zRes = TMath::Sqrt(cov[5]);
1662  fhQAEventsSPDresol[iQAindex]->Fill(zRes);
1663 
1664  AliMultSelection* fMultSelection = 0x0;
1665  fMultSelection = (AliMultSelection*) fEventAOD->FindListObject("MultSelection");
1666  Double_t centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1667  Double_t centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
1668 
1669  fhQAEventsCentralityOutliers[iQAindex]->Fill(centrV0M,centrCL1);
1670 
1671  const Int_t nTracks = fEventAOD->GetNumberOfTracks();
1672  Int_t multEsd = ((AliAODHeader*)fEventAOD->GetHeader())->GetNumberOfESDTracks();
1673  Int_t multTPC = 0;
1674  Int_t multTrk = 0;
1675  Int_t multTrkTOF = 0;
1676  for (Int_t it = 0; it < nTracks; it++) {
1677  AliAODTrack* AODTrk = (AliAODTrack*)fEventAOD->GetTrack(it);
1678  if (!AODTrk){ delete AODTrk; continue; }
1679  if (AODTrk->TestFilterBit(128)) {multTPC++;}
1680  if (AODTrk->TestFilterBit(32)){
1681  multTrk++;
1682  if ( TMath::Abs(AODTrk->GetTOFsignalDz()) <= 10 && AODTrk->GetTOFsignal() >= 12000 && AODTrk->GetTOFsignal() <= 25000) multTrkTOF++;
1683  }
1684 
1685  }// end of for (Int_t it = 0; it < nTracks; it++)
1686  fhQAEventsPileUp[iQAindex]->Fill(multTPC,multEsd);
1687  fhEventsMultTOFFilterbit32[iQAindex]->Fill(multTrk,multTrkTOF);
1688 
1689  return;
1690 }
1691 //_____________________________________________________________________________
1692 Bool_t AliAnalysisTaskFlowModesMC::Filtering() //void
1693 {
1694 
1695  // main (envelope) method for filtering all particles of interests (POIs) in selected events
1696  // All POIs passing selection criteria are saved to relevant TClonesArray for further processing
1697  // return kTRUE if succesfull (no errors in process)
1698  // *************************************************************
1699  if(!fProcessCharged && !fProcessPID) // if neither is ON, filtering is skipped
1700  return kFALSE; //return;
1701 
1702  fVectorCharged->clear();
1703 
1704  FilterCharged();
1705 
1706  if(fAnalType == kAOD){
1707  // estimate centrality & assign indexes (centrality/percentile, ...)
1708  if(fColSystem == kPbPb){
1709 
1710  fIndexCentrality = GetCentralityIndex();
1711  if(fIndexCentrality < 0) return kFALSE; // return; not succesfull estimation
1712  Double_t Mult = fVectorCharged->size();
1713  if(fExtraPileUp && fIndexCentrality< (-1.5*TMath::Power(Mult,0.46)-0.6*TMath::Log(Mult)*TMath::Log(Mult)+81)) return kFALSE;
1714  if(fExtraPileUp && fIndexCentrality>(-2.3*TMath::Power(Mult,0.39)-0.9*TMath::Log(Mult)*TMath::Log(Mult)+110)) return kFALSE;
1715 
1716  }
1717  if(fColSystem == kPP){fIndexCentrality = 1;}
1718  fh2EventCentralityNumSelCharged->Fill(fVectorCharged->size(),fIndexCentrality);
1719  fhEventCentrality->Fill(fIndexCentrality);
1720  fIndexSampling = GetSamplingIndex();
1721  fhEventSampling->Fill(fIndexCentrality,fIndexSampling);
1722  }
1723  if(fAnalType == kMC){
1724 
1725  if(fColSystem == kPbPb){
1726  fIndexImpactParameter = GetImpactParameter();
1727  if(fIndexImpactParameter < 0) return kFALSE; // return; not succesfull estimation
1728 
1729  }
1730  if(fColSystem == kPP){fIndexImpactParameter = 1;}
1731  fh2EventCentralityNumSelCharged->Fill(fVectorCharged->size(),fIndexImpactParameter);
1732  fhEventCentrality->Fill(fIndexImpactParameter);
1733  fIndexSampling = GetSamplingIndex();
1734  fhEventSampling->Fill(fIndexImpactParameter,fIndexSampling);
1735  }
1736 
1737  if(fProcessPID)
1738  {
1739  fVectorPion->clear();
1740  fVectorKaon->clear();
1741  fVectorProton->clear();
1742  FilterPID();
1743 
1744  }
1745 
1746  return kTRUE; //return;
1747 }
1748 //_____________________________________________________________________________
1749 void AliAnalysisTaskFlowModesMC::FilterCharged()
1750 {
1751  // Filtering input charged tracks
1752  // If track passes all requirements as defined in IsChargedSelected(),
1753  // the relevant properties (pT, eta, phi) are stored in FlowPart struct
1754  // and pushed to relevant vector container.
1755  // return kFALSE if any complications occurs
1756  // *************************************************************
1757  AliAODTrack* track = 0x0;
1758  AliMCParticle* particle = 0x0;
1759  Int_t iNumRefs = 0;
1760 
1761  Double_t NUAweight = 0;
1762  Double_t NUEweight = 0;
1763  const Double_t iNumTracks = (fAnalType == kAOD) ? fEventAOD->GetNumberOfTracks() : fEventMC->GetNumberOfTracks();
1764 
1765  if(iNumTracks < 1) return;
1766  /*
1767  if(fAnalType == kAOD){
1768  if(fIndexCentrality>=0 && fIndexCentrality<10){
1769  Int_t c = fIndexCentrality/5;
1770  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
1771  Int_t c = 1 + fIndexCentrality/10;
1772  }else if(fIndexCentrality>=50){ Int_t c = -1;}
1773  }else{Int_t c=-1;}
1774  */
1775 
1776  for(Double_t iTrack(0); iTrack < iNumTracks; iTrack++)
1777  {
1778  if(fAnalType == kMC){
1779  particle = static_cast<AliMCParticle*>(fEventMC->GetTrack(iTrack));
1780  if(fFillQA) FillQACharged(0,particle);
1781  fVectorCharged->emplace_back( FlowPart(particle->Pt(),particle->Phi(),particle->Eta(), particle->Charge(), kCharged) );
1782  if(fFillQA) FillQACharged(1,particle); // QA after selection
1783 
1784  if(fCutFlowRFPsPtMin > 0. && particle->Pt() >= fCutFlowRFPsPtMin && fCutFlowRFPsPtMax > 0. && particle->Pt() <= fCutFlowRFPsPtMax){
1785  iNumRefs++;
1786  FillQARefs(1,particle);
1787  }
1788  }
1789 
1790  if(fAnalType == kAOD){
1791  track = static_cast<AliAODTrack*>(fEventAOD->GetTrack(iTrack));
1792  if(!track) continue;
1793 
1794  if(fFillQA) FillQACharged(0,track); // QA before selection
1795 
1796  if(fNegativelyChargedRef==kTRUE && track->Charge()>0) continue;
1797  if(fPositivelyChargedRef==kTRUE && track->Charge()<0) continue;
1798 
1799  if(!IsChargedSelected(track)) continue;
1800 
1801  fVectorCharged->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kCharged) );
1802  if(fRunMode == kFillWeights || fFlowFillWeights){
1803  fh3BeforeNUAWeightsCharged->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
1804  fhBeforeNUEWeightsCharged->Fill(track->Pt());
1805  }
1806  if(fFlowUseNUAWeights)
1807  {
1808  if(track->Charge()>0){ NUAweight = fh3NUAWeightChargedPlus->GetBinContent( fh3NUAWeightChargedPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1809  if(track->Charge()<0){ NUAweight = fh3NUAWeightChargedMinus->GetBinContent( fh3NUAWeightChargedMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1810 
1811  fh3AfterNUAWeightsCharged->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
1812  }
1813  if(fFlowUseNUEWeights){
1814  /*
1815  if(fIndexCentrality>=0 && fIndexCentrality<10){
1816  Int_t c = fIndexCentrality/5;
1817  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
1818  Int_t c = 1 + fIndexCentrality/10;
1819  }else if(fIndexCentrality>=50){continue;}
1820  */
1821  if(track->Charge()>0 && fIndexCentrality<50){NUEweight = fhNUEWeightChargedPlus[fIndexCentrality/10]->GetBinContent( fhNUEWeightChargedPlus[fIndexCentrality/10]->FindBin(track->Pt()));}
1822  if(track->Charge()<0 && fIndexCentrality<50){NUEweight = fhNUEWeightChargedMinus[fIndexCentrality/10]->GetBinContent( fhNUEWeightChargedMinus[fIndexCentrality/10]->FindBin(track->Pt()));}
1823 
1824  fhAfterNUEWeightsCharged->Fill(track->Pt(),NUEweight);
1825  }
1826 
1827  if(fFillQA) FillQACharged(1,track); // QA after selection
1828 
1829  // Filling refs QA plots
1830  if(fCutFlowRFPsPtMin > 0. && track->Pt() >= fCutFlowRFPsPtMin && fCutFlowRFPsPtMax > 0. && track->Pt() <= fCutFlowRFPsPtMax)
1831  {
1832  iNumRefs++;
1833  if(fRunMode == kFillWeights || fFlowFillWeights){
1834  fh3BeforeNUAWeightsRefs->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
1835  fhBeforeNUEWeightsRefs->Fill(track->Pt());
1836  }
1837  if(fFlowUseNUAWeights)
1838  {
1839  if(track->Charge()>0){ NUAweight = fh3NUAWeightRefsPlus->GetBinContent( fh3NUAWeightRefsPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1840  if(track->Charge()<0){ NUAweight = fh3NUAWeightRefsMinus->GetBinContent( fh3NUAWeightRefsMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1841 
1842  fh3AfterNUAWeightsRefs->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
1843  }
1844  if(fFlowUseNUEWeights)
1845  {
1846  /*if(fIndexCentrality>=0 && fIndexCentrality<10){
1847  Int_t c = fIndexCentrality/5;
1848  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
1849  Int_t c = 1 + fIndexCentrality/10;
1850  }else if(fIndexCentrality>=50){continue;}
1851  */
1852  if(track->Charge()>0 && fIndexCentrality<50){NUEweight = fhNUEWeightRefsPlus[fIndexCentrality/10]->GetBinContent( fhNUEWeightRefsPlus[fIndexCentrality/10]->FindBin(track->Pt()) );}
1853  if(track->Charge()<0 && fIndexCentrality<50){NUEweight = fhNUEWeightRefsMinus[fIndexCentrality/10]->GetBinContent( fhNUEWeightRefsMinus[fIndexCentrality/10]->FindBin(track->Pt()) );}
1854 
1855  fhAfterNUEWeightsRefs->Fill(track->Pt(),NUEweight);
1856  }
1857  FillQARefs(1,track);
1858  }
1859  }
1860 
1861  }
1862 
1863  // fill QA charged multiplicity
1864  if(fAnalType == kAOD) fh2RefsMult->Fill(fIndexCentrality,iNumRefs);
1865  if(fAnalType == kMC) fh2RefsMult->Fill(fIndexImpactParameter,iNumRefs);
1866 
1867  if(fFillQA)
1868  {
1869  fhQAChargedMult[0]->Fill(iNumTracks);
1870  fhQAChargedMult[1]->Fill(fVectorCharged->size());
1871  }
1872 
1873  return;
1874 }
1875 //_____________________________________________________________________________
1876 Bool_t AliAnalysisTaskFlowModesMC::IsChargedSelected(const AliAODTrack* track)
1877 {
1878  // Selection of charged track
1879  // returns kTRUE if track pass all requirements, kFALSE otherwise
1880  // *************************************************************
1881  if(!track) return kFALSE;
1882  fhChargedCounter->Fill("Input",1);
1883 
1884  // filter bit
1885  if( !track->TestFilterBit(fCutChargedTrackFilterBit) ) return kFALSE;
1886  fhChargedCounter->Fill("FB",1);
1887 
1888  // number of TPC clusters (additional check for not ITS-standalone tracks)
1889  if( track->GetTPCNcls() < fCutChargedNumTPCclsMin && fCutChargedTrackFilterBit != 2) return kFALSE;
1890  fhChargedCounter->Fill("#TPC-Cls",1);
1891 
1892  // track chi2 per space points
1893  Double_t chi2TPC =0.;
1894  chi2TPC = track->Chi2perNDF();
1895  if (fMaxChi2perTPCcls > 0. && chi2TPC > fMaxChi2perTPCcls) return kFALSE;
1896  fhChargedCounter->Fill("TPC-Chi2",1);
1897 
1898  // track DCA coordinates
1899  // note AliAODTrack::XYZAtDCA() works only for constrained tracks
1900  Double_t dTrackXYZ[3] = {0};
1901  Double_t dVertexXYZ[3] = {0.};
1902  Double_t dDCAXYZ[3] = {0.};
1903  if( fCutChargedDCAzMax > 0. || fCutChargedDCAxyMax > 0.)
1904  {
1905  const AliAODVertex* vertex = fEventAOD->GetPrimaryVertex();
1906  if(!vertex) return kFALSE; // event does not have a PV
1907 
1908  track->GetXYZ(dTrackXYZ);
1909  vertex->GetXYZ(dVertexXYZ);
1910 
1911  for(Short_t i(0); i < 3; i++)
1912  dDCAXYZ[i] = dTrackXYZ[i] - dVertexXYZ[i];
1913  }
1914 
1915  if(fCutChargedDCAzMax > 0. && TMath::Abs(dDCAXYZ[2]) > fCutChargedDCAzMax) return kFALSE;
1916  fhChargedCounter->Fill("DCA-z",1);
1917 
1918  if(fCutChargedDCAxyMax > 0. && TMath::Sqrt(dDCAXYZ[0]*dDCAXYZ[0] + dDCAXYZ[1]*dDCAXYZ[1]) > fCutChargedDCAxyMax) return kFALSE;
1919  fhChargedCounter->Fill("DCA-xy",1);
1920 
1921  // pseudorapidity (eta)
1922  if(fCutChargedEtaMax > 0. && TMath::Abs(track->Eta()) > fCutChargedEtaMax) return kFALSE;
1923  fhChargedCounter->Fill("Eta",1);
1924 
1925  // track passing all criteria
1926  fhChargedCounter->Fill("Selected",1);
1927  return kTRUE;
1928 }
1929 //_____________________________________________________________________________
1930 void AliAnalysisTaskFlowModesMC::FillQARefs(const Short_t iQAindex, const AliAODTrack* track)
1931 {
1932  // Filling various QA plots related to RFPs subset of charged track selection
1933  // *************************************************************
1934 
1935  if(!track) return;
1936  if(iQAindex == 0) return; // NOTE implemented only for selected RFPs
1937 
1938  fh2RefsPt->Fill(fIndexCentrality,track->Pt());
1939  fh2RefsEta->Fill(fIndexCentrality,track->Eta());
1940  fh2RefsPhi->Fill(fIndexCentrality,track->Phi());
1941 
1942  return;
1943 }
1944 //_____________________________________________________________________________
1945 void AliAnalysisTaskFlowModesMC::FillQARefs(const Short_t iQAindex, const AliMCParticle* particle)
1946 {
1947  if(!particle) return;
1948  if(iQAindex == 0) return; // NOTE implemented only for selected RFPs
1949 
1950  fh2RefsPt->Fill(fIndexImpactParameter,particle->Pt());
1951  fh2RefsEta->Fill(fIndexImpactParameter,particle->Eta());
1952  fh2RefsPhi->Fill(fIndexImpactParameter,particle->Phi());
1953 
1954  return;
1955 }
1956 //_____________________________________________________________________________
1957 void AliAnalysisTaskFlowModesMC::FillQACharged(const Short_t iQAindex, const AliAODTrack* track)
1958 {
1959  // Filling various QA plots related to charged track selection
1960  // *************************************************************
1961  if(!track) return;
1962 
1963  // filter bit testing
1964  for(Short_t i(0); i < 32; i++)
1965  {
1966  if(track->TestFilterBit(TMath::Power(2.,i)))
1967  fhQAChargedFilterBit[iQAindex]->Fill(i);
1968  }
1969 
1970  // track charge
1971  fhQAChargedCharge[iQAindex]->Fill(track->Charge());
1972 
1973  // number of TPC clusters
1974  fhQAChargedNumTPCcls[iQAindex]->Fill(track->GetTPCNcls());
1975 
1976  // track DCA
1977  Double_t dDCAXYZ[3] = {-999., -999., -999.};
1978  const AliAODVertex* vertex = fEventAOD->GetPrimaryVertex();
1979  if(vertex)
1980  {
1981  Double_t dTrackXYZ[3] = {-999., -999., -999.};
1982  Double_t dVertexXYZ[3] = {-999., -999., -999.};
1983 
1984  track->GetXYZ(dTrackXYZ);
1985  vertex->GetXYZ(dVertexXYZ);
1986 
1987  for(Short_t i(0); i < 3; i++)
1988  dDCAXYZ[i] = dTrackXYZ[i] - dVertexXYZ[i];
1989  }
1990  fhQAChargedDCAxy[iQAindex]->Fill(TMath::Sqrt(dDCAXYZ[0]*dDCAXYZ[0] + dDCAXYZ[1]*dDCAXYZ[1]));
1991  fhQAChargedDCAz[iQAindex]->Fill(dDCAXYZ[2]);
1992 
1993  // kinematics
1994  fhQAChargedPt[iQAindex]->Fill(track->Pt());
1995  fhQAChargedPhi[iQAindex]->Fill(track->Phi());
1996  fhQAChargedEta[iQAindex]->Fill(track->Eta());
1997 
1998  return;
1999 }
2000 //_____________________________________________________________________________
2001 void AliAnalysisTaskFlowModesMC::FillQACharged(const Short_t iQAindex, const AliMCParticle* particle)
2002 {
2003  fhQAChargedPt[iQAindex]->Fill(particle->Pt());
2004  fhQAChargedPhi[iQAindex]->Fill(particle->Phi());
2005  fhQAChargedEta[iQAindex]->Fill(particle->Eta());
2006 
2007  return;
2008 }
2009 //_____________________________________________________________________________
2010 void AliAnalysisTaskFlowModesMC::FilterPID()
2011 {
2012  // If track passes all requirements as defined in IsPIDSelected() (and species dependent),
2013  // the relevant properties (pT, eta, phi) are stored in FlowPart struct
2014  // and pushed to relevant vector container.
2015  // return kFALSE if any complications occurs
2016  // *************************************************************
2017 
2018  const Double_t iNumTracks = (fAnalType == kAOD) ? fEventAOD->GetNumberOfTracks() : fEventMC->GetNumberOfTracks();
2019  if(iNumTracks < 1) return;
2020  //if(fAnalType == kAOD) const Short_t iNumTracks = fEventAOD->GetNumberOfTracks();
2021  //if(fAnalType == kMC) const Short_t iNumTracks = fEventMC->GetNumberOfTracks();
2022 
2023 
2024  PartSpecies species = kUnknown;
2025  AliAODTrack* track = 0x0;
2026  AliMCParticle* particle = 0x0;
2027  Double_t NUAweight = 0;
2028  Double_t NUEweight = 0;
2029  /*
2030  if(fAnalType == kAOD){
2031  if(fIndexCentrality>=0 && fIndexCentrality<10){
2032  Int_t c = fIndexCentrality/5;
2033  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
2034  Int_t c = 1 + fIndexCentrality/10;
2035  }else if(fIndexCentrality>=50){ Int_t c = -1;}
2036  }else{Int_t c=-1;}*/
2037 
2038  for(Double_t iTrack(0); iTrack < iNumTracks; iTrack++)
2039  {
2040  if(fAnalType == kAOD) track = static_cast<AliAODTrack*>(fEventAOD->GetTrack(iTrack));
2041  if(fAnalType == kMC) particle = static_cast<AliMCParticle*>(fEventMC->GetTrack(iTrack));
2042  if(fAnalType == kAOD && !track) continue;
2043  if(fAnalType == kMC && !particle) continue;
2044 
2045  // PID tracks are subset of selected charged tracks (same quality requirements)
2046  if(fAnalType == kAOD){
2047  if(!IsChargedSelected(track)) continue;
2048  if(fFillQA) FillPIDQA(0,track,kUnknown); // filling QA for tracks before selection (but after charged criteria applied)
2049  if(fNegativelyChargedPOI==kTRUE && track->Charge()>0) continue;
2050  if(fPositivelyChargedPOI==kTRUE && track->Charge()<0) continue;
2051  // PID track selection (return most favourable species)
2052  PartSpecies species = IsPIDSelected(track);
2053 
2054  // check if only protons should be used
2055  if(fCutPIDUseAntiProtonOnly && species == kProton && track->Charge() == 1) species = kUnknown;
2056 
2057  // selection of PID tracks
2058  switch (species)
2059  {
2060  case kPion:
2061  fVectorPion->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kPion, fPDGMassPion, track->Px(), track->Py(), track->Pz()) );
2062  if(fRunMode == kFillWeights || fFlowFillWeights){
2063  fh3BeforeNUAWeightsPion->Fill(track->Phi(), track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
2064  fhBeforeNUEWeightsPion->Fill(track->Pt());
2065  }
2066  if(fFlowUseNUAWeights)
2067  {
2068  if(track->Charge() > 0){ NUAweight = fh3NUAWeightPionPlus->GetBinContent( fh3NUAWeightPionPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2069  if(track->Charge() < 0){ NUAweight = fh3NUAWeightPionMinus->GetBinContent( fh3NUAWeightPionMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2070  fh3AfterNUAWeightsPion->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
2071  }
2072  if(fFlowUseNUEWeights)
2073  {
2074  /*if(fIndexCentrality>=0 && fIndexCentrality<10){
2075  Int_t c = fIndexCentrality/5;
2076  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
2077  Int_t c = 1 + fIndexCentrality/10;
2078  }else if(fIndexCentrality>=50){continue;}
2079  */
2080  if(track->Charge() > 0 && fIndexCentrality<50){ NUEweight = fhNUEWeightPionPlus[fIndexCentrality/10]->GetBinContent( fhNUEWeightPionPlus[fIndexCentrality/10]->FindBin(track->Pt()) );}
2081  if(track->Charge() < 0 && fIndexCentrality<50){ NUEweight = fhNUEWeightPionMinus[fIndexCentrality/10]->GetBinContent( fhNUEWeightPionMinus[fIndexCentrality/10]->FindBin(track->Pt()) );}
2082  fhAfterNUEWeightsPion->Fill(track->Pt(),NUEweight);
2083  }
2084  break;
2085  case kKaon:
2086  fVectorKaon->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kKaon, fPDGMassKaon, track->Px(), track->Py(), track->Pz()) );
2087  if(fRunMode == kFillWeights || fFlowFillWeights){
2088  fh3BeforeNUAWeightsKaon->Fill(track->Phi(), track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
2089  fhBeforeNUEWeightsKaon->Fill(track->Pt());
2090  }
2091  if(fFlowUseNUAWeights)
2092  {
2093  if(track->Charge() > 0){NUAweight = fh3NUAWeightKaonPlus->GetBinContent( fh3NUAWeightKaonPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2094  if(track->Charge() < 0){NUAweight = fh3NUAWeightKaonMinus->GetBinContent( fh3NUAWeightKaonMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2095  fh3AfterNUAWeightsKaon->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
2096  }
2097  if(fFlowUseNUEWeights)
2098  {
2099  /*
2100  if(fIndexCentrality>=0 && fIndexCentrality<10){
2101  Int_t c = fIndexCentrality/5;
2102  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
2103  Int_t c = 1 + fIndexCentrality/10;
2104  }else if(fIndexCentrality>=50){continue;}
2105  */
2106  if(track->Charge() > 0 && fIndexCentrality<50){ NUEweight = fhNUEWeightKaonPlus[fIndexCentrality/10]->GetBinContent( fhNUEWeightKaonPlus[fIndexCentrality/10]->FindBin(track->Pt()) );}
2107  if(track->Charge() < 0 && fIndexCentrality<50){ NUEweight = fhNUEWeightKaonMinus[fIndexCentrality/10]->GetBinContent( fhNUEWeightKaonMinus[fIndexCentrality/10]->FindBin(track->Pt()) );}
2108  fhAfterNUEWeightsKaon->Fill(track->Pt(),NUEweight);
2109  }
2110  break;
2111  case kProton:
2112  fVectorProton->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kProton, fPDGMassProton, track->Px(), track->Py(), track->Pz()) );
2113  if(fRunMode == kFillWeights || fFlowFillWeights){
2114  fh3BeforeNUAWeightsProton->Fill(track->Phi(), track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
2115  fhBeforeNUEWeightsProton->Fill(track->Pt());
2116  }
2117  if(fFlowUseNUAWeights)
2118  {
2119  if(track->Charge() > 0){NUAweight = fh3NUAWeightProtonPlus->GetBinContent( fh3NUAWeightProtonPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2120  if(track->Charge() < 0){NUAweight = fh3NUAWeightProtonMinus->GetBinContent( fh3NUAWeightProtonMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2121  fh3AfterNUAWeightsProton->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
2122  }
2123  if(fFlowUseNUEWeights)
2124  {
2125  /*if(fIndexCentrality>=0 && fIndexCentrality<10){
2126  Int_t c = fIndexCentrality/5;
2127  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
2128  Int_t c = 1 + fIndexCentrality/10;
2129  }else if(fIndexCentrality>=50){continue;}
2130  */
2131  if(track->Charge() > 0 && fIndexCentrality<50){ NUEweight = fhNUEWeightProtonPlus[fIndexCentrality/10]->GetBinContent( fhNUEWeightProtonPlus[fIndexCentrality/10]->FindBin(track->Pt()) );}
2132  if(track->Charge() < 0 && fIndexCentrality<50){ NUEweight = fhNUEWeightProtonMinus[fIndexCentrality/10]->GetBinContent( fhNUEWeightProtonMinus[fIndexCentrality/10]->FindBin(track->Pt()) );}
2133  fhAfterNUEWeightsProton->Fill(track->Pt(),NUEweight);
2134  }
2135  break;
2136  default:
2137  break;
2138  }//switch
2139  //if(fFillQA) FillPIDQA(1,track,species); // filling QA for tracks AFTER selection
2140  }
2141  if(fAnalType == kMC){
2142  if(fFillQA) FillPIDQA(0,particle,kUnknown); // filling QA for tracks AFTER selection
2143 
2144  switch (TMath::Abs(particle->PdgCode())) {
2145  case 211:
2146  species = kPion;
2147  fVectorPion->emplace_back( FlowPart(particle->Pt(),particle->Phi(),particle->Eta(), particle->Charge(), kPion, fPDGMassPion, particle->Px(), particle->Py(), particle->Pz()) );
2148  break;
2149 
2150  case 321:
2151  species = kKaon;
2152  fVectorKaon->emplace_back( FlowPart(particle->Pt(),particle->Phi(),particle->Eta(), particle->Charge(), kKaon, fPDGMassKaon, particle->Px(), particle->Py(), particle->Pz()) );
2153  break;
2154 
2155  case 2212:
2156  species = kProton;
2157  fVectorProton->emplace_back( FlowPart(particle->Pt(),particle->Phi(),particle->Eta(), particle->Charge(), kProton, fPDGMassProton, particle->Px(), particle->Py(), particle->Pz()) );
2158  break;
2159 
2160  default:
2161  break;
2162  }
2163  if(fFillQA) FillPIDQA(1,particle,species); // filling QA for tracks AFTER selection
2164  }
2165  }
2166  if(fAnalType == kAOD){
2167  fh2PIDPionMult->Fill(fIndexCentrality,fVectorPion->size());
2168  fh2PIDKaonMult->Fill(fIndexCentrality,fVectorKaon->size());
2169  fh2PIDProtonMult->Fill(fIndexCentrality,fVectorProton->size());
2170  }
2171  if(fAnalType == kMC){
2172  fh2PIDPionMult->Fill(fIndexImpactParameter,fVectorPion->size());
2173  fh2PIDKaonMult->Fill(fIndexImpactParameter,fVectorKaon->size());
2174  fh2PIDProtonMult->Fill(fIndexImpactParameter,fVectorProton->size());
2175  }
2176 
2177 
2178  return;
2179 }
2180 //_____________________________________________________________________________
2181 AliAnalysisTaskFlowModesMC::PartSpecies AliAnalysisTaskFlowModesMC::IsPIDSelected(const AliAODTrack* track)
2182 {
2183  // Selection of PID tracks (pi,K,p) - track identification
2184  // nSigma cutting is used
2185  // returns AliAnalysisTaskFlowModesMC::PartSpecies enum : kPion, kKaon, kProton if any of this passed kUnknown otherwise
2186  // *************************************************************
2187 
2188  // checking detector states
2189  AliPIDResponse::EDetPidStatus pidStatusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, track);
2190  AliPIDResponse::EDetPidStatus pidStatusTOF = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
2191 
2192  Bool_t bIsTPCok = (pidStatusTPC == AliPIDResponse::kDetPidOk);
2193  Bool_t bIsTOFok = ((pidStatusTOF == AliPIDResponse::kDetPidOk) && (track->GetStatus()& AliVTrack::kTOFout) && (track->GetStatus()& AliVTrack::kTIME) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000)); // checking TOF
2194  //if(!bIsTPCok) return kUnknown;
2195 
2196  const Double_t dPt = track->Pt();
2197  const Double_t dP = track->P();
2198 
2199  // use nSigma cuts (based on combination of TPC / TOF nSigma cuts)
2200 
2201  Double_t dNumSigmaTPC[5] = {-99,-99,-99,-99,-99}; // TPC nSigma array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2202  Double_t dNumSigmaTOF[5] = {-99,-99,-99,-99,-99}; // TOF nSigma array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2203 
2204  Float_t *probabilities;
2205  Float_t mismProb;
2206  Float_t ProbBayes[5] = {0,0,0,0,0}; //0=el, 1=mu, 2=pi, 3=ka, 4=pr, 5=deuteron, 6=triton, 7=He3
2207 
2208  Double_t probTPC[AliPID::kSPECIES]={0.};
2209  Double_t probTOF[AliPID::kSPECIES]={0.};
2210  Double_t probTPCTOF[AliPID::kSPECIES]={0.};
2211 
2212 
2213  // filling nSigma arrays
2214  if(bIsTPCok) // should be anyway
2215  {
2216  dNumSigmaTPC[0] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron));
2217  dNumSigmaTPC[1] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kMuon));
2218  dNumSigmaTPC[2] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion));
2219  dNumSigmaTPC[3] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
2220  dNumSigmaTPC[4] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton));
2221  }
2222 
2223  if(bIsTOFok) // should be anyway
2224  {
2225  dNumSigmaTOF[0] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron));
2226  dNumSigmaTOF[1] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kMuon));
2227  dNumSigmaTOF[2] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion));
2228  dNumSigmaTOF[3] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon));
2229  dNumSigmaTOF[4] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton));
2230  }
2231 
2232 
2233  if(fPIDbayesian){
2234  if(!TPCTOFagree(track)){return kUnknown;}
2235 
2236  fBayesianResponse->SetDetResponse(fEventAOD, fCurrCentr,AliESDpid::kTOF_T0); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
2237  if(fEventAOD->GetTOFHeader()){
2238  fESDpid.SetTOFResponse(fEventAOD,AliESDpid::kTOF_T0);
2239  }
2240  fBayesianResponse->SetDetAND(1);
2241  }
2242 
2243  Int_t ParticleFlag[]={0,0,0,0};//Unknown (electron), pion, kaon, proton
2244  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
2245  UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
2246  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) { // TPC is available
2247 
2248  // TPC nSigma cuts
2249  if(dP>0.2 && dP <= 0.5)
2250  {
2251  if(fPIDnsigma){
2252  Double_t dMinSigmasTPC = TMath::MinElement(5,dNumSigmaTPC);
2253  // electron rejection
2254  if(dMinSigmasTPC == dNumSigmaTPC[0] && TMath::Abs(dNumSigmaTPC[0]) <= fCutPIDnSigmaTPCRejectElectron) ParticleFlag[0] = 1; //return kUnknown;
2255  if(dMinSigmasTPC == dNumSigmaTPC[2] && TMath::Abs(dNumSigmaTPC[2]) <= fCutPIDnSigmaPionMax) ParticleFlag[1] = 1;//return kPion;
2256  if(dMinSigmasTPC == dNumSigmaTPC[3] && TMath::Abs(dNumSigmaTPC[3]) <= fCutPIDnSigmaKaonMax) ParticleFlag[2] = 1;//return kKaon;
2257  if(dMinSigmasTPC == dNumSigmaTPC[4] && TMath::Abs(dNumSigmaTPC[4]) <= fCutPIDnSigmaProtonMax) ParticleFlag[3] = 1;//return kProton;
2258  }
2259  if(fPIDbayesian){
2260  ProbBayes[0] = probTPC[0];
2261  ProbBayes[1] = probTPC[1];
2262  ProbBayes[2] = probTPC[2];
2263  ProbBayes[3] = probTPC[3];
2264  ProbBayes[4] = probTPC[4];
2265 
2266  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2267  if(dMaxBayesianProb > fParticleProbability){
2268  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaTPC[0]) <= fCutPIDnSigmaTPCRejectElectron)ParticleFlag[0] = 1;// return kUnknown;
2269  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaTPC[2]) <= fCutPIDnSigmaPionMax)ParticleFlag[1] = 1;//return kPion;
2270  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaTPC[3]) <= fCutPIDnSigmaKaonMax)ParticleFlag[2] = 1;//return kKaon;
2271  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaTPC[4]) <= fCutPIDnSigmaProtonMax)ParticleFlag[3] = 1;//return kProton;
2272  }
2273  /*
2274  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
2275  probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
2276  ProbBayes[0] = probabilities[0];
2277  ProbBayes[1] = probabilities[1];
2278  ProbBayes[2] = probabilities[2];
2279  ProbBayes[3] = probabilities[3];
2280  ProbBayes[4] = probabilities[4];
2281 
2282  mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
2283 
2284  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2285  if(dMaxBayesianProb > fParticleProbability && mismProb < 0.5){
2286  // electron rejection
2287  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaTPC[0]) <= fCutPIDnSigmaTPCRejectElectron)ParticleFlag[0] = 1;// return kUnknown;
2288  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaTPC[2]) <= fCutPIDnSigmaPionMax)ParticleFlag[1] = 1;//return kPion;
2289  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaTPC[3]) <= fCutPIDnSigmaKaonMax)ParticleFlag[2] = 1;//return kKaon;
2290  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaTPC[4]) <= fCutPIDnSigmaProtonMax)ParticleFlag[3] = 1;//return kProton;
2291  }*/
2292  }
2293  }
2294  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
2295  detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
2296 
2297  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()){
2298  AliAODPid* pidObj = track->GetDetPid();
2299  if ((detUsed >= AliPIDResponse::kDetTOF) && (pidObj && pidObj->GetTOFsignal() < 99999)){
2300 
2301  // combined TPC + TOF nSigma cuts
2302  if(dP > 0.5) // && < 4 GeV TODO once TPC dEdx parametrisation is available
2303  {
2304  Double_t dNumSigmaCombined[5] = {-99,-99,-99,-99,-99};
2305 
2306  // discard candidates if no TOF is available if cut is on
2307  if(fCutPIDnSigmaCombinedNoTOFrejection && !bIsTOFok) ParticleFlag[0] = 1;// return kUnknown;
2308 
2309  // calculating combined nSigmas
2310  for(Short_t i(0); i < 5; i++)
2311  {
2312  if(bIsTOFok) { dNumSigmaCombined[i] = TMath::Sqrt(dNumSigmaTPC[i]*dNumSigmaTPC[i] + dNumSigmaTOF[i]*dNumSigmaTOF[i]); }
2313  else { dNumSigmaCombined[i] = dNumSigmaTPC[i]; }
2314  }
2315 
2316  if(fPIDnsigma){
2317  Double_t dMinSigmasCombined = TMath::MinElement(5,dNumSigmaCombined);
2318 
2319  // electron rejection
2320  if(dMinSigmasCombined == dNumSigmaCombined[0] && TMath::Abs(dNumSigmaCombined[0]) <= fCutPIDnSigmaPionMax) ParticleFlag[0] = 1;// return kUnknown;
2321 
2322  switch (fPIDnsigmaCombination) {
2323  case 1:
2324  //combination 1
2325  if(dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 3.) ParticleFlag[1] = 1;//return kPion;
2326  if(dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.5) ParticleFlag[2] = 1; //return kKaon;
2327  if(dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 3.) ParticleFlag[3] = 1;//return kProton;
2328  break;
2329  case 2:
2330  //combination 2
2331  if(dP < 2.5 && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 3.) ParticleFlag[1] = 1; //return kPion;
2332  if(dP > 2.5 && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 2.) ParticleFlag[1] = 1; //return kPion;
2333 
2334  if(dP < 2. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.5) ParticleFlag[2] = 1; //return kKaon;
2335  if(dP > 2. && dP < 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.) ParticleFlag[2] = 1; //return kKaon;
2336  if(dP > 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 1.5) ParticleFlag[2] = 1; //return kKaon;
2337 
2338  if(dP < 3. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 3.) ParticleFlag[3] = 1; //return kProton;
2339  if(dP > 3. && dP < 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 2.) ParticleFlag[3] = 1; //return kProton;
2340  if(dP > 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 1.5) ParticleFlag[3] = 1; //return kProton;
2341 
2342  break;
2343  case 3:
2344  //combination 3
2345  if(dP < 2.5 && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 3.) ParticleFlag[1] = 1; //return kPion;
2346  if(dP > 2.5 && dP < 4. && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 1.5) ParticleFlag[1] = 1; //return kPion;
2347  if(dP > 4. && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 1.) ParticleFlag[1] = 1; //return kPion;
2348 
2349  if(dP < 2. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.5) ParticleFlag[2] = 1; //return kKaon;
2350  if(dP > 2. && dP < 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 1.5)ParticleFlag[2] = 1; //return kKaon;
2351  if(dP > 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 1.) ParticleFlag[2] = 1; //return kKaon;
2352 
2353  if(dP < 3. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 3.) ParticleFlag[3] = 1; //return kProton;
2354  if(dP > 3. && dP < 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 2.) ParticleFlag[3] = 1; //return kProton;
2355  if(dP > 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 1.) ParticleFlag[3] = 1; //return kProton;
2356  break;
2357  default:
2358  break;
2359  }
2360  }
2361  if(fPIDbayesian){
2362  ProbBayes[0] = probTPCTOF[0];
2363  ProbBayes[1] = probTPCTOF[1];
2364  ProbBayes[2] = probTPCTOF[2];
2365  ProbBayes[3] = probTPCTOF[3];
2366  ProbBayes[4] = probTPCTOF[4];
2367  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2368  if(dMaxBayesianProb > fParticleProbability){
2369  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaCombined[0]) <= fCutPIDnSigmaPionMax) ParticleFlag[0] = 1; //return kUnknown;
2370  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaCombined[2]) <= fCutPIDnSigmaPionMax) ParticleFlag[1] = 1; //return kPion;
2371  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaCombined[3]) <= fCutPIDnSigmaKaonMax) ParticleFlag[2] = 1; //return kKaon;
2372  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaCombined[4]) <= fCutPIDnSigmaProtonMax) ParticleFlag[3] = 1; //return kProton;
2373  }else{ParticleFlag[0] = 1;}
2374  /*
2375  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
2376  probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
2377  ProbBayes[0] = probabilities[0];
2378  ProbBayes[1] = probabilities[1];
2379  ProbBayes[2] = probabilities[2];
2380  ProbBayes[3] = probabilities[3];
2381  ProbBayes[4] = probabilities[4];
2382 
2383  mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
2384 
2385  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2386  if(dMaxBayesianProb > fParticleProbability && mismProb < 0.5){
2387  // electron rejection
2388  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaCombined[0]) <= fCutPIDnSigmaPionMax) ParticleFlag[0] = 1; //return kUnknown;
2389  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaCombined[2]) <= fCutPIDnSigmaPionMax) ParticleFlag[1] = 1; //return kPion;
2390  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaCombined[3]) <= fCutPIDnSigmaKaonMax) ParticleFlag[2] = 1; //return kKaon;
2391  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaCombined[4]) <= fCutPIDnSigmaProtonMax) ParticleFlag[3] = 1; //return kProton;
2392  }*/
2393  }
2394  }//if(dP > 0.5)
2395  }//if ((detUsed >= AliPIDResponse::kDetTOF) && (pidObj && pidObj->GetTOFsignal() < 99999))
2396  }//TPC or TOF if (detUsed)
2397  }//TPC if (detUsed)
2398 
2399  PartSpecies species;
2400  if(!ParticleFlag[0] && ParticleFlag[1] && !ParticleFlag[2] && !ParticleFlag[3]){
2401  species = kPion;
2402  }else if(!ParticleFlag[0] && !ParticleFlag[1] && ParticleFlag[2] && !ParticleFlag[3]){
2403  species = kKaon;
2404  }else if(!ParticleFlag[0] && !ParticleFlag[1] && !ParticleFlag[2] && ParticleFlag[3]){
2405  species = kProton;
2406  }else{species = kUnknown;}
2407 
2408  if(fFillQA) FillPIDQA(1,track,species); // filling QA for tracks AFTER selection
2409 
2410  return species;
2411 }
2412 //_____________________________________________________________________________
2413 void AliAnalysisTaskFlowModesMC::FillPIDQA(const Short_t iQAindex, const AliAODTrack* track, const PartSpecies species)
2414 {
2415  // Filling various QA plots related to PID (pi,K,p) track selection
2416  // *************************************************************
2417  if(!track) return;
2418  if(!fPIDResponse || !fPIDCombined)
2419  {
2420  ::Error("FillPIDQA","AliPIDResponse or AliPIDCombined object not found!");
2421  return;
2422  }
2423 
2424  // TPC & TOF statuses & measures
2425  AliPIDResponse::EDetPidStatus pidStatusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, track);
2426  AliPIDResponse::EDetPidStatus pidStatusTOF = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
2427 
2428  fhQAPIDTOFstatus[iQAindex]->Fill((Int_t) pidStatusTOF );
2429  fhQAPIDTPCstatus[iQAindex]->Fill((Int_t) pidStatusTPC );
2430 
2431  Bool_t bIsTPCok = (pidStatusTPC == AliPIDResponse::kDetPidOk);
2432  //Bool_t bIsTOFok = ((pidStatusTOF == AliPIDResponse::kDetPidOk) && (track->GetStatus()& AliVTrack::kTOFout) && (track->GetStatus()& AliVTrack::kTIME));
2433 
2434  Bool_t bIsTOFok = ((pidStatusTOF == AliPIDResponse::kDetPidOk) && (track->GetStatus()& AliVTrack::kTOFout) && (track->GetStatus()& AliVTrack::kTIME) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000)); // checking TOF
2435 
2436  Double_t dNumSigmaTPC[5] = {-11}; // array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2437  Double_t dNumSigmaTOF[5] = {-11}; // array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2438 
2439  Double_t dTPCdEdx = -5; // TPC dEdx for selected particle
2440  Double_t dTOFbeta = -0.05; //TOF beta for selected particle
2441 
2442  Double_t dP = track->P();
2443  Double_t dPt = track->Pt();
2444 
2445  // detector status dependent
2446  if(bIsTPCok)
2447  {
2448  dNumSigmaTPC[0] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2449  dNumSigmaTPC[1] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kMuon);
2450  dNumSigmaTPC[2] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
2451  dNumSigmaTPC[3] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
2452  dNumSigmaTPC[4] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2453 
2454  dTPCdEdx = track->GetTPCsignal();
2455  fhQAPIDTPCdEdx[iQAindex]->Fill(track->P(), dTPCdEdx);
2456  }
2457  else // TPC status not OK
2458  {
2459  dNumSigmaTPC[0] = -11.;
2460  dNumSigmaTPC[1] = -11.;
2461  dNumSigmaTPC[2] = -11.;
2462  dNumSigmaTPC[3] = -11.;
2463  dNumSigmaTPC[4] = -11.;
2464 
2465  fhQAPIDTPCdEdx[iQAindex]->Fill(track->P(), -5.);
2466  }
2467 
2468  if(bIsTOFok)
2469  {
2470  dNumSigmaTOF[0] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
2471  dNumSigmaTOF[1] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kMuon);
2472  dNumSigmaTOF[2] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
2473  dNumSigmaTOF[3] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
2474  dNumSigmaTOF[4] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
2475 
2476  Double_t dTOF[5];
2477  track->GetIntegratedTimes(dTOF);
2478  dTOFbeta = dTOF[0] / track->GetTOFsignal();
2479  fhQAPIDTOFbeta[iQAindex]->Fill(dP,dTOFbeta);
2480  }
2481  else // TOF status not OK
2482  {
2483  dNumSigmaTOF[0] = -11.;
2484  dNumSigmaTOF[1] = -11.;
2485  dNumSigmaTOF[2] = -11.;
2486  dNumSigmaTOF[3] = -11.;
2487  dNumSigmaTOF[4] = -11.;
2488 
2489  fhQAPIDTOFbeta[iQAindex]->Fill(track->P(),-0.05);
2490  }
2491 
2492 
2493  // species dependent QA
2494  switch (species)
2495  {
2496  case kPion:
2497  fh2PIDPionPt->Fill(fIndexCentrality,track->Pt());
2498  fh2PIDPionPhi->Fill(fIndexCentrality,track->Phi());
2499  fh2PIDPionEta->Fill(fIndexCentrality,track->Eta());
2500  fhPIDPionCharge->Fill(track->Charge());
2501  fh2PIDPionTPCdEdx->Fill(dPt,dTPCdEdx);
2502  fh2PIDPionTOFbeta->Fill(dPt,dTOFbeta);
2503  fh2PIDPionTPCnSigmaPion->Fill(dPt,dNumSigmaTPC[2]);
2504  fh2PIDPionTOFnSigmaPion->Fill(dPt,dNumSigmaTOF[2]);
2505  fh2PIDPionTPCnSigmaKaon->Fill(dPt,dNumSigmaTPC[3]);
2506  fh2PIDPionTOFnSigmaKaon->Fill(dPt,dNumSigmaTOF[3]);
2507  fh2PIDPionTPCnSigmaProton->Fill(dPt,dNumSigmaTPC[4]);
2508  fh2PIDPionTOFnSigmaProton->Fill(dPt,dNumSigmaTOF[4]);
2509 
2510  fh3PIDPionTPCTOFnSigmaPion[iQAindex]->Fill(dNumSigmaTPC[2],dNumSigmaTOF[2],dPt);
2511  fh3PIDPionTPCTOFnSigmaKaon[iQAindex]->Fill(dNumSigmaTPC[3],dNumSigmaTOF[3],dPt);
2512  fh3PIDPionTPCTOFnSigmaProton[iQAindex]->Fill(dNumSigmaTPC[4],dNumSigmaTOF[4],dPt);
2513 
2514  break;
2515 
2516  case kKaon:
2517  fh2PIDKaonPt->Fill(fIndexCentrality,track->Pt());
2518  fh2PIDKaonPhi->Fill(fIndexCentrality,track->Phi());
2519  fh2PIDKaonEta->Fill(fIndexCentrality,track->Eta());
2520  fhPIDKaonCharge->Fill(track->Charge());
2521  fh2PIDKaonTPCdEdx->Fill(dP,dTPCdEdx);
2522  fh2PIDKaonTOFbeta->Fill(dP,dTOFbeta);
2523  fh2PIDKaonTPCnSigmaPion->Fill(dPt,dNumSigmaTPC[2]);
2524  fh2PIDKaonTOFnSigmaPion->Fill(dPt,dNumSigmaTOF[2]);
2525  fh2PIDKaonTPCnSigmaKaon->Fill(dPt,dNumSigmaTPC[3]);
2526  fh2PIDKaonTOFnSigmaKaon->Fill(dPt,dNumSigmaTOF[3]);
2527  fh2PIDKaonTPCnSigmaProton->Fill(dPt,dNumSigmaTPC[4]);
2528  fh2PIDKaonTOFnSigmaProton->Fill(dPt,dNumSigmaTOF[4]);
2529 
2530  fh3PIDKaonTPCTOFnSigmaPion[iQAindex]->Fill(dNumSigmaTPC[2],dNumSigmaTOF[2],dPt);
2531  fh3PIDKaonTPCTOFnSigmaKaon[iQAindex]->Fill(dNumSigmaTPC[3],dNumSigmaTOF[3],dPt);
2532  fh3PIDKaonTPCTOFnSigmaProton[iQAindex]->Fill(dNumSigmaTPC[4],dNumSigmaTOF[4],dPt);
2533 
2534  break;
2535 
2536  case kProton:
2537  fh2PIDProtonPt->Fill(fIndexCentrality,track->Pt());
2538  fh2PIDProtonPhi->Fill(fIndexCentrality,track->Phi());
2539  fh2PIDProtonEta->Fill(fIndexCentrality,track->Eta());
2540  fhPIDProtonCharge->Fill(track->Charge());
2541  fh2PIDProtonTPCdEdx->Fill(dP,dTPCdEdx);
2542  fh2PIDProtonTOFbeta->Fill(dP,dTOFbeta);
2543  fh2PIDProtonTPCnSigmaPion->Fill(dPt,dNumSigmaTPC[2]);
2544  fh2PIDProtonTOFnSigmaPion->Fill(dPt,dNumSigmaTOF[2]);
2545  fh2PIDProtonTPCnSigmaKaon->Fill(dPt,dNumSigmaTPC[3]);
2546  fh2PIDProtonTOFnSigmaKaon->Fill(dPt,dNumSigmaTOF[3]);
2547  fh2PIDProtonTPCnSigmaProton->Fill(dPt,dNumSigmaTPC[4]);
2548  fh2PIDProtonTOFnSigmaProton->Fill(dPt,dNumSigmaTOF[4]);
2549 
2550  fh3PIDProtonTPCTOFnSigmaPion[iQAindex]->Fill(dNumSigmaTPC[2],dNumSigmaTOF[2],dPt);
2551  fh3PIDProtonTPCTOFnSigmaKaon[iQAindex]->Fill(dNumSigmaTPC[3],dNumSigmaTOF[3],dPt);
2552  fh3PIDProtonTPCTOFnSigmaProton[iQAindex]->Fill(dNumSigmaTPC[4],dNumSigmaTOF[4],dPt);
2553 
2554  break;
2555 
2556  default:
2557  break;
2558  }
2559 
2560 
2561  //
2562 
2563 
2564  return;
2565 }
2566 
2567 //_____________________________________________________________________________
2568 void AliAnalysisTaskFlowModesMC::FillPIDQA(const Short_t iQAindex, const AliMCParticle* particle, const PartSpecies species)
2569 {
2570  switch (species)
2571  {
2572  case kPion:
2573  fh2PIDPionPt->Fill(fIndexImpactParameter,particle->Pt());
2574  fh2PIDPionPhi->Fill(fIndexImpactParameter,particle->Phi());
2575  fh2PIDPionEta->Fill(fIndexImpactParameter,particle->Eta());
2576  fhPIDPionCharge->Fill(particle->Charge());
2577  break;
2578  case kKaon:
2579  fh2PIDKaonPt->Fill(fIndexImpactParameter,particle->Pt());
2580  fh2PIDKaonPhi->Fill(fIndexImpactParameter,particle->Phi());
2581  fh2PIDKaonEta->Fill(fIndexImpactParameter,particle->Eta());
2582  fhPIDKaonCharge->Fill(particle->Charge());
2583  break;
2584  case kProton:
2585  fh2PIDProtonPt->Fill(fIndexImpactParameter,particle->Pt());
2586  fh2PIDProtonPhi->Fill(fIndexImpactParameter,particle->Phi());
2587  fh2PIDProtonEta->Fill(fIndexImpactParameter,particle->Eta());
2588  fhPIDProtonCharge->Fill(particle->Charge());
2589  break;
2590  default:
2591  break;
2592  }
2593 
2594 }
2595 //_____________________________________________________________________________
2596 Bool_t AliAnalysisTaskFlowModesMC::ProcessEvent()
2597 {
2598  // main method for processing of (selected) event:
2599  // - Filtering of tracks / particles for flow calculations
2600  // - Phi,eta,pt weights for generic framework are calculated if specified
2601  // - Flow calculations
2602  // returns kTRUE if succesfull
2603  // *************************************************************
2604 
2605  // printf("======= EVENT ================\n");
2606  if(fAnalType==kAOD){
2607  // checking the run number for applying weights & loading TList with weights
2608  if(fRunNumber < 0 || fRunNumber != fEventAOD->GetRunNumber() )
2609  {
2610  fRunNumber = fEventAOD->GetRunNumber();
2611 
2612  if(fFlowUseNUAWeights && fFlowNUAWeightsFile)
2613  {
2614  TDirectory* dirFlowNUAWeights = (TDirectory*) fFlowNUAWeightsFile->Get(Form("000%d",fRunNumber));
2615  if(!dirFlowNUAWeights) {::Error("ProcessEvent","TList from flow weights not found."); return kFALSE; }
2616  fh3NUAWeightRefsPlus = (TH3D*) dirFlowNUAWeights->Get("ChargedPlus"); if(!fh3NUAWeightRefsPlus) { ::Error("ProcessEvent","Positive Refs NUA weights not found"); return kFALSE; }
2617  fh3NUAWeightRefsMinus = (TH3D*) dirFlowNUAWeights->Get("ChargedMinus"); if(!fh3NUAWeightRefsMinus) { ::Error("ProcessEvent","Negative Refs NUA weights not found"); return kFALSE; }
2618 
2619  fh3NUAWeightChargedPlus = (TH3D*) dirFlowNUAWeights->Get("ChargedPlus"); if(!fh3NUAWeightChargedPlus) { ::Error("ProcessEvent","Positive Charged NUA weights not found"); return kFALSE; }
2620  fh3NUAWeightChargedMinus = (TH3D*) dirFlowNUAWeights->Get("ChargedMinus"); if(!fh3NUAWeightChargedMinus) { ::Error("ProcessEvent","Nagative Charged NUA weights not found"); return kFALSE; }
2621 
2622  fh3NUAWeightPionPlus = (TH3D*) dirFlowNUAWeights->Get("PionPlus"); if(!fh3NUAWeightPionPlus) { ::Error("ProcessEvent","Positive Pion NUA weights not found"); return kFALSE; }
2623  fh3NUAWeightPionMinus = (TH3D*) dirFlowNUAWeights->Get("PionMinus"); if(!fh3NUAWeightPionMinus) { ::Error("ProcessEvent","Negative Pion NUA weights not found"); return kFALSE; }
2624 
2625 
2626  fh3NUAWeightKaonPlus = (TH3D*) dirFlowNUAWeights->Get("KaonPlus"); if(!fh3NUAWeightKaonPlus) { ::Error("ProcessEvent","Positive Kaon NUA weights not found"); return kFALSE; }
2627  fh3NUAWeightKaonMinus = (TH3D*) dirFlowNUAWeights->Get("KaonMinus"); if(!fh3NUAWeightKaonMinus) { ::Error("ProcessEvent","Negative Kaon NUA weights not found"); return kFALSE; }
2628 
2629 
2630  fh3NUAWeightProtonPlus = (TH3D*) dirFlowNUAWeights->Get("ProtonPlus"); if(!fh3NUAWeightProtonPlus) { ::Error("ProcessEvent","Positive Proton NUA weights not found"); return kFALSE; }
2631  fh3NUAWeightProtonMinus = (TH3D*) dirFlowNUAWeights->Get("ProtonMinus"); if(!fh3NUAWeightProtonMinus) { ::Error("ProcessEvent","Negative Proton NUA weights not found"); return kFALSE; }
2632 
2633 
2634  }
2635  }
2636 
2637  if(fColSystem == kPbPb){
2638  fIndexCentrality = GetCentralityIndex();
2639  if(fIndexCentrality < 0) { return kFALSE;}
2640  }
2641 
2642  if(fFlowUseNUEWeights && fFlowNUEWeightsFile)
2643  {
2644  TDirectory* dirFlowNUEWeights[5] = {0x0};
2645 
2646  const char* gCentrality[] = {"0-10cc","10-20cc","20-30cc","30-40cc","40-50cc"};
2647  for(int iCentrality = 0; iCentrality<fNumCentralities;iCentrality++){
2648  dirFlowNUEWeights[iCentrality] = (TDirectory*) fFlowNUEWeightsFile->Get(gCentrality[iCentrality]);
2649  if(!dirFlowNUEWeights[iCentrality]) {::Error("ProcessEvent","TDirectoy from NUE weights not found."); return kFALSE; }
2650  fhNUEWeightRefsPlus[iCentrality] = dynamic_cast<TH1F*>( dirFlowNUEWeights[iCentrality]->FindObject("ChargedPlus")); if(!fhNUEWeightRefsPlus[iCentrality]) { ::Error("ProcessEvent","Positive Refs NUE weights not found"); return kFALSE; }
2651  fhNUEWeightRefsMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ChargedMinus")); if(!fhNUEWeightRefsMinus[iCentrality]) { ::Error("ProcessEvent","Negative Refs NUE weights not found"); return kFALSE; }
2652  fhNUEWeightChargedPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ChargedPlus")); if(!fhNUEWeightChargedPlus[iCentrality]) { ::Error("ProcessEvent","Positive Charged NUE weights not found"); return kFALSE; }
2653  fhNUEWeightChargedMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ChargedMinus")); if(!fhNUEWeightChargedMinus[iCentrality]) { ::Error("ProcessEvent","Negative Charged NUE weights not found"); return kFALSE; }
2654  fhNUEWeightPionPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("PionPlus")); if(!fhNUEWeightPionPlus[iCentrality]) { ::Error("ProcessEvent","Positive Pion NUE weights not found"); return kFALSE; }
2655  fhNUEWeightPionMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("PionMinus")); if(!fhNUEWeightPionMinus[iCentrality]) { ::Error("ProcessEvent","Negative Pion NUE weights not found"); return kFALSE; }
2656  fhNUEWeightKaonPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("KaonPlus")); if(!fhNUEWeightKaonPlus[iCentrality]) { ::Error("ProcessEvent","Positive Kaon NUE weights not found"); return kFALSE; }
2657  fhNUEWeightKaonMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("KaonMinus")); if(!fhNUEWeightKaonMinus[iCentrality]) { ::Error("ProcessEvent","Negative Kaon NUE weights not found"); return kFALSE; }
2658  fhNUEWeightProtonPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ProtonPlus")); if(!fhNUEWeightProtonPlus[iCentrality]) { ::Error("ProcessEvent","Positive Proton NUE weights not found"); return kFALSE; }
2659  fhNUEWeightProtonMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ProtonMinus")); if(!fhNUEWeightProtonMinus[iCentrality]) { ::Error("ProcessEvent","Negative Proton NUE weights not found"); return kFALSE; }
2660 
2661  }
2662  }
2663  }
2664 
2665  // filtering particles
2666  if(!Filtering()) return kFALSE;
2667  // at this point, centrality index (percentile) should be properly estimated, if not, skip event
2668 
2669  // if running in kFillWeights mode, skip the remaining part
2670  if(fRunMode == kFillWeights) { fEventCounter++; return kTRUE; }
2671 
2672  // checking if there is at least one charged track selected;
2673  // if not, event is skipped: unable to compute Reference flow (and thus any differential flow)
2674  if(fVectorCharged->size() < 1) return kFALSE;
2675 
2676  // at this point, all particles fullfiling relevant POIs (and REFs) criteria are filled in TClonesArrays
2677 
2678  // >>>> flow starts here <<<<
2679  // >>>> Flow a la General Framework <<<<
2680  for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
2681  {
2682 
2683  // Reference (pT integrated) flow
2684  DoFlowRefs(iGap);
2685 
2686  // pT differential
2687  if(fProcessCharged)
2688  {
2689  // charged track flow
2690  DoFlowCharged(iGap);
2691 
2692  }
2693 
2694  if(fProcessPID)
2695  {
2696  const Int_t iSizePion = fVectorPion->size();
2697  const Int_t iSizeKaon = fVectorKaon->size();
2698  const Int_t iSizeProton = fVectorProton->size();
2699 
2700  // pi,K,p flow
2701  if(iSizePion > 0) DoFlowPID(iGap,kPion);
2702 
2703  if(iSizeKaon > 0) DoFlowPID(iGap,kKaon);
2704 
2705  if(iSizeProton > 0) DoFlowPID(iGap,kProton);
2706 
2707  }
2708  } // endfor {iGap} eta gaps
2709 
2710  fEventCounter++; // counter of processed events
2711  //printf("event %d\n",fEventCounter);
2712 
2713  return kTRUE;
2714 }
2715 //_____________________________________________________________________________
2716 void AliAnalysisTaskFlowModesMC::DoFlowRefs(const Short_t iEtaGapIndex)
2717 {
2718 
2719  // Estimate <2>, <4> and <6> for reference flow for all harmonics based on relevant flow vectors
2720  // *************************************************************
2721 
2722  //Float_t dEtaGap = fEtaGap[iEtaGapIndex];
2723  Short_t iHarmonics = 0;
2724  Double_t Cn2 = 0;
2725  Double_t Dn4GapP = 0;
2726  Double_t Dn6GapP = 0;
2727  TComplex vector = TComplex(0,0,kFALSE);
2728  Double_t dValue = 999;
2729 
2730  FillRefsVectors(iEtaGapIndex); // filling RFPs (Q) flow vectors
2731  if(!fDoOnlyMixedCorrelations){
2732  // estimating <2>
2733  Cn2 = TwoGap(0,0).Re();
2734  if(Cn2 != 0)
2735  {
2736  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2737  {
2738  iHarmonics = fHarmonics[iHarm];
2739  vector = TwoGap(iHarmonics,-iHarmonics);
2740  dValue = vector.Re()/Cn2;
2741  // printf("Gap (RFPs): %g Harm %d | Dn2: %g | fFlowVecQpos[0][0]: %g | fFlowVecQneg[0][0]: %g | fIndexCentrality %d\n\n", dEtaGap,iHarmonics,Cn2,fFlowVecQpos[0][0].Re(),fFlowVecQneg[0][0].Re(),fIndexCentrality);
2742  if( TMath::Abs(dValue < 1) )
2743  if(fAnalType == kAOD) fpRefsCor2[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexCentrality, dValue, Cn2);
2744  if(fAnalType == kMC) fpRefsCor2[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexImpactParameter, dValue, Cn2);
2745 
2746  }
2747  }
2748  }
2749  //for mixed harmonics
2750  if(fDoOnlyMixedCorrelations){
2751  //estimating <4>
2752  Dn4GapP = FourGapPos(0,0,0,0).Re();
2753  if(Dn4GapP != 0)
2754  {
2755  // (2,2 | 2,2)_gap , referece flow for v4/psi2
2756  TComplex Four_2222_GapP = FourGapPos(2, 2, -2, -2);
2757  double c4_2222_GapP = Four_2222_GapP.Re()/Dn4GapP;
2758  if(fAnalType == kAOD) fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][0]->Fill(fIndexCentrality, c4_2222_GapP, Dn4GapP);
2759  if(fAnalType == kMC) fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][0]->Fill(fIndexImpactParameter, c4_2222_GapP, Dn4GapP);
2760 
2761  // (3,3 | 3,3)_gap, referece flow for v6/psi3
2762  TComplex Four_3333_GapP = FourGapPos(3, 3, -3, -3);
2763  double c4_3333_GapP = Four_3333_GapP.Re()/Dn4GapP;
2764  if(fAnalType == kAOD) fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][1]->Fill(fIndexCentrality, c4_3333_GapP, Dn4GapP);
2765  if(fAnalType == kMC) fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][1]->Fill(fIndexImpactParameter, c4_3333_GapP, Dn4GapP);
2766 
2767  // (3,2 | 3,2)_gap, reference flow for v5/psi23
2768  TComplex Four_3232_GapP = FourGapPos(3, 2, -3, -2);
2769  double c4_3232_GapP = Four_3232_GapP.Re()/Dn4GapP;
2770  if(fAnalType == kAOD) fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][2]->Fill(fIndexCentrality, c4_3232_GapP, Dn4GapP);
2771  if(fAnalType == kMC) fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][2]->Fill(fIndexImpactParameter, c4_3232_GapP, Dn4GapP);
2772  }
2773  //estimating <6>
2774  Dn6GapP = SixGapPos(0,0,0,0,0,0).Re();
2775  if(Dn6GapP != 0)
2776  {
2777  // (2,2,2 | 2,2,2)_gap, reference flow for v6/psi2
2778  TComplex Six_222222_GapP = SixGapPos(2, 2, 2, -2, -2, -2);
2779  double c6_222222_GapP = Six_222222_GapP.Re()/Dn6GapP;
2780  if(fAnalType == kAOD) fpMixedRefsCor6[fIndexSampling][iEtaGapIndex]->Fill(fIndexCentrality, c6_222222_GapP, Dn6GapP);
2781  if(fAnalType == kMC) fpMixedRefsCor6[fIndexSampling][iEtaGapIndex]->Fill(fIndexImpactParameter, c6_222222_GapP, Dn6GapP);
2782  }
2783  }
2784  return;
2785 }
2786 //_____________________________________________________________________________
2787 void AliAnalysisTaskFlowModesMC::DoFlowCharged(const Short_t iEtaGapIndex)
2788 {
2789  // Estimate <2'>, <3'> and <4'> for pT diff flow of charged tracks for all harmonics based on relevant flow vectors
2790  // *************************************************************
2791 
2792  FillPOIsVectors(iEtaGapIndex,kCharged); // filling POIs (P,S) flow vectors
2793 
2794  const Double_t dPtBinWidth = (fFlowPOIsPtMax - fFlowPOIsPtMin) / fFlowPOIsPtNumBins;
2795 
2796  //Float_t dEtaGap = fEtaGap[iEtaGapIndex];
2797  Short_t iHarmonics = 0;
2798  Double_t Dn2 = 0;
2799  Double_t DDn3GapP=0;
2800  Double_t DDn4GapP=0;
2801  TComplex vector = TComplex(0,0,kFALSE);
2802  Double_t dValue = 999;
2803 
2804 
2805  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
2806  {
2807  if(!fDoOnlyMixedCorrelations){
2808  // estimating <2'>
2809  // POIs in positive eta
2810  Dn2 = TwoDiffGapPos(0,0,iPt).Re();
2811  if(Dn2 != 0)
2812  {
2813  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2814  {
2815  iHarmonics = fHarmonics[iHarm];
2816  vector = TwoDiffGapPos(iHarmonics,-iHarmonics,iPt);
2817  dValue = vector.Re()/Dn2;
2818  if( TMath::Abs(dValue < 1) )
2819  if(fAnalType == kAOD) fp2ChargedCor2Pos[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
2820  if(fAnalType == kMC) fp2ChargedCor2Pos[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexImpactParameter, iPt*dPtBinWidth, dValue, Dn2);
2821  }
2822  }
2823 
2824  // POIs in negative eta
2825  Dn2 = TwoDiffGapNeg(0,0,iPt).Re();
2826  if(Dn2 != 0)
2827  {
2828  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2829  {
2830  iHarmonics = fHarmonics[iHarm];
2831  vector = TwoDiffGapNeg(iHarmonics,-iHarmonics,iPt);
2832  dValue = vector.Re()/Dn2;
2833  if( TMath::Abs(dValue < 1) )
2834  if(fAnalType == kAOD) fp2ChargedCor2Neg[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
2835  if(fAnalType == kMC) fp2ChargedCor2Neg[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexImpactParameter, iPt*dPtBinWidth, dValue, Dn2);
2836  }
2837  }
2838  }
2839  if(fDoOnlyMixedCorrelations){
2840  // estimating <3'>
2841  // POIs in positive eta
2842  DDn3GapP = ThreeDiffGapPos(0,0,0,iPt).Re();
2843  if(DDn3GapP!=0)
2844  {
2845  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
2846  {
2847  if(iMixedHarm==0){ vector = ThreeDiffGapPos(4,-2,-2,iPt);}
2848  if(iMixedHarm==1){ vector = ThreeDiffGapPos(6,-3,-3,iPt);}
2849  if(iMixedHarm==2){ vector = ThreeDiffGapPos(5,-3,-2,iPt);}
2850  dValue = vector.Re()/DDn3GapP;
2851  if( TMath::Abs(dValue < 1) ){
2852  if(fAnalType == kAOD) fpMixedChargedCor3Pos[fIndexSampling][iEtaGapIndex][iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
2853  if(fAnalType == kMC) fpMixedChargedCor3Pos[fIndexSampling][iEtaGapIndex][iMixedHarm]->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn3GapP);
2854  }
2855  }
2856  }
2857  // POIs in negative eta
2858  DDn3GapP = ThreeDiffGapNeg(0,0,0,iPt).Re();
2859  if(DDn3GapP!=0)
2860  {
2861  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
2862  {
2863  if(iMixedHarm==0){ vector = ThreeDiffGapNeg(4,-2,-2,iPt);}
2864  if(iMixedHarm==1){ vector = ThreeDiffGapNeg(6,-3,-3,iPt);}
2865  if(iMixedHarm==2){ vector = ThreeDiffGapNeg(5,-2,-3,iPt);}
2866  dValue = vector.Re()/DDn3GapP;
2867  if( TMath::Abs(dValue < 1) ){
2868  if(fAnalType == kAOD) fpMixedChargedCor3Neg[fIndexSampling][iEtaGapIndex][iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
2869  if(fAnalType == kMC) fpMixedChargedCor3Neg[fIndexSampling][iEtaGapIndex][iMixedHarm]->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn3GapP);
2870  }
2871  }
2872  }
2873  // estimating <4'>
2874  // POIs in positive eta
2875  DDn4GapP = Four13DiffGapPos(0,0,0,0,iPt).Re();
2876  if(DDn4GapP!=0)
2877  {
2878  vector = Four13DiffGapPos(6,-2,-2,-2,iPt);
2879  if( TMath::Abs(dValue < 1) ){
2880  if(fAnalType == kAOD) fpMixedChargedCor4Pos[fIndexSampling][iEtaGapIndex]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
2881  if(fAnalType == kMC) fpMixedChargedCor4Pos[fIndexSampling][iEtaGapIndex]->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn4GapP);
2882  }
2883  }
2884  // POIs in negative eta
2885  DDn4GapP = Four13DiffGapNeg(0,0,0,0,iPt).Re();
2886  if(DDn4GapP!=0)
2887  {
2888  vector = Four13DiffGapNeg(6,-2,-2,-2,iPt);
2889  if( TMath::Abs(dValue < 1) ){
2890  if(fAnalType == kAOD) fpMixedChargedCor4Neg[fIndexSampling][iEtaGapIndex]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
2891  if(fAnalType == kMC) fpMixedChargedCor4Neg[fIndexSampling][iEtaGapIndex]->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn4GapP);
2892  }
2893  }
2894  }
2895  } // endfor {iPt}
2896  return;
2897 }
2898 //_____________________________________________________________________________
2899 void AliAnalysisTaskFlowModesMC::DoFlowPID(const Short_t iEtaGapIndex, const PartSpecies species)
2900 {
2901  // Estimate <2'>, <3'> and <4'> for pT diff flow of pi/K/p tracks for all harmonics based on relevant flow vectors
2902  // *************************************************************
2903 
2904  TProfile2D** profile2Pos = 0x0;
2905  TProfile2D** profile2Neg = 0x0;
2906  //TProfile2D** profile4 = 0x0;
2907  TProfile2D** profile3Pos = 0x0;
2908  TProfile2D** profile3Neg = 0x0;
2909  TProfile2D* profile4Pos = 0x0;
2910  TProfile2D* profile4Neg = 0x0;
2911 
2912 
2913  if(!fDoOnlyMixedCorrelations){
2914 
2915  switch (species)
2916  {
2917  case kPion:
2918  profile2Pos = fp2PionCor2Pos[fIndexSampling][iEtaGapIndex];
2919  profile2Neg = fp2PionCor2Neg[fIndexSampling][iEtaGapIndex];
2920  //profile4 = fp2PionCor4;
2921  break;
2922 
2923  case kKaon:
2924  profile2Pos = fp2KaonCor2Pos[fIndexSampling][iEtaGapIndex];
2925  profile2Neg = fp2KaonCor2Neg[fIndexSampling][iEtaGapIndex];
2926  //profile4 = fp2KaonCor4;
2927  break;
2928 
2929  case kProton:
2930  profile2Pos = fp2ProtonCor2Pos[fIndexSampling][iEtaGapIndex];
2931  profile2Neg = fp2ProtonCor2Neg[fIndexSampling][iEtaGapIndex];
2932  //profile4 = fp2ProtonCor4;
2933  break;
2934 
2935  default:
2936  ::Error("DoFlowPID","Unexpected species! Terminating!");
2937  return;
2938  }
2939  }
2940 
2941  if(fDoOnlyMixedCorrelations){
2942 
2943  switch (species)
2944  {
2945  case kPion:
2946  profile3Pos = fpMixedPionCor3Pos[fIndexSampling][iEtaGapIndex];
2947  profile3Neg = fpMixedPionCor3Neg[fIndexSampling][iEtaGapIndex];
2948 
2949  profile4Pos = fpMixedPionCor4Pos[fIndexSampling][iEtaGapIndex];
2950  profile4Neg = fpMixedPionCor4Neg[fIndexSampling][iEtaGapIndex];
2951  break;
2952 
2953  case kKaon:
2954  profile3Pos = fpMixedKaonCor3Pos[fIndexSampling][iEtaGapIndex];
2955  profile3Neg = fpMixedKaonCor3Neg[fIndexSampling][iEtaGapIndex];
2956 
2957  profile4Pos = fpMixedKaonCor4Pos[fIndexSampling][iEtaGapIndex];
2958  profile4Neg = fpMixedKaonCor4Neg[fIndexSampling][iEtaGapIndex];
2959  break;
2960 
2961  case kProton:
2962  profile3Pos = fpMixedProtonCor3Pos[fIndexSampling][iEtaGapIndex];
2963  profile3Neg = fpMixedProtonCor3Neg[fIndexSampling][iEtaGapIndex];
2964 
2965  profile4Pos = fpMixedProtonCor4Pos[fIndexSampling][iEtaGapIndex];
2966  profile4Neg = fpMixedProtonCor4Neg[fIndexSampling][iEtaGapIndex];
2967  break;
2968 
2969  default:
2970  ::Error("DoFlowPID with mixed harmonics","Unexpected species! Terminating!");
2971  return;
2972  }
2973  }
2974 
2975  FillPOIsVectors(iEtaGapIndex,species); // Filling POIs vectors
2976 
2977  const Double_t dPtBinWidth = (fFlowPOIsPtMax - fFlowPOIsPtMin) / fFlowPOIsPtNumBins;
2978 
2979  //Float_t dEtaGap = fEtaGap[iEtaGapIndex];
2980  Short_t iHarmonics = 0;
2981  Double_t Dn2 = 0;
2982  TComplex vector = TComplex(0,0,kFALSE);
2983  Double_t dValue = 999;
2984  Double_t DDn3GapP=0;
2985  Double_t DDn4GapP=0;
2986 
2987 
2988  // filling POIs (P,S) flow vectors
2989 
2990  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
2991  {
2992  if(!fDoOnlyMixedCorrelations){
2993  // estimating <2'>
2994  // POIs in positive eta
2995  Dn2 = TwoDiffGapPos(0,0,iPt).Re();
2996  if(Dn2 != 0)
2997  {
2998  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2999  {
3000  iHarmonics = fHarmonics[iHarm];
3001  vector = TwoDiffGapPos(iHarmonics,-iHarmonics,iPt);
3002  dValue = vector.Re()/Dn2;
3003  if( TMath::Abs(dValue < 1) ){
3004  if(fAnalType == kAOD) profile2Pos[iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
3005  if(fAnalType == kMC) profile2Pos[iHarm]->Fill(fIndexImpactParameter, iPt*dPtBinWidth, dValue, Dn2);
3006  }
3007  }
3008  }
3009  // POIs in negative eta
3010  Dn2 = TwoDiffGapNeg(0,0,iPt).Re();
3011  if(Dn2 != 0)
3012  {
3013  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
3014  {
3015  iHarmonics = fHarmonics[iHarm];
3016  vector = TwoDiffGapNeg(iHarmonics,-iHarmonics,iPt);
3017  dValue = vector.Re()/Dn2;
3018  if( TMath::Abs(dValue < 1) ){
3019  if(fAnalType == kAOD) profile2Neg[iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
3020  if(fAnalType == kMC) profile2Neg[iHarm]->Fill(fIndexImpactParameter, iPt*dPtBinWidth, dValue, Dn2);
3021  }
3022  }
3023  }
3024  }
3025  if(fDoOnlyMixedCorrelations){
3026  // estimating <3'>
3027  // POIs in positive eta
3028  DDn3GapP = ThreeDiffGapPos(0,0,0,iPt).Re();
3029  if(DDn3GapP!=0)
3030  {
3031  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
3032  {
3033  if(iMixedHarm==0){ vector = ThreeDiffGapPos(4,-2,-2,iPt);}
3034  if(iMixedHarm==1){ vector = ThreeDiffGapPos(6,-3,-3,iPt);}
3035  if(iMixedHarm==2){ vector = ThreeDiffGapPos(5,-3,-2,iPt);}
3036  dValue = vector.Re()/DDn3GapP;
3037  if( TMath::Abs(dValue < 1) ){
3038  if(fAnalType == kAOD) profile3Pos[iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
3039  if(fAnalType == kMC) profile3Pos[iMixedHarm]->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn3GapP);
3040  }
3041  }
3042  }
3043  // POIs in negative eta
3044  DDn3GapP = ThreeDiffGapNeg(0,0,0,iPt).Re();
3045  if(DDn3GapP!=0)
3046  {
3047  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
3048  {
3049  if(iMixedHarm==0){ vector = ThreeDiffGapNeg(4,-2,-2,iPt);}
3050  if(iMixedHarm==1){ vector = ThreeDiffGapNeg(6,-3,-3,iPt);}
3051  if(iMixedHarm==2){ vector = ThreeDiffGapNeg(5,-2,-3,iPt);}
3052  dValue = vector.Re()/DDn3GapP;
3053  if( TMath::Abs(dValue < 1) ){
3054  if(fAnalType == kAOD) profile3Neg[iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
3055  if(fAnalType == kMC) profile3Neg[iMixedHarm]->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn3GapP);
3056  }
3057  }
3058  }
3059 
3060  //estimating <4'>
3061  //POIs in positive eta
3062  DDn4GapP = Four13DiffGapPos(0,0,0,0,iPt).Re();
3063  if(DDn4GapP!=0)
3064  {
3065  vector = Four13DiffGapPos(6,-2,-2,-2,iPt);
3066  dValue = vector.Re()/DDn4GapP;
3067  if( TMath::Abs(dValue < 1) ){
3068  if(fAnalType == kAOD) profile4Pos->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
3069  if(fAnalType == kMC) profile4Pos->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn4GapP);
3070  }
3071  }
3072  //POIs in negative eta
3073  DDn4GapP = Four13DiffGapNeg(0,0,0,0,iPt).Re();
3074  if(DDn4GapP!=0)
3075  {
3076  vector = Four13DiffGapNeg(6,-2,-2,-2,iPt);
3077  dValue = vector.Re()/DDn4GapP;
3078  if( TMath::Abs(dValue < 1) ){
3079  if(fAnalType == kAOD) profile4Neg->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
3080  if(fAnalType == kMC) profile4Neg->Fill(fIndexImpactParameter,iPt*dPtBinWidth, dValue, DDn4GapP);
3081  }
3082  }
3083  }
3084  } // endfor {iPt}
3085  return;
3086 }
3087 //_____________________________________________________________________________
3088 //Stopped
3089 
3090 void AliAnalysisTaskFlowModesMC::FillRefsVectors(const Short_t iEtaGapIndex)
3091 {
3092  // Filling Q flow vector with RFPs
3093  // return kTRUE if succesfull (i.e. no error occurs), kFALSE otherwise
3094  // *************************************************************
3095  const Float_t dEtaGap = fEtaGap[iEtaGapIndex];
3096  TH3D* h3NUAWeights = 0x0;
3097  TH1F* hNUEWeights = 0x0;
3098  Double_t dNUAWeight = 1.;
3099  Double_t dNUEWeight = 1.;
3100 
3101  /*if(fAnalType == kAOD){
3102  if(fIndexCentrality>=0 && fIndexCentrality<10){
3103  Int_t c = fIndexCentrality/5;
3104  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
3105  Int_t c = 1 + fIndexCentrality/10;
3106  }
3107  }*/
3108 
3109 
3110 
3111 
3112  // clearing output (global) flow vectors
3113  ResetRFPsVector(fFlowVecQpos);
3114  ResetRFPsVector(fFlowVecQneg);
3115 
3116  Double_t dQcosPos, dQcosNeg, dQsinPos, dQsinNeg;
3117 
3118  for (auto part = fVectorCharged->begin(); part != fVectorCharged->end(); part++)
3119  {
3120  // checking species of used particles (just for double checking purpose)
3121  if( part->species != kCharged)
3122  {
3123  ::Warning("FillRefsVectors","Unexpected part. species (%d) in selected sample (expected %d)",part->species,kCharged);
3124  continue;
3125  }
3126  if(fAnalType == kAOD && fFlowUseNUAWeights)
3127  {
3128  if(part->charge >0) h3NUAWeights = fh3NUAWeightRefsPlus;
3129  if(part->charge <0) h3NUAWeights = fh3NUAWeightRefsMinus;
3130 
3131  if(!h3NUAWeights) { ::Error("FillRefsVectors","Histogram with NUA weights not found."); continue; }
3132  }
3133 
3134  if(fAnalType == kAOD && fFlowUseNUEWeights)
3135  {
3136  if(part->charge >0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightRefsPlus[fIndexCentrality/10];
3137  if(part->charge <0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightRefsMinus[fIndexCentrality/10];
3138 
3139  if(!hNUEWeights) { ::Error("FillRefsVectors","Histogram with NUE weights not found."); return; }
3140  }
3141  // RFPs pT check
3142  if(fCutFlowRFPsPtMin > 0. && part->pt < fCutFlowRFPsPtMin)
3143  continue;
3144 
3145  if(fCutFlowRFPsPtMax > 0. && part->pt > fCutFlowRFPsPtMax)
3146  continue;
3147 
3148  // 0-ing variables
3149  dQcosPos = 0;
3150  dQcosNeg = 0;
3151  dQsinPos = 0;
3152  dQsinNeg = 0;
3153 
3154  // loading weights if needed
3155  if(fAnalType == kAOD && fFlowUseNUAWeights && h3NUAWeights)
3156  {
3157  dNUAWeight = h3NUAWeights->GetBinContent(h3NUAWeights->FindBin(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ()));
3158  if(dNUAWeight <= 0) dNUAWeight = 1.;
3159  if(iEtaGapIndex == 0) fh3AfterNUAWeightsRefs->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3160  }
3161  if(fAnalType == kAOD && fFlowUseNUEWeights && hNUEWeights)
3162  {
3163  dNUEWeight = hNUEWeights->GetBinContent(hNUEWeights->FindBin(part->pt));
3164  if(dNUEWeight <= 0) dNUEWeight = 1.;
3165  if(iEtaGapIndex == 0) fhAfterNUEWeightsRefs->Fill(part->pt, dNUEWeight);
3166  }
3167 
3168  // RPF candidate passing all criteria: start filling flow vectors
3169 
3170  if(part->eta > dEtaGap / 2)
3171  {
3172  // RFP in positive eta acceptance
3173  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3174  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3175  {
3176  dQcosPos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3177  dQsinPos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3178  fFlowVecQpos[iHarm][iPower] += TComplex(dQcosPos,dQsinPos,kFALSE);
3179  }
3180  }
3181  }
3182  if(part->eta < -dEtaGap / 2 )
3183  {
3184  // RFP in negative eta acceptance
3185  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3186  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3187  {
3188  dQcosNeg = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3189  dQsinNeg = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3190  fFlowVecQneg[iHarm][iPower] += TComplex(dQcosNeg,dQsinNeg,kFALSE);
3191  }
3192  }
3193  }
3194  } // endfor {tracks} particle loop
3195 
3196  // // filling local flow vectors to global flow vector arrays
3197  // for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3198  // for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3199  // {
3200  // fFlowVecQpos[iHarm][iPower] = TComplex(dQcosPos[iHarm][iPower],dQsinPos[iHarm][iPower],kFALSE);
3201  // if(dEtaGap > -1)
3202  // fFlowVecQneg[iHarm][iPower] = TComplex(dQcosNeg[iHarm][iPower],dQsinNeg[iHarm][iPower],kFALSE);
3203  // }
3204 
3205  // printf("RFPs EtaGap %g : number %g (pos) %g (neg) \n", dEtaGap,fFlowVecQpos[0][0].Re(),fFlowVecQneg[0][0].Re());
3206  return;
3207 }
3208 //_____________________________________________________________________________
3209 void AliAnalysisTaskFlowModesMC::FillPOIsVectors(const Short_t iEtaGapIndex, const PartSpecies species, const Short_t iMassIndex)
3210 {
3211  // Filling p,q and s flow vectors with POIs (given by species) for differential flow calculation
3212  // *************************************************************
3213 
3214  if(species == kUnknown) return;
3215 
3216  Double_t dNUAWeight = 1.; // for generic framework != 1
3217  Double_t dNUEWeight = 1.; // for generic framework != 1
3218 
3219  // clearing output (global) flow vectors
3220  ResetPOIsVector(fFlowVecPpos);
3221  ResetPOIsVector(fFlowVecPneg);
3222  ResetPOIsVector(fFlowVecS);
3223 
3224  std::vector<FlowPart>* vector = 0x0;
3225  //TH3D* hist = 0x0;
3226  //Double_t dMassLow = 0, dMassHigh = 0;
3227  TH3D* h3NUAWeights = 0x0;
3228  TH1F* hNUEWeights = 0x0;
3229 
3230  /*if(fAnalType == kAOD){
3231  if(fIndexCentrality>=0 && fIndexCentrality<10){
3232  Int_t c = fIndexCentrality/5;
3233  }else if(fIndexCentrality>=10 && fIndexCentrality<50){
3234  Int_t c = 1 + fIndexCentrality/10;
3235  }
3236  }*/
3237 
3238  // swich based on species
3239  switch (species)
3240  {
3241  case kCharged:
3242  vector = fVectorCharged;
3243  break;
3244 
3245  case kPion:
3246  vector = fVectorPion;
3247  break;
3248 
3249  case kKaon:
3250  vector = fVectorKaon;
3251  break;
3252 
3253  case kProton:
3254  vector = fVectorProton;
3255  break;
3256 
3257  default:
3258  ::Error("FillPOIsVectors","Selected species unknown.");
3259  return;
3260  }//switch(species)
3261 
3262  const Double_t dEtaGap = fEtaGap[iEtaGapIndex];
3263  //const Double_t dMass = (dMassLow+dMassHigh)/2;
3264 
3265  Short_t iPtBin = 0;
3266  Double_t dCos = 0, dSin = 0;
3267 
3268  for (auto part = vector->begin(); part != vector->end(); part++)
3269  {
3270  if(fAnalType == kAOD){
3271  // swich based on species
3272  switch (species)
3273  {
3274  case kCharged:
3275  if(fFlowUseNUAWeights) {
3276  if(part->charge >0) h3NUAWeights = fh3NUAWeightChargedPlus;
3277  if(part->charge <0) h3NUAWeights = fh3NUAWeightChargedMinus;
3278  }
3279  if(fFlowUseNUEWeights) {
3280  if(part->charge >0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightChargedPlus[fIndexCentrality/10];
3281  if(part->charge <0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightChargedMinus[fIndexCentrality/10];
3282  }
3283  break;
3284 
3285  case kPion:
3286  if(fFlowUseNUAWeights) {
3287  if(part->charge >0) h3NUAWeights = fh3NUAWeightPionPlus;
3288  if(part->charge <0) h3NUAWeights = fh3NUAWeightPionMinus;
3289  }
3290  if(fFlowUseNUEWeights) {
3291  if(part->charge >0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightPionPlus[fIndexCentrality/10];
3292  if(part->charge <0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightPionMinus[fIndexCentrality/10];
3293  }
3294  break;
3295 
3296  case kKaon:
3297  if(fFlowUseNUAWeights) {
3298  if(part->charge >0) h3NUAWeights = fh3NUAWeightKaonPlus;
3299  if(part->charge <0) h3NUAWeights = fh3NUAWeightKaonMinus;
3300  }
3301  if(fFlowUseNUEWeights) {
3302  if(part->charge >0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightKaonPlus[fIndexCentrality/10];
3303  if(part->charge <0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightKaonMinus[fIndexCentrality/10];
3304  }
3305  break;
3306 
3307  case kProton:
3308  if(fFlowUseNUAWeights) {
3309  if(part->charge >0) h3NUAWeights = fh3NUAWeightProtonPlus;
3310  if(part->charge <0) h3NUAWeights = fh3NUAWeightProtonMinus;
3311  }
3312  if(fFlowUseNUEWeights) {
3313  if(part->charge >0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightProtonPlus[fIndexCentrality/10];
3314  if(part->charge <0 && fIndexCentrality<50) hNUEWeights = fhNUEWeightProtonMinus[fIndexCentrality/10];
3315  }
3316  break;
3317 
3318  default:
3319  ::Error("FillPOIsVectors","Selected species unknown.");
3320  return;
3321  }//switch(species)
3322 
3323  if(fFlowUseNUAWeights && !h3NUAWeights) { ::Error("FillPOIsVectors","Histogram with NUA weights not found."); continue; }
3324  if(fFlowUseNUEWeights && !hNUEWeights) { ::Error("FillPOIsVectors","Histogram with NUE weights not found."); return; }
3325  }
3326  // checking species of used particles (just for double checking purpose)
3327  if( part->species != species)
3328  {
3329  ::Warning("FillPOIsVectors","Unexpected part. species (%d) in selected sample (expected %d)",part->species,species);
3330  continue;
3331  }//endif{part->species != species}
3332 
3333  // assign iPtBin based on particle momenta
3334  iPtBin = GetPOIsPtBinIndex(part->pt);
3335 
3336  // 0-ing variables
3337  dCos = 0;
3338  dSin = 0;
3339 
3340  // POIs candidate passing all criteria: start filling flow vectors
3341 
3342  // loading weights if needed
3343  if(fAnalType == kAOD && fFlowUseNUAWeights && h3NUAWeights)
3344  {
3345  dNUAWeight = h3NUAWeights->GetBinContent(h3NUAWeights->FindBin(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ()));
3346  if(dNUAWeight <= 0){ dNUAWeight = 1.;}
3347  if(iEtaGapIndex == 0){
3348  switch (species)
3349  {
3350  case kCharged:
3351  fh3AfterNUAWeightsCharged->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3352  break;
3353  case kPion:
3354  fh3AfterNUAWeightsPion->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3355  break;
3356  case kKaon:
3357  fh3AfterNUAWeightsKaon->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3358  break;
3359  case kProton:
3360  fh3AfterNUAWeightsProton->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3361  break;
3362  default:
3363  ::Error("Fill fh3AfterNUAWeights","Selected species unknown.");
3364  return;
3365  }
3366  }
3367  }//endif{fFlowUseNUAWeights && h3NUAWeights}
3368 
3369  if(fAnalType == kAOD && fFlowUseNUEWeights && hNUEWeights)
3370  {
3371  dNUEWeight = hNUEWeights->GetBinContent(hNUEWeights->FindBin(part->pt));
3372  if(dNUEWeight <= 0) dNUEWeight = 1.;
3373  if(iEtaGapIndex == 0){
3374  switch (species)
3375  {
3376  case kCharged:
3377  fhAfterNUEWeightsCharged->Fill(part->pt, dNUEWeight);
3378  break;
3379  case kPion:
3380  fhAfterNUEWeightsPion->Fill(part->pt, dNUEWeight);
3381  break;
3382  case kKaon:
3383  fhAfterNUEWeightsKaon->Fill(part->pt, dNUEWeight);
3384  break;
3385  case kProton:
3386  fhAfterNUEWeightsProton->Fill(part->pt, dNUEWeight);
3387  break;
3388  default:
3389  ::Error("Fill fhAfterNUEWeights","Selected species unknown.");
3390  return;
3391  }
3392  }
3393  }//endif{fFlowUseNUEWeights && hNUEWeights}
3394 
3395  if(part->eta > dEtaGap / 2 )
3396  {
3397  // particle in positive eta acceptance
3398  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3399  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3400  {
3401  dCos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3402  dSin = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3403  fFlowVecPpos[iHarm][iPower][iPtBin] += TComplex(dCos,dSin,kFALSE);
3404  }//endfor{Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++}
3405  }//endfor{Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++}
3406  }//endif{part->eta > dEtaGap / 2 }
3407  if(part->eta < -dEtaGap / 2 ){
3408  // particle in negative eta acceptance
3409  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3410  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3411  {
3412  dCos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3413  dSin = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3414  fFlowVecPneg[iHarm][iPower][iPtBin] += TComplex(dCos,dSin,kFALSE);
3415  }//endfor{Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++}
3416  }//endfor{Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++}
3417  }//endif{part->eta < -dEtaGap / 2 }
3418  } // endfor {tracks}
3419  return;
3420 }
3421 //_____________________________________________________________________________
3422 Short_t AliAnalysisTaskFlowModesMC::GetPOIsPtBinIndex(const Double_t pt)
3423 {
3424  // Return POIs pT bin index based on pT value
3425  // *************************************************************
3426  const Double_t dPtBinWidth = (fFlowPOIsPtMax - fFlowPOIsPtMin) / fFlowPOIsPtNumBins;
3427  // printf("Pt %g | index %d\n",pt,(Short_t) (pt / dPtBinWidth) );
3428  return (Short_t) (pt / dPtBinWidth);
3429 }
3430 //_____________________________________________________________________________
3431 void AliAnalysisTaskFlowModesMC::ResetRFPsVector(TComplex (&array)[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax])
3432 {
3433  // Reset RFPs (Q) array values to TComplex(0,0,kFALSE) for given array
3434  // *************************************************************
3435  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3436  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3437  array[iHarm][iPower] = TComplex(0,0,kFALSE);
3438  return;
3439 }
3440 //_____________________________________________________________________________
3441 void AliAnalysisTaskFlowModesMC::ResetPOIsVector(TComplex (&array)[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax][fFlowPOIsPtNumBins])
3442 {
3443  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3444  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3445  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
3446  array[iHarm][iPower][iPt] = TComplex(0,0,kFALSE);
3447  return;
3448 }
3449 //_____________________________________________________________________________
3450 void AliAnalysisTaskFlowModesMC::ListFlowVector(TComplex (&array)[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax])
3451 {
3452  // List all values of given flow vector TComplex array
3453  // *************************************************************
3454  printf(" ### Listing (TComplex) flow vector array ###########################\n");
3455  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3456  {
3457  printf("Harm %d (power):",iHarm);
3458  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3459  {
3460  printf("|(%d) %g+%g(i)",iPower, array[iHarm][iPower].Re(), array[iHarm][iPower].Im());
3461  }
3462  printf("\n");
3463  }
3464  return;
3465 }
3466 //_____________________________________________________________________________
3467 Short_t AliAnalysisTaskFlowModesMC::GetSamplingIndex()
3468 {
3469  // Assessing sampling index based on generated random number
3470  // returns centrality index
3471  // *************************************************************
3472 
3473  Short_t index = 0x0;
3474 
3475  if(fSampling && fNumSamples > 1)
3476  {
3477  TRandom3 rr(0);
3478  Double_t ranNum = rr.Rndm(); // getting random number in (0,1)
3479  Double_t generated = ranNum * fNumSamples; // getting random number in range (0, fNumSamples)
3480  // finding right index for sampling based on generated number and total number of samples
3481  for(Short_t i(0); i < fNumSamples; i++)
3482  {
3483  if(generated < (i+1) )
3484  {
3485  index = i;
3486  break;
3487  }
3488  }
3489  }
3490 
3491  return index;
3492 }
3493 
3494 //_____________________________________________________________________________
3495 Short_t AliAnalysisTaskFlowModesMC::GetCentralityIndex()
3496 {
3497  // Estimating centrality percentile based on selected estimator.
3498  // (Default) If no multiplicity estimator is specified (fMultEstimator == '' || Charged), percentile is estimated as number of selected / filtered charged tracks.
3499  // If a valid multiplicity estimator is specified, centrality percentile is estimated via AliMultSelection
3500  // otherwise -1 is returned (and event is skipped)
3501  // *************************************************************
3502  if(fAnalType == kAOD){
3503  fMultEstimator.ToUpper();
3504  //cout<<"++++++++++++++++++GetCentralityIndex+++++++++++++++++"<<endl;
3505  if(
3506  fMultEstimator.EqualTo("V0A") || fMultEstimator.EqualTo("V0C") || fMultEstimator.EqualTo("V0M") ||
3507  fMultEstimator.EqualTo("CL0") || fMultEstimator.EqualTo("CL1") || fMultEstimator.EqualTo("ZNA") ||
3508  fMultEstimator.EqualTo("ZNC") || fMultEstimator.EqualTo("TRK") ){
3509  //cout<<"++++++++++++++++++CentralityIndex+++++++++++++++++"<<endl;
3510 
3511  // some of supported AliMultSelection estimators (listed above)
3512  Float_t dPercentile = 300;
3513 
3514  // checking AliMultSelection
3515  AliMultSelection* multSelection = 0x0;
3516  multSelection = (AliMultSelection*) fEventAOD->FindListObject("MultSelection");
3517  if(!multSelection) { AliError("AliMultSelection object not found! Returning -1"); return -1;}
3518 
3519  dPercentile = multSelection->GetMultiplicityPercentile(fMultEstimator.Data());
3520  //cout<<"fMultEstimator.Data()"<<fMultEstimator.Data()<<endl;
3521  //cout<<"dPercentile = "<<dPercentile<<endl;
3522  if(fFullCentralityRange){
3523  if(dPercentile > 100 || dPercentile < 0) { AliWarning("Centrality percentile estimated not within 0-100 range. Returning -1"); return -1;}
3524  else { return dPercentile;}
3525  }
3526  if(!fFullCentralityRange){
3527  if(dPercentile > 100 || dPercentile <50) { AliWarning("Centrality percentile estimated not within 50-100 range. Returning -1"); return -1;}
3528  else { return dPercentile;}
3529  }
3530  }
3531  else if(fMultEstimator.EqualTo("") || fMultEstimator.EqualTo("CHARGED"))
3532  {
3533  // assigning centrality based on number of selected charged tracks
3534  return fVectorCharged->size();
3535  }
3536  else
3537  {
3538  AliWarning(Form("Multiplicity estimator '%s' not supported. Returning -1\n",fMultEstimator.Data()));
3539  return -1;
3540  }
3541  }
3542  return -1;
3543 }
3544 //_____________________________________________________________________________
3545 Double_t AliAnalysisTaskFlowModesMC::GetImpactParameter()
3546 {
3547  if(fAnalType == kMC){
3548  if(fCutImpactParameter) {
3549  Double_t gImpactParameter = 0.;
3550  if(fEventMC) {
3551  AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(fEventMC)->GenEventHeader());
3552  if(headerH) gImpactParameter = headerH->ImpactParameter();
3553  }
3554  if ((gImpactParameter > fImpactParameterMin) && (gImpactParameter < fImpactParameterMax )) return gImpactParameter;
3555  if ((gImpactParameter < fImpactParameterMin) || (gImpactParameter >= fImpactParameterMax )) return -1;
3556  }
3557  }
3558  return -1;
3559 }
3560 //____________________________
3561 Double_t AliAnalysisTaskFlowModesMC::GetWDist(const AliAODVertex* v0, const AliAODVertex* v1) {
3562  // calculate sqrt of weighted distance to other vertex
3563  if (!v0 || !v1) {
3564  printf("One of vertices is not valid\n");
3565  return kFALSE;
3566  }
3567  static TMatrixDSym vVb(3);
3568  double dist = -1;
3569  double dx = v0->GetX()-v1->GetX();
3570  double dy = v0->GetY()-v1->GetY();
3571  double dz = v0->GetZ()-v1->GetZ();
3572  double cov0[6],cov1[6];
3573  v0->GetCovarianceMatrix(cov0);
3574  v1->GetCovarianceMatrix(cov1);
3575  vVb(0,0) = cov0[0]+cov1[0];
3576  vVb(1,1) = cov0[2]+cov1[2];
3577  vVb(2,2) = cov0[5]+cov1[5];
3578  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
3579  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
3580  vVb.InvertFast();
3581  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
3582  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz+2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
3583  return dist>0 ? TMath::Sqrt(dist) : -1;
3584 }
3585 //_________________________________________________
3586 void AliAnalysisTaskFlowModesMC::Terminate(Option_t* option)
3587 {
3588  // called on end of task, after all events are processed
3589  // *************************************************************
3590 
3591  return;
3592 }
3593 //_____________________________________________________________________________
3594 // Set of methods returning given complex flow vector based on flow harmonics (n) and weight power indexes (p)
3595 // a la General Framework implementation.
3596 // Q: flow vector of RFPs (with/out eta gap)
3597 // P: flow vector of POIs (with/out eta gap) (in usual notation p)
3598 // S: flow vector of overlaping RFPs and POIs (in usual notation q)
3599 
3600 TComplex AliAnalysisTaskFlowModesMC::Q(const Short_t n, const Short_t p)
3601 {
3602  if (n < 0) return TComplex::Conjugate(fFlowVecQpos[-n][p]);
3603  else return fFlowVecQpos[n][p];
3604 }
3605 //____________________________________________________________________
3606 TComplex AliAnalysisTaskFlowModesMC::QGapPos(const Short_t n, const Short_t p)
3607 {
3608  if (n < 0) return TComplex::Conjugate(fFlowVecQpos[-n][p]);
3609  else return fFlowVecQpos[n][p];
3610 }
3611 //____________________________________________________________________
3612 TComplex AliAnalysisTaskFlowModesMC::QGapNeg(const Short_t n, const Short_t p)
3613 {
3614  if(n < 0) return TComplex::Conjugate(fFlowVecQneg[-n][p]);
3615  else return fFlowVecQneg[n][p];
3616 }
3617 //____________________________________________________________________
3618 TComplex AliAnalysisTaskFlowModesMC::P(const Short_t n, const Short_t p, const Short_t pt)
3619 {
3620  if(n < 0) return TComplex::Conjugate(fFlowVecPpos[-n][p][pt]);
3621  else return fFlowVecPpos[n][p][pt];
3622 }
3623 //____________________________________________________________________
3624 TComplex AliAnalysisTaskFlowModesMC::PGapPos(const Short_t n, const Short_t p, const Short_t pt)
3625 {
3626  if(n < 0) return TComplex::Conjugate(fFlowVecPpos[-n][p][pt]);
3627  else return fFlowVecPpos[n][p][pt];
3628 }
3629 //____________________________________________________________________
3630 TComplex AliAnalysisTaskFlowModesMC::PGapNeg(const Short_t n, const Short_t p, const Short_t pt)
3631 {
3632  if(n < 0) return TComplex::Conjugate(fFlowVecPneg[-n][p][pt]);
3633  else return fFlowVecPneg[n][p][pt];
3634 }
3635 //____________________________________________________________________
3636 TComplex AliAnalysisTaskFlowModesMC::S(const Short_t n, const Short_t p, const Short_t pt)
3637 {
3638  if(n < 0) return TComplex::Conjugate(fFlowVecS[-n][p][pt]);
3639  else return fFlowVecS[n][p][pt];
3640 }
3641 //____________________________________________________________________
3642 
3643 // Set of flow calculation methods for cumulants of different orders with/out eta gap
3644 
3645 TComplex AliAnalysisTaskFlowModesMC::Two(const Short_t n1, const Short_t n2)
3646 {
3647  TComplex formula = Q(n1,1)*Q(n2,1) - Q(n1+n2,2);
3648  return formula;
3649 }
3650 //____________________________________________________________________
3651 TComplex AliAnalysisTaskFlowModesMC::TwoGap(const Short_t n1, const Short_t n2)
3652 {
3653  TComplex formula = QGapPos(n1,1)*QGapNeg(n2,1);
3654  return formula;
3655 }
3656 //____________________________________________________________________
3657 TComplex AliAnalysisTaskFlowModesMC::TwoDiff(const Short_t n1, const Short_t n2, const Short_t pt)
3658 {
3659  TComplex formula = P(n1,1,pt)*Q(n2,1) - S(n1+n2,1,pt);
3660  return formula;
3661 }
3662 //____________________________________________________________________
3663 TComplex AliAnalysisTaskFlowModesMC::TwoDiffGapPos(const Short_t n1, const Short_t n2, const Short_t pt)
3664 {
3665  TComplex formula = PGapPos(n1,1,pt)*QGapNeg(n2,1);
3666  return formula;
3667 }
3668 //____________________________________________________________________
3669 TComplex AliAnalysisTaskFlowModesMC::TwoDiffGapNeg(const Short_t n1, const Short_t n2, const Short_t pt)
3670 {
3671  TComplex formula = PGapNeg(n1,1,pt)*QGapPos(n2,1);
3672  return formula;
3673 }
3674 //____________________________________________________________________
3675 TComplex AliAnalysisTaskFlowModesMC::ThreeDiffGapPos(const Short_t n1, const Short_t n2, const Short_t n3,const Short_t pt)
3676 {
3677  TComplex formula = PGapPos(n1,1,pt)*QGapNeg(n2,1)*QGapNeg(n3,1)- PGapPos(n1,1,pt)*QGapNeg(n2+n3,2);
3678  return formula;
3679 }
3680 //____________________________________________________________________
3681 TComplex AliAnalysisTaskFlowModesMC::ThreeDiffGapNeg(const Short_t n1, const Short_t n2, const Short_t n3,const Short_t pt)
3682 {
3683  TComplex formula = PGapNeg(n1,1,pt)*QGapPos(n2,1)*QGapPos(n3,1)- PGapNeg(n1,1,pt)*QGapPos(n2+n3,2);
3684  return formula;
3685 }
3686 //____________________________________________________________________
3687 TComplex AliAnalysisTaskFlowModesMC::Four(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4)
3688 {
3689  TComplex formula = Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4,1)-Q(n1+n2,2)*Q(n3,1)*Q(n4,1)-Q(n2,1)*Q(n1+n3,2)*Q(n4,1)
3690  - Q(n1,1)*Q(n2+n3,2)*Q(n4,1)+2.*Q(n1+n2+n3,3)*Q(n4,1)-Q(n2,1)*Q(n3,1)*Q(n1+n4,2)
3691  + Q(n2+n3,2)*Q(n1+n4,2)-Q(n1,1)*Q(n3,1)*Q(n2+n4,2)+Q(n1+n3,2)*Q(n2+n4,2)
3692  + 2.*Q(n3,1)*Q(n1+n2+n4,3)-Q(n1,1)*Q(n2,1)*Q(n3+n4,2)+Q(n1+n2,2)*Q(n3+n4,2)
3693  + 2.*Q(n2,1)*Q(n1+n3+n4,3)+2.*Q(n1,1)*Q(n2+n3+n4,3)-6.*Q(n1+n2+n3+n4,4);
3694  return formula;
3695 }
3696 // //____________________________________________________________________
3697 TComplex AliAnalysisTaskFlowModesMC::FourDiff(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t pt)
3698 {
3699  TComplex formula = P(n1,1,pt)*Q(n2,1)*Q(n3,1)*Q(n4,1)-S(n1+n2,2,pt)*Q(n3,1)*Q(n4,1)-Q(n2,1)*S(n1+n3,2,pt)*Q(n4,1)
3700  - P(n1,1,pt)*Q(n2+n3,2)*Q(n4,1)+2.*S(n1+n2+n3,3,pt)*Q(n4,1)-Q(n2,1)*Q(n3,1)*S(n1+n4,2,pt)
3701  + Q(n2+n3,2)*S(n1+n4,2,pt)-P(n1,1,pt)*Q(n3,1)*Q(n2+n4,2)+S(n1+n3,2,pt)*Q(n2+n4,2)
3702  + 2.*Q(n3,1)*S(n1+n2+n4,3,pt)-P(n1,1,pt)*Q(n2,1)*Q(n3+n4,2)+S(n1+n2,2,pt)*Q(n3+n4,2)
3703  + 2.*Q(n2,1)*S(n1+n3+n4,3,pt)+2.*P(n1,1,pt)*Q(n2+n3+n4,3)-6.*S(n1+n2+n3+n4,4,pt);
3704  return formula;
3705 }
3706 //____________________________________________________________________
3707 TComplex AliAnalysisTaskFlowModesMC::FourGapPos(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4) // n1+n2 = n3+n4; n1, n2 from P & n3, n4 from M
3708 {
3709  TComplex formula = QGapPos(n1,1)*QGapPos(n2,1)*QGapNeg(n3,1)*QGapNeg(n4,1)-QGapPos(n1+n2,2)*QGapNeg(n3,1)*QGapNeg(n4,1)-QGapPos(n1,1)*QGapPos(n2,1)*QGapNeg(n3+n4,2)+QGapPos(n1+n2,2)*QGapNeg(n3+n4,2);
3710  //TComplex *out = (TComplex*) &formula;
3711  return formula;
3712 }
3713 //____________________________________________________________________
3714 TComplex AliAnalysisTaskFlowModesMC::FourGapNeg(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4) // n1+n2 = n3+n4; n1, n2 from M & n3, n4 from P
3715 {
3716  TComplex formula = QGapNeg(n1,1)*QGapNeg(n2,1)*QGapPos(n3,1)*QGapPos(n4,1)-QGapNeg(n1+n2,2)*QGapPos(n3,1)*QGapPos(n4,1)-QGapNeg(n1,1)*QGapNeg(n2,1)*QGapPos(n3+n4,2)+QGapNeg(n1+n2,2)*QGapPos(n3+n4,2);
3717  //TComplex *out = (TComplex*) &formula;
3718  return formula;
3719 }
3720 //____________________________________________________________________
3721 TComplex AliAnalysisTaskFlowModesMC::Four13DiffGapPos(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4,const Short_t pt) // n1 = n2 + n3 + n4
3722 {
3723  TComplex formula = PGapPos(n1,1,pt)*QGapNeg(n2,1)*QGapNeg(n3,1)*QGapNeg(n4,1)- PGapPos(n1,1,pt)*QGapNeg(n2+n3,2)*QGapNeg(n4,1) - PGapPos(n1,1,pt)*QGapNeg(n2+n4,2)*QGapNeg(n3,1) - PGapPos(n1,1,pt)*QGapNeg(n3+n4,2)*QGapNeg(n2,1) + 2.*PGapPos(n1,1,pt)*QGapNeg(n2+n3+n4,3);
3724  return formula;
3725 }
3726 //____________________________________________________________________
3727 TComplex AliAnalysisTaskFlowModesMC::Four13DiffGapNeg(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4,const Short_t pt) // n1 = n2 + n3 + n4
3728 {
3729  TComplex formula = PGapNeg(n1,1,pt)*QGapPos(n2,1)*QGapPos(n3,1)*QGapPos(n4,1)- PGapNeg(n1,1,pt)*QGapPos(n2+n3,2)*QGapPos(n4,1) - PGapNeg(n1,1,pt)*QGapPos(n2+n4,2)*QGapPos(n3,1) - PGapNeg(n1,1,pt)*QGapPos(n3+n4,2)*QGapPos(n2,1) + 2.*PGapNeg(n1,1,pt)*QGapPos(n2+n3+n4,3);
3730  return formula;
3731 }
3732 //____________________________________________________________________
3733 TComplex AliAnalysisTaskFlowModesMC::SixGapPos(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t n5, const Short_t n6) // n1 + n2 + n3 = n4 + n5 + n6
3734 {
3735  TComplex formula = QGapPos(n1,1)*QGapPos(n2,1)*QGapPos(n3,1)*QGapNeg(n4,1)*QGapNeg(n5,1)*QGapNeg(n6,1) - QGapPos(n1,1)*QGapPos(n2,1)*QGapPos(n3,1)*QGapNeg(n4+n5,2)*QGapNeg(n6,1) - QGapPos(n1,1)*QGapPos(n2,1)*QGapPos(n3,1)*QGapNeg(n4+n6,2)*QGapNeg(n5,1) - QGapPos(n1,1)*QGapPos(n2,1)*QGapPos(n3,1)*QGapNeg(n5+n6,2)*QGapNeg(n4,1) + 2.*QGapPos(n1,1)*QGapPos(n2,1)*QGapPos(n3,1)*QGapNeg(n4+n5+n6,3) - QGapPos(n1+n2,2)*QGapPos(n3,1)*QGapNeg(n4,1)*QGapNeg(n5,1)*QGapNeg(n6,1) + QGapPos(n1+n2,2)*QGapPos(n3,1)*QGapNeg(n4+n5,2)*QGapNeg(n6,1) + QGapPos(n1+n2,2)*QGapPos(n3,1)*QGapNeg(n4+n6,2)*QGapNeg(n5,1) + QGapPos(n1+n2,2)*QGapPos(n3,1)*QGapNeg(n5+n6,2)*QGapNeg(n4,1) - 2.*QGapPos(n1+n2,2)*QGapPos(n3,1)*QGapNeg(n4+n5+n6,3) - QGapPos(n1+n3,2)*QGapPos(n2,1)*QGapNeg(n4,1)*QGapNeg(n5,1)*QGapNeg(n6,1) + QGapPos(n1+n3,2)*QGapPos(n2,1)*QGapNeg(n4+n5,2)*QGapNeg(n6,1) + QGapPos(n1+n3,2)*QGapPos(n2,1)*QGapNeg(n4+n6,2)*QGapNeg(n5,1) + QGapPos(n1+n3,2)*QGapPos(n2,1)*QGapNeg(n5+n6,2)*QGapNeg(n4,1) - 2.*QGapPos(n1+n3,2)*QGapPos(n2,1)*QGapNeg(n4+n5+n6,3) - QGapPos(n2+n3,2)*QGapPos(n1,1)*QGapNeg(n4,1)*QGapNeg(n5,1)*QGapNeg(n6,1) + QGapPos(n2+n3,2)*QGapPos(n1,1)*QGapNeg(n4+n5,2)*QGapNeg(n6,1) + QGapPos(n2+n3,2)*QGapPos(n1,1)*QGapNeg(n4+n6,2)*QGapNeg(n5,1) + QGapPos(n2+n3,2)*QGapPos(n1,1)*QGapNeg(n5+n6,2)*QGapNeg(n4,1) - 2.*QGapPos(n2+n3,2)*QGapPos(n1,1)*QGapNeg(n4+n5+n6,3)+ 2.*QGapPos(n1+n2+n3,3)*QGapNeg(n4,1)*QGapNeg(n5,1)*QGapNeg(n6,1) - 2.*QGapPos(n1+n2+n3,3)*QGapNeg(n4+n5,2)*QGapNeg(n6,1) - 2.*QGapPos(n1+n2+n3,3)*QGapNeg(n4+n6,2)*QGapNeg(n5,1) - 2.*QGapPos(n1+n2+n3,3)*QGapNeg(n5+n6,2)*QGapNeg(n4,1) + 4.*QGapPos(n1+n2+n3,3)*QGapNeg(n4+n5+n6,3);
3736  return formula;
3737 }
3738 //____________________________________________________________________
3739 TComplex AliAnalysisTaskFlowModesMC::SixGapNeg(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t n5, const Short_t n6) // n1 + n2 + n3 = n4 + n5 + n6
3740 {
3741 
3742  TComplex formula = QGapNeg(n1,1)*QGapNeg(n2,1)*QGapNeg(n3,1)*QGapPos(n4,1)*QGapPos(n5,1)*QGapPos(n6,1) - QGapNeg(n1,1)*QGapNeg(n2,1)*QGapNeg(n3,1)*QGapPos(n4+n5,2)*QGapPos(n6,1) - QGapNeg(n1,1)*QGapNeg(n2,1)*QGapNeg(n3,1)*QGapPos(n4+n6,2)*QGapPos(n5,1) - QGapNeg(n1,1)*QGapNeg(n2,1)*QGapNeg(n3,1)*QGapPos(n5+n6,2)*QGapPos(n4,1) + 2.*QGapNeg(n1,1)*QGapNeg(n2,1)*QGapNeg(n3,1)*QGapPos(n4+n5+n6,3)- QGapNeg(n1+n2,2)*QGapNeg(n3,1)*QGapPos(n4,1)*QGapPos(n5,1)*QGapPos(n6,1) + QGapNeg(n1+n2,2)*QGapNeg(n3,1)*QGapPos(n4+n5,2)*QGapPos(n6,1) + QGapNeg(n1+n2,2)*QGapNeg(n3,1)*QGapPos(n4+n6,2)*QGapPos(n5,1) + QGapNeg(n1+n2,2)*QGapNeg(n3,1)*QGapPos(n5+n6,2)*QGapPos(n4,1) - 2.*QGapNeg(n1+n2,2)*QGapNeg(n3,1)*QGapPos(n4+n5+n6,3) - QGapNeg(n1+n3,2)*QGapNeg(n2,1)*QGapPos(n4,1)*QGapPos(n5,1)*QGapPos(n6,1) + QGapNeg(n1+n3,2)*QGapNeg(n2,1)*QGapPos(n4+n5,2)*QGapPos(n6,1) + QGapNeg(n1+n3,2)*QGapNeg(n2,1)*QGapPos(n4+n6,2)*QGapPos(n5,1) + QGapNeg(n1+n3,2)*QGapNeg(n2,1)*QGapPos(n5+n6,2)*QGapPos(n4,1) - 2.*QGapNeg(n1+n3,2)*QGapNeg(n2,1)*QGapPos(n4+n5+n6,3) - QGapNeg(n2+n3,2)*QGapNeg(n1,1)*QGapPos(n4,1)*QGapPos(n5,1)*QGapPos(n6,1) + QGapNeg(n2+n3,2)*QGapNeg(n1,1)*QGapPos(n4+n5,2)*QGapPos(n6,1) + QGapNeg(n2+n3,2)*QGapNeg(n1,1)*QGapPos(n4+n6,2)*QGapPos(n5,1) + QGapNeg(n2+n3,2)*QGapNeg(n1,1)*QGapPos(n5+n6,2)*QGapPos(n4,1) - 2.*QGapNeg(n2+n3,2)*QGapNeg(n1,1)*QGapPos(n4+n5+n6,3)+ 2.*QGapNeg(n1+n2+n3,3)*QGapPos(n4,1)*QGapPos(n5,1)*QGapPos(n6,1) -2.*QGapNeg(n1+n2+n3,3)*QGapPos(n4+n5,2)*QGapPos(n6,1) - 2.*QGapNeg(n1+n2+n3,3)*QGapPos(n4+n6,2)*QGapPos(n5,1) - 2.*QGapNeg(n1+n2+n3,3)*QGapPos(n5+n6,2)*QGapPos(n4,1) + 4.*QGapNeg(n1+n2+n3,3)*QGapPos(n4+n5+n6,3);
3743  return formula;
3744 }
3745 
3746 //-----------------------------------------------------------------------
3747 void AliAnalysisTaskFlowModesMC::SetPriors(Float_t centrCur){
3748  //set priors for the bayesian pid selection
3749  fCurrCentr = centrCur;
3750 
3751  fBinLimitPID[0] = 0.300000;
3752  fBinLimitPID[1] = 0.400000;
3753  fBinLimitPID[2] = 0.500000;
3754  fBinLimitPID[3] = 0.600000;
3755  fBinLimitPID[4] = 0.700000;
3756  fBinLimitPID[5] = 0.800000;
3757  fBinLimitPID[6] = 0.900000;
3758  fBinLimitPID[7] = 1.000000;
3759  fBinLimitPID[8] = 1.200000;
3760  fBinLimitPID[9] = 1.400000;
3761  fBinLimitPID[10] = 1.600000;
3762  fBinLimitPID[11] = 1.800000;
3763  fBinLimitPID[12] = 2.000000;
3764  fBinLimitPID[13] = 2.200000;
3765  fBinLimitPID[14] = 2.400000;
3766  fBinLimitPID[15] = 2.600000;
3767  fBinLimitPID[16] = 2.800000;
3768  fBinLimitPID[17] = 3.000000;
3769 
3770  for(Int_t i=18;i<fgkPIDptBin;i++){
3771  fBinLimitPID[i] = 3.0 + 0.2 * (i-17);
3772  }
3773 
3774  // 0-10%
3775  if(centrCur < 10){
3776  fC[0][0] = 0.005;
3777  fC[0][1] = 0.005;
3778  fC[0][2] = 1.0000;
3779  fC[0][3] = 0.0010;
3780  fC[0][4] = 0.0010;
3781 
3782  fC[1][0] = 0.005;
3783  fC[1][1] = 0.005;
3784  fC[1][2] = 1.0000;
3785  fC[1][3] = 0.0168;
3786  fC[1][4] = 0.0122;
3787 
3788  fC[2][0] = 0.005;
3789  fC[2][1] = 0.005;
3790  fC[2][2] = 1.0000;
3791  fC[2][3] = 0.0272;
3792  fC[2][4] = 0.0070;
3793 
3794  fC[3][0] = 0.005;
3795  fC[3][1] = 0.005;
3796  fC[3][2] = 1.0000;
3797  fC[3][3] = 0.0562;
3798  fC[3][4] = 0.0258;
3799 
3800  fC[4][0] = 0.005;
3801  fC[4][1] = 0.005;
3802  fC[4][2] = 1.0000;
3803  fC[4][3] = 0.0861;
3804  fC[4][4] = 0.0496;
3805 
3806  fC[5][0] = 0.005;
3807  fC[5][1] = 0.005;
3808  fC[5][2] = 1.0000;
3809  fC[5][3] = 0.1168;
3810  fC[5][4] = 0.0740;
3811 
3812  fC[6][0] = 0.005;
3813  fC[6][1] = 0.005;
3814  fC[6][2] = 1.0000;
3815  fC[6][3] = 0.1476;
3816  fC[6][4] = 0.0998;
3817 
3818  fC[7][0] = 0.005;
3819  fC[7][1] = 0.005;
3820  fC[7][2] = 1.0000;
3821  fC[7][3] = 0.1810;
3822  fC[7][4] = 0.1296;
3823 
3824  fC[8][0] = 0.005;
3825  fC[8][1] = 0.005;
3826  fC[8][2] = 1.0000;
3827  fC[8][3] = 0.2240;
3828  fC[8][4] = 0.1827;
3829 
3830  fC[9][0] = 0.005;
3831  fC[9][1] = 0.005;
3832  fC[9][2] = 1.0000;
3833  fC[9][3] = 0.2812;
3834  fC[9][4] = 0.2699;
3835 
3836  fC[10][0] = 0.005;
3837  fC[10][1] = 0.005;
3838  fC[10][2] = 1.0000;
3839  fC[10][3] = 0.3328;
3840  fC[10][4] = 0.3714;
3841 
3842  fC[11][0] = 0.005;
3843  fC[11][1] = 0.005;
3844  fC[11][2] = 1.0000;
3845  fC[11][3] = 0.3780;
3846  fC[11][4] = 0.4810;
3847 
3848  fC[12][0] = 0.005;
3849  fC[12][1] = 0.005;
3850  fC[12][2] = 1.0000;
3851  fC[12][3] = 0.4125;
3852  fC[12][4] = 0.5771;
3853 
3854  fC[13][0] = 0.005;
3855  fC[13][1] = 0.005;
3856  fC[13][2] = 1.0000;
3857  fC[13][3] = 0.4486;
3858  fC[13][4] = 0.6799;
3859 
3860  fC[14][0] = 0.005;
3861  fC[14][1] = 0.005;
3862  fC[14][2] = 1.0000;
3863  fC[14][3] = 0.4840;
3864  fC[14][4] = 0.7668;
3865 
3866  fC[15][0] = 0.005;
3867  fC[15][1] = 0.005;
3868  fC[15][2] = 1.0000;
3869  fC[15][3] = 0.4971;
3870  fC[15][4] = 0.8288;
3871 
3872  fC[16][0] = 0.005;
3873  fC[16][1] = 0.005;
3874  fC[16][2] = 1.0000;
3875  fC[16][3] = 0.4956;
3876  fC[16][4] = 0.8653;
3877 
3878  fC[17][0] = 0.005;
3879  fC[17][1] = 0.005;
3880  fC[17][2] = 1.0000;
3881  fC[17][3] = 0.5173;
3882  fC[17][4] = 0.9059;
3883 
3884  for(Int_t i=18;i<fgkPIDptBin;i++){
3885  fC[i][0] = fC[17][0];
3886  fC[i][1] = fC[17][1];
3887  fC[i][2] = fC[17][2];
3888  fC[i][3] = fC[17][3];
3889  fC[i][4] = fC[17][4];
3890  }
3891  }
3892  // 10-20%
3893  else if(centrCur < 20){
3894  fC[0][0] = 0.005;
3895  fC[0][1] = 0.005;
3896  fC[0][2] = 1.0000;
3897  fC[0][3] = 0.0010;
3898  fC[0][4] = 0.0010;
3899 
3900  fC[1][0] = 0.005;
3901  fC[1][1] = 0.005;
3902  fC[1][2] = 1.0000;
3903  fC[1][3] = 0.0132;
3904  fC[1][4] = 0.0088;
3905 
3906  fC[2][0] = 0.005;
3907  fC[2][1] = 0.005;
3908  fC[2][2] = 1.0000;
3909  fC[2][3] = 0.0283;
3910  fC[2][4] = 0.0068;
3911 
3912  fC[3][0] = 0.005;
3913  fC[3][1] = 0.005;
3914  fC[3][2] = 1.0000;
3915  fC[3][3] = 0.0577;
3916  fC[3][4] = 0.0279;
3917 
3918  fC[4][0] = 0.005;
3919  fC[4][1] = 0.005;
3920  fC[4][2] = 1.0000;
3921  fC[4][3] = 0.0884;
3922  fC[4][4] = 0.0534;
3923 
3924  fC[5][0] = 0.005;
3925  fC[5][1] = 0.005;
3926  fC[5][2] = 1.0000;
3927  fC[5][3] = 0.1179;
3928  fC[5][4] = 0.0794;
3929 
3930  fC[6][0] = 0.005;
3931  fC[6][1] = 0.005;
3932  fC[6][2] = 1.0000;
3933  fC[6][3] = 0.1480;
3934  fC[6][4] = 0.1058;
3935 
3936  fC[7][0] = 0.005;
3937  fC[7][1] = 0.005;
3938  fC[7][2] = 1.0000;
3939  fC[7][3] = 0.1807;
3940  fC[7][4] = 0.1366;
3941 
3942  fC[8][0] = 0.005;
3943  fC[8][1] = 0.005;
3944  fC[8][2] = 1.0000;
3945  fC[8][3] = 0.2219;
3946  fC[8][4] = 0.1891;
3947 
3948  fC[9][0] = 0.005;
3949  fC[9][1] = 0.005;
3950  fC[9][2] = 1.0000;
3951  fC[9][3] = 0.2804;
3952  fC[9][4] = 0.2730;
3953 
3954  fC[10][0] = 0.005;
3955  fC[10][1] = 0.005;
3956  fC[10][2] = 1.0000;
3957  fC[10][3] = 0.3283;
3958  fC[10][4] = 0.3660;
3959 
3960  fC[11][0] = 0.005;
3961  fC[11][1] = 0.005;
3962  fC[11][2] = 1.0000;
3963  fC[11][3] = 0.3710;
3964  fC[11][4] = 0.4647;
3965 
3966  fC[12][0] = 0.005;
3967  fC[12][1] = 0.005;
3968  fC[12][2] = 1.0000;
3969  fC[12][3] = 0.4093;
3970  fC[12][4] = 0.5566;
3971 
3972  fC[13][0] = 0.005;
3973  fC[13][1] = 0.005;
3974  fC[13][2] = 1.0000;
3975  fC[13][3] = 0.4302;
3976  fC[13][4] = 0.6410;
3977 
3978  fC[14][0] = 0.005;
3979  fC[14][1] = 0.005;
3980  fC[14][2] = 1.0000;
3981  fC[14][3] = 0.4649;
3982  fC[14][4] = 0.7055;
3983 
3984  fC[15][0] = 0.005;
3985  fC[15][1] = 0.005;
3986  fC[15][2] = 1.0000;
3987  fC[15][3] = 0.4523;
3988  fC[15][4] = 0.7440;
3989 
3990  fC[16][0] = 0.005;
3991  fC[16][1] = 0.005;
3992  fC[16][2] = 1.0000;
3993  fC[16][3] = 0.4591;
3994  fC[16][4] = 0.7799;
3995 
3996  fC[17][0] = 0.005;
3997  fC[17][1] = 0.005;
3998  fC[17][2] = 1.0000;
3999  fC[17][3] = 0.4804;
4000  fC[17][4] = 0.8218;
4001 
4002  for(Int_t i=18;i<fgkPIDptBin;i++){
4003  fC[i][0] = fC[17][0];
4004  fC[i][1] = fC[17][1];
4005  fC[i][2] = fC[17][2];
4006  fC[i][3] = fC[17][3];
4007  fC[i][4] = fC[17][4];
4008  }
4009  }
4010  // 20-30%
4011  else if(centrCur < 30){
4012  fC[0][0] = 0.005;
4013  fC[0][1] = 0.005;
4014  fC[0][2] = 1.0000;
4015  fC[0][3] = 0.0010;
4016  fC[0][4] = 0.0010;
4017 
4018  fC[1][0] = 0.005;
4019  fC[1][1] = 0.005;
4020  fC[1][2] = 1.0000;
4021  fC[1][3] = 0.0102;
4022  fC[1][4] = 0.0064;
4023 
4024  fC[2][0] = 0.005;
4025  fC[2][1] = 0.005;
4026  fC[2][2] = 1.0000;
4027  fC[2][3] = 0.0292;
4028  fC[2][4] = 0.0066;
4029 
4030  fC[3][0] = 0.005;
4031  fC[3][1] = 0.005;
4032  fC[3][2] = 1.0000;
4033  fC[3][3] = 0.0597;
4034  fC[3][4] = 0.0296;
4035 
4036  fC[4][0] = 0.005;
4037  fC[4][1] = 0.005;
4038  fC[4][2] = 1.0000;
4039  fC[4][3] = 0.0900;
4040  fC[4][4] = 0.0589;
4041 
4042  fC[5][0] = 0.005;
4043  fC[5][1] = 0.005;
4044  fC[5][2] = 1.0000;
4045  fC[5][3] = 0.1199;
4046  fC[5][4] = 0.0859;
4047 
4048  fC[6][0] = 0.005;
4049  fC[6][1] = 0.005;
4050  fC[6][2] = 1.0000;
4051  fC[6][3] = 0.1505;
4052  fC[6][4] = 0.1141;
4053 
4054  fC[7][0] = 0.005;
4055  fC[7][1] = 0.005;
4056  fC[7][2] = 1.0000;
4057  fC[7][3] = 0.1805;
4058  fC[7][4] = 0.1454;
4059 
4060  fC[8][0] = 0.005;
4061  fC[8][1] = 0.005;
4062  fC[8][2] = 1.0000;
4063  fC[8][3] = 0.2221;
4064  fC[8][4] = 0.2004;
4065 
4066  fC[9][0] = 0.005;
4067  fC[9][1] = 0.005;
4068  fC[9][2] = 1.0000;
4069  fC[9][3] = 0.2796;
4070  fC[9][4] = 0.2838;
4071 
4072  fC[10][0] = 0.005;
4073  fC[10][1] = 0.005;
4074  fC[10][2] = 1.0000;
4075  fC[10][3] = 0.3271;
4076  fC[10][4] = 0.3682;
4077 
4078  fC[11][0] = 0.005;
4079  fC[11][1] = 0.005;
4080  fC[11][2] = 1.0000;
4081  fC[11][3] = 0.3648;
4082  fC[11][4] = 0.4509;
4083 
4084  fC[12][0] = 0.005;
4085  fC[12][1] = 0.005;
4086  fC[12][2] = 1.0000;
4087  fC[12][3] = 0.3988;
4088  fC[12][4] = 0.5339;
4089 
4090  fC[13][0] = 0.005;
4091  fC[13][1] = 0.005;
4092  fC[13][2] = 1.0000;
4093  fC[13][3] = 0.4315;
4094  fC[13][4] = 0.5995;
4095 
4096  fC[14][0] = 0.005;
4097  fC[14][1] = 0.005;
4098  fC[14][2] = 1.0000;
4099  fC[14][3] = 0.4548;
4100  fC[14][4] = 0.6612;
4101 
4102  fC[15][0] = 0.005;
4103  fC[15][1] = 0.005;
4104  fC[15][2] = 1.0000;
4105  fC[15][3] = 0.4744;
4106  fC[15][4] = 0.7060;
4107 
4108  fC[16][0] = 0.005;
4109  fC[16][1] = 0.005;
4110  fC[16][2] = 1.0000;
4111  fC[16][3] = 0.4899;
4112  fC[16][4] = 0.7388;
4113 
4114  fC[17][0] = 0.005;
4115  fC[17][1] = 0.005;
4116  fC[17][2] = 1.0000;
4117  fC[17][3] = 0.4411;
4118  fC[17][4] = 0.7293;
4119 
4120  for(Int_t i=18;i<fgkPIDptBin;i++){
4121  fC[i][0] = fC[17][0];
4122  fC[i][1] = fC[17][1];
4123  fC[i][2] = fC[17][2];
4124  fC[i][3] = fC[17][3];
4125  fC[i][4] = fC[17][4];
4126  }
4127  }
4128  // 30-40%
4129  else if(centrCur < 40){
4130  fC[0][0] = 0.005;
4131  fC[0][1] = 0.005;
4132  fC[0][2] = 1.0000;
4133  fC[0][3] = 0.0010;
4134  fC[0][4] = 0.0010;
4135 
4136  fC[1][0] = 0.005;
4137  fC[1][1] = 0.005;
4138  fC[1][2] = 1.0000;
4139  fC[1][3] = 0.0102;
4140  fC[1][4] = 0.0048;
4141 
4142  fC[2][0] = 0.005;
4143  fC[2][1] = 0.005;
4144  fC[2][2] = 1.0000;
4145  fC[2][3] = 0.0306;
4146  fC[2][4] = 0.0079;
4147 
4148  fC[3][0] = 0.005;
4149  fC[3][1] = 0.005;
4150  fC[3][2] = 1.0000;
4151  fC[3][3] = 0.0617;
4152  fC[3][4] = 0.0338;
4153 
4154  fC[4][0] = 0.005;
4155  fC[4][1] = 0.005;
4156  fC[4][2] = 1.0000;
4157  fC[4][3] = 0.0920;
4158  fC[4][4] = 0.0652;
4159 
4160  fC[5][0] = 0.005;
4161  fC[5][1] = 0.005;
4162  fC[5][2] = 1.0000;
4163  fC[5][3] = 0.1211;
4164  fC[5][4] = 0.0955;
4165 
4166  fC[6][0] = 0.005;
4167  fC[6][1] = 0.005;
4168  fC[6][2] = 1.0000;
4169  fC[6][3] = 0.1496;
4170  fC[6][4] = 0.1242;
4171 
4172  fC[7][0] = 0.005;
4173  fC[7][1] = 0.005;
4174  fC[7][2] = 1.0000;
4175  fC[7][3] = 0.1807;
4176  fC[7][4] = 0.1576;
4177 
4178  fC[8][0] = 0.005;
4179  fC[8][1] = 0.005;
4180  fC[8][2] = 1.0000;
4181  fC[8][3] = 0.2195;
4182  fC[8][4] = 0.2097;
4183 
4184  fC[9][0] = 0.005;
4185  fC[9][1] = 0.005;
4186  fC[9][2] = 1.0000;
4187  fC[9][3] = 0.2732;
4188  fC[9][4] = 0.2884;
4189 
4190  fC[10][0] = 0.005;
4191  fC[10][1] = 0.005;
4192  fC[10][2] = 1.0000;
4193  fC[10][3] = 0.3204;
4194  fC[10][4] = 0.3679;<