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