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