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