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