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