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