AliPhysics  a6017e1 (a6017e1)
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 
1367  double fEstimatorsCorrelationCoef[2];
1368  double fEstimatorsSigmaPars[4];
1369  double fDeltaEstimatorNsigma[2];
1370 
1371  fEstimatorsCorrelationCoef[0] = 0.0157497;
1372  fEstimatorsCorrelationCoef[1] = 0.973488;
1373  fEstimatorsSigmaPars[0] = 0.673612;
1374  fEstimatorsSigmaPars[1] = 0.0290718;
1375  fEstimatorsSigmaPars[2] = -0.000546728;
1376  fEstimatorsSigmaPars[3] = 5.82749e-06;
1377  fDeltaEstimatorNsigma[0] = 5.;
1378  fDeltaEstimatorNsigma[1] = 5.5;
1379 
1380  const double center = centrCL1 * fEstimatorsCorrelationCoef[1] + fEstimatorsCorrelationCoef[0];
1381  const double sigma = fEstimatorsSigmaPars[0] + fEstimatorsSigmaPars[1] * centrCL1 + fEstimatorsSigmaPars[2] * centrCL1 * centrCL1 + fEstimatorsSigmaPars[3] * centrCL1 * centrCL1 * centrCL1;
1382  if (centrV0M < center - fDeltaEstimatorNsigma[0] * sigma && centrV0M > center + fDeltaEstimatorNsigma[1] * sigma) {return kFALSE;}
1383  //if(fabs(centrV0M-centrCL1)>7.5) return kFALSE;
1384  fhEventCounter->Fill("Centr. Est. Consis. OK",1);
1385 
1386 
1387  // primary vertex selection
1388  const AliAODVertex* vtx = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertex());
1389  if(!vtx || vtx->GetNContributors() < 1){ return kFALSE;}
1390  fhEventCounter->Fill("PV OK",1);
1391 
1392  // SPD vertex selection
1393  const AliAODVertex* vtxSPD = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertexSPD());
1394 
1395 
1396  if (((AliAODHeader*)fEventAOD->GetHeader())->GetRefMultiplicityComb08() < 0) return kFALSE;
1397 
1398  if (fEventAOD->IsIncompleteDAQ()) return kFALSE;
1399 
1400  if (vtx->GetNContributors() < 2 || vtxSPD->GetNContributors()<1) return kFALSE;
1401  double cov[6]={0}; double covSPD[6]={0};
1402  vtx->GetCovarianceMatrix(cov);
1403  vtxSPD->GetCovarianceMatrix(covSPD);
1404 
1405  //Special selections for SPD vertex
1406  double zRes = TMath::Sqrt(covSPD[5]);
1407  double fMaxResol=0.25;
1408  if ( vtxSPD->IsFromVertexerZ() && (zRes>fMaxResol)) return kFALSE;
1409  fhEventCounter->Fill("SPD Vtx OK",1);
1410 
1411 
1412  // check for multi-vertexer pile-up
1413  /*
1414  const int kMinPileUpContrib = 5;
1415  const double kMaxPileUpChi2 = 5.0;
1416  const double kMinWDist = 15;
1417  const AliAODVertex* vtxPileUp = 0;
1418  int nPileUp = 0;
1419 
1420  nPileUp=fEventAOD->GetNumberOfPileupVerticesTracks();
1421 
1422  if (nPileUp) {
1423  if (vtx == vtxSPD) return kTRUE; // there are pile-up vertices but no primary
1424  Int_t bcPrim = vtPrm->GetBC();
1425  for (int iPileUp=0;iPileUp<nPileUp;iPileUp++) {
1426  vtxPileUp = (const AliAODVertex*)fEventAOD->GetPileupVertexTracks(iPileUp);
1427  if (vtxPileUp->GetNContributors() < kMinPileUpContrib) continue;
1428  if (vtxPileUp->GetChi2perNDF() > kMaxPileUpChi2) continue;
1429  int bcPlp = vtxPileUp->GetBC(); ///newly added
1430  if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other Bunch crossing (BC)
1431  double wDst = GetWDist(vtx,vtxPileUp);
1432  if (wDst<kMinWDist) continue;
1433  return kTRUE; // pile-up: well separated vertices
1434  }
1435  }
1436 */
1437 // fhEventCounter->Fill("Pileup MV OK",1);
1439  AliAnalysisUtils utils;
1440  utils.SetMinPlpContribMV(5);
1441  utils.SetMaxPlpChi2MV(5);
1442  utils.SetMinWDistMV(15);
1443  utils.SetCheckPlpFromDifferentBCMV(kTRUE);
1444 
1445  Bool_t isPileupFromMV = utils.IsPileUpMV(fEventAOD);
1446 
1447  if(isPileupFromMV) return kFALSE;
1448  fhEventCounter->Fill("Pileup MV OK",1);
1449 
1450  //Bool_t fRejectOutOfBunchPileUp = kFALSE;
1451  //if(fRejectOutOfBunchPileUp) // out-of-bunch rejection (provided by Christian)
1452  //{
1453  //out-of-bunch
1454  // if (utils.IsOutOfBunchPileUp(fEventAOD)){ return kFALSE; }
1455  //}
1456  fhEventCounter->Fill("Out-of-bunch Pileup OK",1);
1458  // check vertex consistency
1459  double dz = vtx->GetZ() - vtxSPD->GetZ();
1460  double errTot = TMath::Sqrt(cov[5]+covSPD[5]);
1461  double err = TMath::Sqrt(cov[5]);
1462  double nsigTot = dz/errTot;
1463  double nsig = dz/err;
1464  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsig)>20) return kFALSE;
1465  fhEventCounter->Fill("Vtx Consis. OK",1);
1466 
1467 
1468  const Double_t aodVtxZ = vtx->GetZ();
1469  if( TMath::Abs(aodVtxZ) > fPVtxCutZ ){return kFALSE;}
1470  fhEventCounter->Fill("PV #it{z} OK",1);
1471 
1472  // cut on # ESD tracks vs # TPC only tracks
1473  const Int_t nTracks = fEventAOD->GetNumberOfTracks();
1474  Int_t multEsd = ((AliAODHeader*)fEventAOD->GetHeader())->GetNumberOfESDTracks();
1475  Int_t multTrk = 0;
1476  //Int_t multTrkBefC = 0;
1477  Int_t multTrkTOF = 0;
1478  Int_t multTPC = 0;
1479  for (Int_t it = 0; it < nTracks; it++) {
1480  AliAODTrack* AODTrk = (AliAODTrack*)fEventAOD->GetTrack(it);
1481  if (!AODTrk){ delete AODTrk; continue; }
1482  if (AODTrk->TestFilterBit(128)) {multTPC++;}
1483  if (AODTrk->TestFilterBit(32)){
1484  multTrk++;
1485  if ( TMath::Abs(AODTrk->GetTOFsignalDz()) <= 10 && AODTrk->GetTOFsignal() >= 12000 && AODTrk->GetTOFsignal() <= 25000) multTrkTOF++;
1486  }
1487  } // end of for (Int_t it = 0; it < nTracks; it++)
1488  Double_t multTPCn = multTPC;
1489  Double_t multEsdn = multEsd;
1490  Double_t fESDvsTPConlyLinearCut[2];
1491 
1492  fESDvsTPConlyLinearCut[0] = 700.;
1493  fESDvsTPConlyLinearCut[1] = 3.38;
1494 
1495  if(multEsdn > fESDvsTPConlyLinearCut[0] + fESDvsTPConlyLinearCut[1] * multTPCn) return kFALSE;
1496  fhEventCounter->Fill("ESD TPC Mult. Diff. OK",1);
1497 
1498  Int_t fTOFvsFB32nSigmaCut[2];
1499  fTOFvsFB32nSigmaCut[0] = 4.;
1500  fTOFvsFB32nSigmaCut[1] = 4.;
1501 
1502  Double_t multTrkn = multTrk;
1503  Double_t multTrkTOFn = multTrkTOF;
1504 
1505  Double_t fTOFvsFB32correlationPars[4];
1506  Double_t fTOFvsFB32sigmaPars[6];
1507 
1508  fTOFvsFB32correlationPars[0] = -1.0178;
1509  fTOFvsFB32correlationPars[1] = 0.333132;
1510  fTOFvsFB32correlationPars[2] = 9.10282e-05;
1511  fTOFvsFB32correlationPars[3] = -1.61861e-08;
1512 
1513  fTOFvsFB32sigmaPars[0] = 1.47848;
1514  fTOFvsFB32sigmaPars[1] = 0.0385923;
1515  fTOFvsFB32sigmaPars[2] = -5.06153e-05;
1516  fTOFvsFB32sigmaPars[3] = 4.37641e-08;
1517  fTOFvsFB32sigmaPars[4] = -1.69082e-11;
1518  fTOFvsFB32sigmaPars[5] = 2.35085e-15;
1519 
1520  //Double_t mu32tof = PolN(multTrkn,fTOFvsFB32correlationPars,3);
1521  //Double_t sigma32tof = PolN(multTrkn,fTOFvsFB32sigmaPars, 5);
1522 
1523  Double_t mu32tof = fTOFvsFB32correlationPars[0] + fTOFvsFB32correlationPars[1]* multTrkn + fTOFvsFB32correlationPars[2]* pow(multTrkn,2) + fTOFvsFB32correlationPars[3]* pow(multTrkn,3);
1524  Double_t sigma32tof = fTOFvsFB32sigmaPars[0] + fTOFvsFB32sigmaPars[1]* multTrkn + fTOFvsFB32sigmaPars[2]* pow(multTrkn,2) + fTOFvsFB32sigmaPars[3]* pow(multTrkn,3) + fTOFvsFB32sigmaPars[4]* pow(multTrkn,4) + fTOFvsFB32sigmaPars[5]* pow(multTrkn,5);
1525 
1526 
1527  if ((multTrkTOFn > mu32tof + fTOFvsFB32nSigmaCut[0] * sigma32tof || multTrkTOFn < mu32tof - fTOFvsFB32nSigmaCut[1] * sigma32tof)) return kFALSE;
1528  //if(fExtraPileUp && multTrkTOFn< (-32+ 0.32*multTrkn+0.000037*multTrkn*multTrkn)) return kFALSE;
1529  //if(fExtraPileUp && multTrkTOFn> (13+0.46*multTrkn+0.000018*multTrkn*multTrkn)) return kFALSE;
1530  fhEventCounter->Fill("TOF fb32 Mult. correlation OK",1);
1531 
1532  fhEventCounter->Fill("Selected",1);
1533 
1534  // Fill event QA AFTER cuts
1535  if(fFillQA) FillEventsQA(1);
1536 
1537  return kTRUE;
1538 
1539 }
1540 //_____________________________________________________________________________
1542 {
1543  // Event selection for small system collision recorder in Run 2 year 2016
1544  // pp (LHC16kl), pPb (LHC16rqts)
1545  // return kTRUE if event passes all criteria, kFALSE otherwise
1546  // *************************************************************
1547 
1548  fhEventCounter->Fill("Input",1);
1549 
1550  // Physics selection (trigger)
1551  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1552  AliInputEventHandler* inputHandler = (AliInputEventHandler*) mgr->GetInputEventHandler();
1553  UInt_t fSelectMask = inputHandler->IsEventSelected();
1554 
1555  Bool_t isTriggerSelected = kFALSE;
1556  switch(fTrigger) // check for high multiplicity trigger
1557  {
1558  case 0:
1559  isTriggerSelected = fSelectMask& AliVEvent::kINT7;
1560  break;
1561 
1562  case 1:
1563  isTriggerSelected = fSelectMask& AliVEvent::kHighMultV0;
1564  break;
1565 
1566  case 2:
1567  isTriggerSelected = fSelectMask& AliVEvent::kHighMultSPD;
1568  break;
1569 
1570  default: isTriggerSelected = kFALSE;
1571  }
1572 
1573  if(!isTriggerSelected)
1574  return kFALSE;
1575 
1576  // events passing physics selection
1577  fhEventCounter->Fill("Physics selection OK",1);
1578 
1579  // primary vertex selection
1580  const AliAODVertex* vtx = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertex());
1581  if(!vtx || vtx->GetNContributors() < 1)
1582  return kFALSE;
1583  fhEventCounter->Fill("PV OK",1);
1584 
1585  // SPD vertex selection
1586  const AliAODVertex* vtxSPD = dynamic_cast<const AliAODVertex*>(fEventAOD->GetPrimaryVertexSPD());
1587 
1588  Double_t dMaxResol = 0.25; // suggested from DPG
1589  Double_t cov[6] = {0};
1590  vtxSPD->GetCovarianceMatrix(cov);
1591  Double_t zRes = TMath::Sqrt(cov[5]);
1592  if ( vtxSPD->IsFromVertexerZ() && (zRes > dMaxResol)) return kFALSE;
1593  fhEventCounter->Fill("SPD Vtx OK",1);
1594 
1595  // PileUp rejection included in Physics selection
1596  // but with values for high mult pp (> 5 contrib) => for low ones: do manually (> 3 contrib)
1597 
1598  /*
1599  if(fTrigger == 0 && fAOD->IsPileupFromSPD(3,0.8) )
1600  {
1601  return kFALSE;
1602  }
1603  */
1604 
1605  //fhEventCounter->Fill("Pileup SPD OK",1);
1606 
1607  // pileup rejection from multivertexer
1608 
1609  AliAnalysisUtils utils;
1610  utils.SetMinPlpContribMV(5);
1611  utils.SetMaxPlpChi2MV(5);
1612  utils.SetMinWDistMV(15);
1613  utils.SetCheckPlpFromDifferentBCMV(kFALSE);
1614  Bool_t isPileupFromMV = utils.IsPileUpMV(fEventAOD);
1615  utils.SetCheckPlpFromDifferentBCMV(kTRUE);
1616 
1617  if(isPileupFromMV) return kFALSE;
1618  fhEventCounter->Fill("Pileup MV OK",1);
1619 
1620  //Bool_t fRejectOutOfBunchPileUp = kTRUE; //you have to add this to the task and set from header
1621  //if(fRejectOutOfBunchPileUp) // out-of-bunch rejection (provided by Christian)
1622  //{
1623  //out-of-bunch
1624  // if (utils.IsOutOfBunchPileUp(fEventAOD)){ return kFALSE; }
1625  //}
1626  fhEventCounter->Fill("Out-of-bunch Pileup OK",1);
1627 
1628  // if (utils.IsSPDClusterVsTrackletBG(fEventAOD))
1629  // {
1630  // return kFALSE;
1631  // }
1632  //
1633  // fhEventCounter->Fill("SPDClTrBG OK",1);
1634  //
1635  // // SPD pileup
1636  // if (utils.IsPileUpSPD(fEventAOD))
1637  // {
1638  // return kFALSE;
1639  // }
1640  //
1641  // fhEventCounter->Fill("SPDPU OK",1);
1642  // }
1643 
1644  //fhEventCounter->Fill("Utils OK",1);
1645 
1646  // cutting on PV z-distance
1647  const Double_t aodVtxZ = vtx->GetZ();
1648  if( TMath::Abs(aodVtxZ) > fPVtxCutZ )
1649  {
1650  return kFALSE;
1651  }
1652  fhEventCounter->Fill("PV #it{z} OK",1);
1653 
1654  fhEventCounter->Fill("Selected",1);
1655  return kTRUE;
1656 }
1657 //_____________________________________________________________________________
1659 {
1660  // Filling various QA plots related with event selection
1661  // *************************************************************
1662 
1663  const AliAODVertex* aodVtx = fEventAOD->GetPrimaryVertex();
1664  const Double_t dVtxZ = aodVtx->GetZ();
1665  const Int_t iNumContr = aodVtx->GetNContributors();
1666  const AliAODVertex* spdVtx = fEventAOD->GetPrimaryVertexSPD();
1667  const Int_t iNumContrSPD = spdVtx->GetNContributors();
1668  const Double_t spdVtxZ = spdVtx->GetZ();
1669 
1670  fhQAEventsPVz[iQAindex]->Fill(dVtxZ);
1671  fhQAEventsNumContrPV[iQAindex]->Fill(iNumContr);
1672  fhQAEventsNumSPDContrPV[iQAindex]->Fill(iNumContrSPD);
1673  fhQAEventsDistPVSPD[iQAindex]->Fill(TMath::Abs(dVtxZ - spdVtxZ));
1674 
1675  // SPD vertexer resolution
1676  Double_t cov[6] = {0};
1677  spdVtx->GetCovarianceMatrix(cov);
1678  Double_t zRes = TMath::Sqrt(cov[5]);
1679  fhQAEventsSPDresol[iQAindex]->Fill(zRes);
1680 
1681  AliMultSelection* fMultSelection = 0x0;
1682  fMultSelection = (AliMultSelection*) fEventAOD->FindListObject("MultSelection");
1683  Double_t centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1684  Double_t centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
1685 
1686  fhQAEventsCentralityOutliers[iQAindex]->Fill(centrV0M,centrCL1);
1687 
1688  const Int_t nTracks = fEventAOD->GetNumberOfTracks();
1689  Int_t multEsd = ((AliAODHeader*)fEventAOD->GetHeader())->GetNumberOfESDTracks();
1690  Int_t multTPC = 0;
1691  Int_t multTrk = 0;
1692  Int_t multTrkTOF = 0;
1693  for (Int_t it = 0; it < nTracks; it++) {
1694  AliAODTrack* AODTrk = (AliAODTrack*)fEventAOD->GetTrack(it);
1695  if (!AODTrk){ delete AODTrk; continue; }
1696  if (AODTrk->TestFilterBit(128)) {multTPC++;}
1697  if (AODTrk->TestFilterBit(32)){
1698  multTrk++;
1699  if ( TMath::Abs(AODTrk->GetTOFsignalDz()) <= 10 && AODTrk->GetTOFsignal() >= 12000 && AODTrk->GetTOFsignal() <= 25000) multTrkTOF++;
1700  }
1701 
1702  }// end of for (Int_t it = 0; it < nTracks; it++)
1703  if(iQAindex==1){
1704  Double_t multTPCn = multTPC;
1705  Double_t multEsdn = multEsd;
1706 
1707  //Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
1708  //if (multESDTPCDif > 700.) return;//15000
1709 
1710  Double_t fESDvsTPConlyLinearCut[2];
1711 
1712  fESDvsTPConlyLinearCut[0] = 700.;
1713  fESDvsTPConlyLinearCut[1] = 3.38;
1714 
1715  if(multEsdn > fESDvsTPConlyLinearCut[0] + fESDvsTPConlyLinearCut[1] * multTPCn) return;
1716  }
1717  fhQAEventsPileUp[iQAindex]->Fill(multTPC,multEsd);
1718 
1719 
1720  if(iQAindex==1){
1721  Double_t multTrkn = multTrk;
1722  Double_t multTrkTOFn = multTrkTOF;
1723 
1724  Int_t fTOFvsFB32nSigmaCut[2];
1725  fTOFvsFB32nSigmaCut[0] = 4.;
1726  fTOFvsFB32nSigmaCut[1] = 4.;
1727 
1728  Double_t fTOFvsFB32correlationPars[4];
1729  Double_t fTOFvsFB32sigmaPars[6];
1730 
1731  fTOFvsFB32correlationPars[0] = -1.0178;
1732  fTOFvsFB32correlationPars[1] = 0.333132;
1733  fTOFvsFB32correlationPars[2] = 9.10282e-05;
1734  fTOFvsFB32correlationPars[3] = -1.61861e-08;
1735 
1736  fTOFvsFB32sigmaPars[0] = 1.47848;
1737  fTOFvsFB32sigmaPars[1] = 0.0385923;
1738  fTOFvsFB32sigmaPars[2] = -5.06153e-05;
1739  fTOFvsFB32sigmaPars[3] = 4.37641e-08;
1740  fTOFvsFB32sigmaPars[4] = -1.69082e-11;
1741  fTOFvsFB32sigmaPars[5] = 2.35085e-15;
1742 
1743  //Double_t mu32tof = PolN(multTrkn,fTOFvsFB32correlationPars,3);
1744  //Double_t sigma32tof = PolN(multTrkn,fTOFvsFB32sigmaPars, 5);
1745 
1746  Double_t mu32tof = fTOFvsFB32correlationPars[0] + fTOFvsFB32correlationPars[1]* multTrkn + fTOFvsFB32correlationPars[2]* pow(multTrkn,2) + fTOFvsFB32correlationPars[3]* pow(multTrkn,3);
1747  Double_t sigma32tof = fTOFvsFB32sigmaPars[0] + fTOFvsFB32sigmaPars[1]* multTrkn + fTOFvsFB32sigmaPars[2]* pow(multTrkn,2) + fTOFvsFB32sigmaPars[3]* pow(multTrkn,3) + fTOFvsFB32sigmaPars[4]* pow(multTrkn,4) + fTOFvsFB32sigmaPars[5]* pow(multTrkn,5);
1748 
1749  if ((multTrkTOFn > mu32tof + fTOFvsFB32nSigmaCut[0] * sigma32tof || multTrkTOFn < mu32tof - fTOFvsFB32nSigmaCut[1] * sigma32tof)) return;
1750 
1751  //if(fExtraPileUp && multTrkTOFn< (-32+ 0.32*multTrkn+0.000037*multTrkn*multTrkn)) return;
1752  //if(fExtraPileUp && multTrkTOFn> (13+0.46*multTrkn+0.000018*multTrkn*multTrkn)) return;
1753  }
1754  fhEventsMultTOFFilterbit32[iQAindex]->Fill(multTrk,multTrkTOF);
1755 
1756  return;
1757 }
1758 //_____________________________________________________________________________
1760 {
1761 
1762  // main (envelope) method for filtering all particles of interests (POIs) in selected events
1763  // All POIs passing selection criteria are saved to relevant TClonesArray for further processing
1764  // return kTRUE if succesfull (no errors in process)
1765  // *************************************************************
1766  if(!fProcessCharged && !fProcessPID) // if neither is ON, filtering is skipped
1767  return kFALSE; //return;
1768 
1769  fVectorCharged->clear();
1770  FilterCharged();
1771 
1772  // estimate centrality & assign indexes (centrality/percentile, ...)
1773  if(fColSystem == kPbPb){
1774 
1776  if(fIndexCentrality < 0) return kFALSE; // return; not succesfull estimation
1777  Double_t Mult = fVectorCharged->size();
1778  if(fExtraPileUp && fIndexCentrality< (-1.5*TMath::Power(Mult,0.46)-0.6*TMath::Log(Mult)*TMath::Log(Mult)+81)) return kFALSE;
1779  if(fExtraPileUp && fIndexCentrality>(-2.3*TMath::Power(Mult,0.39)-0.9*TMath::Log(Mult)*TMath::Log(Mult)+110)) return kFALSE;
1780 
1781  }
1782  if(fColSystem == kPP){fIndexCentrality = 1;}
1783 
1785 
1787 
1789 
1791 
1792  if(fProcessPID)
1793  {
1794 
1795  fVectorPion->clear();
1796  fVectorKaon->clear();
1797  fVectorProton->clear();
1798  FilterPID();
1799  }
1800 
1801  return kTRUE; //return;
1802 }
1803 //_____________________________________________________________________________
1805 {
1806  // Filtering input charged tracks
1807  // If track passes all requirements as defined in IsChargedSelected(),
1808  // the relevant properties (pT, eta, phi) are stored in FlowPart struct
1809  // and pushed to relevant vector container.
1810  // return kFALSE if any complications occurs
1811  // *************************************************************
1812 
1813  const Short_t iNumTracks = fEventAOD->GetNumberOfTracks();
1814  if(iNumTracks < 1) return;
1815 
1816  AliAODTrack* track = 0x0;
1817  Int_t iNumRefs = 0;
1818  Double_t NUAweight = 0;
1819  Double_t NUEweight = 0;
1820 
1821  Int_t c = -1;
1822  if(fIndexCentrality>=0 && fIndexCentrality<5) c = 0;
1823  if(fIndexCentrality>=5 && fIndexCentrality<10) c = 1;
1824  if(fIndexCentrality>=10 && fIndexCentrality<20) c = 2;
1825  if(fIndexCentrality>=20 && fIndexCentrality<30) c = 3;
1826  if(fIndexCentrality>=30 && fIndexCentrality<40) c = 4;
1827  if(fIndexCentrality>=40 && fIndexCentrality<50) c = 5;
1828 
1829  for(Short_t iTrack(0); iTrack < iNumTracks; iTrack++)
1830  {
1831  track = static_cast<AliAODTrack*>(fEventAOD->GetTrack(iTrack));
1832  if(!track) continue;
1833 
1834  if(fFillQA) FillQACharged(0,track); // QA before selection
1835 
1836  if(fNegativelyChargedRef==kTRUE && track->Charge()>0) continue;
1837  if(fPositivelyChargedRef==kTRUE && track->Charge()<0) continue;
1838 
1839  if(IsChargedSelected(track))
1840  {
1841  fVectorCharged->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kCharged) );
1843  fh3BeforeNUAWeightsCharged->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
1844  fhBeforeNUEWeightsCharged->Fill(track->Pt());
1845  }
1846  if(fFlowUseNUAWeights)
1847  {
1848  if(track->Charge()>0){ NUAweight = fh3NUAWeightChargedPlus->GetBinContent( fh3NUAWeightChargedPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1849  if(track->Charge()<0){ NUAweight = fh3NUAWeightChargedMinus->GetBinContent( fh3NUAWeightChargedMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1850 
1851  fh3AfterNUAWeightsCharged->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
1852  }
1853  if(fFlowUseNUEWeights){
1854 
1855  if(track->Charge()>0 && c>=0 && c<6){NUEweight = fhNUEWeightChargedPlus[c]->GetBinContent( fhNUEWeightChargedPlus[c]->FindBin(track->Pt()));}
1856  if(track->Charge()<0 && c>=0 && c<6){NUEweight = fhNUEWeightChargedMinus[c]->GetBinContent( fhNUEWeightChargedMinus[c]->FindBin(track->Pt()));}
1857 
1858  fhAfterNUEWeightsCharged->Fill(track->Pt(),NUEweight);
1859  }
1860 
1861  if(fFillQA) FillQACharged(1,track); // QA after selection
1862 
1863  // Filling refs QA plots
1864  if(fCutFlowRFPsPtMin > 0. && track->Pt() >= fCutFlowRFPsPtMin && fCutFlowRFPsPtMax > 0. && track->Pt() <= fCutFlowRFPsPtMax)
1865  {
1866  iNumRefs++;
1868  fh3BeforeNUAWeightsRefs->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
1869  fhBeforeNUEWeightsRefs->Fill(track->Pt());
1870  }
1871  if(fFlowUseNUAWeights)
1872  {
1873  if(track->Charge()>0){ NUAweight = fh3NUAWeightRefsPlus->GetBinContent( fh3NUAWeightRefsPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1874  if(track->Charge()<0){ NUAweight = fh3NUAWeightRefsMinus->GetBinContent( fh3NUAWeightRefsMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
1875 
1876  fh3AfterNUAWeightsRefs->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
1877  }
1878  if(fFlowUseNUEWeights)
1879  {
1880  if(track->Charge()>0 && c>=0 && c<6){NUEweight = fhNUEWeightRefsPlus[c]->GetBinContent( fhNUEWeightRefsPlus[c]->FindBin(track->Pt()) );}
1881  if(track->Charge()<0 && c>=0 && c<6){NUEweight = fhNUEWeightRefsMinus[c]->GetBinContent( fhNUEWeightRefsMinus[c]->FindBin(track->Pt()) );}
1882 
1883  fhAfterNUEWeightsRefs->Fill(track->Pt(),NUEweight);
1884  }
1885  FillQARefs(1,track);
1886  }
1887  }
1888  }
1889 
1890  // fill QA charged multiplicity
1891  fh2RefsMult->Fill(fIndexCentrality,iNumRefs);
1892  if(fFillQA)
1893  {
1894  fhQAChargedMult[0]->Fill(fEventAOD->GetNumberOfTracks());
1895  fhQAChargedMult[1]->Fill(fVectorCharged->size());
1896  }
1897 
1898  return;
1899 }
1900 //_____________________________________________________________________________
1902 {
1903  // Selection of charged track
1904  // returns kTRUE if track pass all requirements, kFALSE otherwise
1905  // *************************************************************
1906  if(!track) return kFALSE;
1907  fhChargedCounter->Fill("Input",1);
1908 
1909  // filter bit
1910  if( !track->TestFilterBit(fCutChargedTrackFilterBit) ) return kFALSE;
1911  fhChargedCounter->Fill("FB",1);
1912 
1913  // number of TPC clusters (additional check for not ITS-standalone tracks)
1914  if( track->GetTPCNcls() < fCutChargedNumTPCclsMin && fCutChargedTrackFilterBit != 2) return kFALSE;
1915  fhChargedCounter->Fill("#TPC-Cls",1);
1916 
1917  // track chi2 per space points
1918  Double_t chi2TPC =0.;
1919  chi2TPC = track->Chi2perNDF();
1920  if (fMaxChi2perTPCcls > 0. && chi2TPC > fMaxChi2perTPCcls) return kFALSE;
1921  fhChargedCounter->Fill("TPC-Chi2",1);
1922 
1923  // track DCA coordinates
1924  // note AliAODTrack::XYZAtDCA() works only for constrained tracks
1925  Double_t dTrackXYZ[3] = {0};
1926  Double_t dVertexXYZ[3] = {0.};
1927  Double_t dDCAXYZ[3] = {0.};
1928  if( fCutChargedDCAzMax > 0. || fCutChargedDCAxyMax > 0.)
1929  {
1930  const AliAODVertex* vertex = fEventAOD->GetPrimaryVertex();
1931  if(!vertex) return kFALSE; // event does not have a PV
1932 
1933  track->GetXYZ(dTrackXYZ);
1934  vertex->GetXYZ(dVertexXYZ);
1935 
1936  for(Short_t i(0); i < 3; i++)
1937  dDCAXYZ[i] = dTrackXYZ[i] - dVertexXYZ[i];
1938  }
1939 
1940  if(fCutChargedDCAzMax > 0. && TMath::Abs(dDCAXYZ[2]) > fCutChargedDCAzMax) return kFALSE;
1941  fhChargedCounter->Fill("DCA-z",1);
1942 
1943  if(fCutChargedDCAxyMax > 0. && TMath::Sqrt(dDCAXYZ[0]*dDCAXYZ[0] + dDCAXYZ[1]*dDCAXYZ[1]) > fCutChargedDCAxyMax) return kFALSE;
1944  fhChargedCounter->Fill("DCA-xy",1);
1945 
1946  // pseudorapidity (eta)
1947  if(fCutChargedEtaMax > 0. && TMath::Abs(track->Eta()) > fCutChargedEtaMax) return kFALSE;
1948  fhChargedCounter->Fill("Eta",1);
1949 
1950  // track passing all criteria
1951  fhChargedCounter->Fill("Selected",1);
1952  return kTRUE;
1953 }
1954 //_____________________________________________________________________________
1955 void AliAnalysisTaskFlowModes::FillQARefs(const Short_t iQAindex, const AliAODTrack* track)
1956 {
1957  // Filling various QA plots related to RFPs subset of charged track selection
1958  // *************************************************************
1959 
1960  if(!track) return;
1961  if(iQAindex == 0) return; // NOTE implemented only for selected RFPs
1962 
1963  fh2RefsPt->Fill(fIndexCentrality,track->Pt());
1964  fh2RefsEta->Fill(fIndexCentrality,track->Eta());
1965  fh2RefsPhi->Fill(fIndexCentrality,track->Phi());
1966 
1967  return;
1968 }
1969 //_____________________________________________________________________________
1970 void AliAnalysisTaskFlowModes::FillQACharged(const Short_t iQAindex, const AliAODTrack* track)
1971 {
1972  // Filling various QA plots related to charged track selection
1973  // *************************************************************
1974  if(!track) return;
1975 
1976  // filter bit testing
1977  for(Short_t i(0); i < 32; i++)
1978  {
1979  if(track->TestFilterBit(TMath::Power(2.,i)))
1980  fhQAChargedFilterBit[iQAindex]->Fill(i);
1981  }
1982 
1983  // track charge
1984  fhQAChargedCharge[iQAindex]->Fill(track->Charge());
1985 
1986  // number of TPC clusters
1987  fhQAChargedNumTPCcls[iQAindex]->Fill(track->GetTPCNcls());
1988 
1989  // track DCA
1990  Double_t dDCAXYZ[3] = {-999., -999., -999.};
1991  const AliAODVertex* vertex = fEventAOD->GetPrimaryVertex();
1992  if(vertex)
1993  {
1994  Double_t dTrackXYZ[3] = {-999., -999., -999.};
1995  Double_t dVertexXYZ[3] = {-999., -999., -999.};
1996 
1997  track->GetXYZ(dTrackXYZ);
1998  vertex->GetXYZ(dVertexXYZ);
1999 
2000  for(Short_t i(0); i < 3; i++)
2001  dDCAXYZ[i] = dTrackXYZ[i] - dVertexXYZ[i];
2002  }
2003  fhQAChargedDCAxy[iQAindex]->Fill(TMath::Sqrt(dDCAXYZ[0]*dDCAXYZ[0] + dDCAXYZ[1]*dDCAXYZ[1]));
2004  fhQAChargedDCAz[iQAindex]->Fill(dDCAXYZ[2]);
2005 
2006  // kinematics
2007  fhQAChargedPt[iQAindex]->Fill(track->Pt());
2008  fhQAChargedPhi[iQAindex]->Fill(track->Phi());
2009  fhQAChargedEta[iQAindex]->Fill(track->Eta());
2010 
2011  return;
2012 }
2013 //_____________________________________________________________________________
2015 {
2016  // If track passes all requirements as defined in IsPIDSelected() (and species dependent),
2017  // the relevant properties (pT, eta, phi) are stored in FlowPart struct
2018  // and pushed to relevant vector container.
2019  // return kFALSE if any complications occurs
2020  // *************************************************************
2021 
2022  const Short_t iNumTracks = fEventAOD->GetNumberOfTracks();
2023  if(iNumTracks < 1) return;
2024 
2025  PartSpecies species = kUnknown;
2026  AliAODTrack* track = 0x0;
2027  Double_t NUAweight = 0;
2028  Double_t NUEweight = 0;
2029 
2030 
2031  Int_t c = -1;
2032  if(fIndexCentrality>=0 && fIndexCentrality<5) c = 0;
2033  if(fIndexCentrality>=5 && fIndexCentrality<10) c = 1;
2034  if(fIndexCentrality>=10 && fIndexCentrality<20) c = 2;
2035  if(fIndexCentrality>=20 && fIndexCentrality<30) c = 3;
2036  if(fIndexCentrality>=30 && fIndexCentrality<40) c = 4;
2037  if(fIndexCentrality>=40 && fIndexCentrality<50) c = 5;
2038 
2039  for(Short_t iTrack(0); iTrack < iNumTracks; iTrack++)
2040  {
2041  track = static_cast<AliAODTrack*>(fEventAOD->GetTrack(iTrack));
2042  if(!track) continue;
2043 
2044  // PID tracks are subset of selected charged tracks (same quality requirements)
2045  if(!IsChargedSelected(track)) continue;
2046 
2047  if(fFillQA) FillPIDQA(0,track,kUnknown); // filling QA for tracks before selection (but after charged criteria applied)
2048 
2049  if(fNegativelyChargedPOI==kTRUE && track->Charge()>0) continue;
2050  if(fPositivelyChargedPOI==kTRUE && track->Charge()<0) continue;
2051 
2052  // PID track selection (return most favourable species)
2053  PartSpecies species = IsPIDSelected(track);
2054  // check if only protons should be used
2055  if(fCutPIDUseAntiProtonOnly && species == kProton && track->Charge() == 1) species = kUnknown;
2056 
2057  // selection of PID tracks
2058  switch (species)
2059  {
2060  case kPion:
2061  fVectorPion->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kPion, fPDGMassPion, track->Px(), track->Py(), track->Pz()) );
2063  fh3BeforeNUAWeightsPion->Fill(track->Phi(), track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
2064  fhBeforeNUEWeightsPion->Fill(track->Pt());
2065  }
2066  if(fFlowUseNUAWeights)
2067  {
2068  if(track->Charge() > 0){ NUAweight = fh3NUAWeightPionPlus->GetBinContent( fh3NUAWeightPionPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2069  if(track->Charge() < 0){ NUAweight = fh3NUAWeightPionMinus->GetBinContent( fh3NUAWeightPionMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2070 
2071  fh3AfterNUAWeightsPion->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
2072  }
2073  if(fFlowUseNUEWeights)
2074  {
2075 
2076  if(track->Charge() > 0 && c>=0 && c<6){ NUEweight = fhNUEWeightPionPlus[c]->GetBinContent( fhNUEWeightPionPlus[c]->FindBin(track->Pt()) );}
2077  if(track->Charge() < 0 && c>=0 && c<6){ NUEweight = fhNUEWeightPionMinus[c]->GetBinContent( fhNUEWeightPionMinus[c]->FindBin(track->Pt()) );}
2078 
2079  fhAfterNUEWeightsPion->Fill(track->Pt(),NUEweight);
2080  }
2081  break;
2082  case kKaon:
2083  fVectorKaon->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kKaon, fPDGMassKaon, track->Px(), track->Py(), track->Pz()) );
2085  fh3BeforeNUAWeightsKaon->Fill(track->Phi(), track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
2086  fhBeforeNUEWeightsKaon->Fill(track->Pt());
2087  }
2088  if(fFlowUseNUAWeights)
2089  {
2090  if(track->Charge() > 0){NUAweight = fh3NUAWeightKaonPlus->GetBinContent( fh3NUAWeightKaonPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2091  if(track->Charge() < 0){NUAweight = fh3NUAWeightKaonMinus->GetBinContent( fh3NUAWeightKaonMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2092 
2093  fh3AfterNUAWeightsKaon->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
2094  }
2095  if(fFlowUseNUEWeights)
2096  {
2097  if(track->Charge() > 0 && c>=0 && c<6){ NUEweight = fhNUEWeightKaonPlus[c]->GetBinContent( fhNUEWeightKaonPlus[c]->FindBin(track->Pt()) );}
2098  if(track->Charge() < 0 && c>=0 && c<6){ NUEweight = fhNUEWeightKaonMinus[c]->GetBinContent( fhNUEWeightKaonMinus[c]->FindBin(track->Pt()) );}
2099 
2100  fhAfterNUEWeightsKaon->Fill(track->Pt(),NUEweight);
2101  }
2102  break;
2103  case kProton:
2104  fVectorProton->emplace_back( FlowPart(track->Pt(),track->Phi(),track->Eta(), track->Charge(), kProton, fPDGMassProton, track->Px(), track->Py(), track->Pz()) );
2106  fh3BeforeNUAWeightsProton->Fill(track->Phi(), track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ());
2107  fhBeforeNUEWeightsProton->Fill(track->Pt());
2108  }
2109  if(fFlowUseNUAWeights)
2110  {
2111  if(track->Charge() > 0){NUAweight = fh3NUAWeightProtonPlus->GetBinContent( fh3NUAWeightProtonPlus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2112  if(track->Charge() < 0){NUAweight = fh3NUAWeightProtonMinus->GetBinContent( fh3NUAWeightProtonMinus->FindBin(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ()) );}
2113 
2114  fh3AfterNUAWeightsProton->Fill(track->Phi(),track->Eta(),fEventAOD->GetPrimaryVertex()->GetZ(),NUAweight);
2115  }
2116  if(fFlowUseNUEWeights)
2117  {
2118  if(track->Charge() > 0 && c>=0 && c<6){ NUEweight = fhNUEWeightProtonPlus[c]->GetBinContent( fhNUEWeightProtonPlus[c]->FindBin(track->Pt()) );}
2119  if(track->Charge() < 0 && c>=0 && c<6){ NUEweight = fhNUEWeightProtonMinus[c]->GetBinContent( fhNUEWeightProtonMinus[c]->FindBin(track->Pt()) );}
2120 
2121  fhAfterNUEWeightsProton->Fill(track->Pt(),NUEweight);
2122  }
2123  break;
2124  default:
2125  break;
2126  }
2127 
2128  //if(fFillQA) FillPIDQA(1,track,species); // filling QA for tracks AFTER selection
2129  }
2130 
2134 
2135  return;
2136 }
2137 //_____________________________________________________________________________
2139 {
2140  // Selection of PID tracks (pi,K,p) - track identification
2141  // nSigma cutting is used
2142  // returns AliAnalysisTaskFlowModes::PartSpecies enum : kPion, kKaon, kProton if any of this passed kUnknown otherwise
2143  // *************************************************************
2144 
2145  // checking detector states
2146  AliPIDResponse::EDetPidStatus pidStatusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, track);
2147  AliPIDResponse::EDetPidStatus pidStatusTOF = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
2148 
2149  Bool_t bIsTPCok = (pidStatusTPC == AliPIDResponse::kDetPidOk);
2150  Bool_t bIsTOFok = ((pidStatusTOF == AliPIDResponse::kDetPidOk) && (track->GetStatus()& AliVTrack::kTOFout) && (track->GetStatus()& AliVTrack::kTIME) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000)); // checking TOF
2151  //if(!bIsTPCok) return kUnknown;
2152 
2153  const Double_t dPt = track->Pt();
2154  const Double_t dP = track->P();
2155 
2156  // use nSigma cuts (based on combination of TPC / TOF nSigma cuts)
2157 
2158  Double_t dNumSigmaTPC[5] = {-99,-99,-99,-99,-99}; // TPC nSigma array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2159  Double_t dNumSigmaTOF[5] = {-99,-99,-99,-99,-99}; // TOF nSigma array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2160 
2161  Float_t *probabilities;
2162  Float_t mismProb;
2163  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
2164 
2165  Double_t probTPC[AliPID::kSPECIES]={0.};
2166  Double_t probTOF[AliPID::kSPECIES]={0.};
2167  Double_t probTPCTOF[AliPID::kSPECIES]={0.};
2168 
2169 
2170  // filling nSigma arrays
2171  if(bIsTPCok) // should be anyway
2172  {
2173  dNumSigmaTPC[0] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron));
2174  dNumSigmaTPC[1] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kMuon));
2175  dNumSigmaTPC[2] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion));
2176  dNumSigmaTPC[3] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
2177  dNumSigmaTPC[4] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton));
2178  }
2179 
2180  if(bIsTOFok) // should be anyway
2181  {
2182  dNumSigmaTOF[0] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron));
2183  dNumSigmaTOF[1] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kMuon));
2184  dNumSigmaTOF[2] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion));
2185  dNumSigmaTOF[3] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon));
2186  dNumSigmaTOF[4] = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton));
2187  }
2188 
2189 
2190  if(fPIDbayesian){
2191  if(!TPCTOFagree(track)){return kUnknown;}
2192 
2193  fBayesianResponse->SetDetResponse(fEventAOD, fCurrCentr,AliESDpid::kTOF_T0); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
2194  if(fEventAOD->GetTOFHeader()){
2195  fESDpid.SetTOFResponse(fEventAOD,AliESDpid::kTOF_T0);
2196  }
2198  }
2199 
2200  Int_t ParticleFlag[]={0,0,0,0};//Unknown (electron), pion, kaon, proton
2201  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
2202  UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
2203  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) { // TPC is available
2204 
2205  // TPC nSigma cuts
2206  if(dP>0.2 && dP <= 0.5)
2207  {
2208  if(fPIDnsigma){
2209  Double_t dMinSigmasTPC = TMath::MinElement(5,dNumSigmaTPC);
2210  // electron rejection
2211  if(dMinSigmasTPC == dNumSigmaTPC[0] && TMath::Abs(dNumSigmaTPC[0]) <= fCutPIDnSigmaTPCRejectElectron) ParticleFlag[0] = 1; //return kUnknown;
2212  if(dMinSigmasTPC == dNumSigmaTPC[2] && TMath::Abs(dNumSigmaTPC[2]) <= fCutPIDnSigmaPionMax) ParticleFlag[1] = 1;//return kPion;
2213  if(dMinSigmasTPC == dNumSigmaTPC[3] && TMath::Abs(dNumSigmaTPC[3]) <= fCutPIDnSigmaKaonMax) ParticleFlag[2] = 1;//return kKaon;
2214  if(dMinSigmasTPC == dNumSigmaTPC[4] && TMath::Abs(dNumSigmaTPC[4]) <= fCutPIDnSigmaProtonMax) ParticleFlag[3] = 1;//return kProton;
2215  }
2216  if(fPIDbayesian){
2217  ProbBayes[0] = probTPC[0];
2218  ProbBayes[1] = probTPC[1];
2219  ProbBayes[2] = probTPC[2];
2220  ProbBayes[3] = probTPC[3];
2221  ProbBayes[4] = probTPC[4];
2222 
2223  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2224  if(dMaxBayesianProb > fParticleProbability){
2225  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaTPC[0]) <= fCutPIDnSigmaTPCRejectElectron)ParticleFlag[0] = 1;// return kUnknown;
2226  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaTPC[2]) <= fCutPIDnSigmaPionMax)ParticleFlag[1] = 1;//return kPion;
2227  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaTPC[3]) <= fCutPIDnSigmaKaonMax)ParticleFlag[2] = 1;//return kKaon;
2228  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaTPC[4]) <= fCutPIDnSigmaProtonMax)ParticleFlag[3] = 1;//return kProton;
2229  }
2230  /*
2231  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
2232  probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
2233  ProbBayes[0] = probabilities[0];
2234  ProbBayes[1] = probabilities[1];
2235  ProbBayes[2] = probabilities[2];
2236  ProbBayes[3] = probabilities[3];
2237  ProbBayes[4] = probabilities[4];
2238 
2239  mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
2240 
2241  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2242  if(dMaxBayesianProb > fParticleProbability && mismProb < 0.5){
2243  // electron rejection
2244  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaTPC[0]) <= fCutPIDnSigmaTPCRejectElectron)ParticleFlag[0] = 1;// return kUnknown;
2245  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaTPC[2]) <= fCutPIDnSigmaPionMax)ParticleFlag[1] = 1;//return kPion;
2246  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaTPC[3]) <= fCutPIDnSigmaKaonMax)ParticleFlag[2] = 1;//return kKaon;
2247  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaTPC[4]) <= fCutPIDnSigmaProtonMax)ParticleFlag[3] = 1;//return kProton;
2248  }*/
2249  }
2250  }
2251  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
2252  detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
2253 
2254  if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()){
2255  AliAODPid* pidObj = track->GetDetPid();
2256  if ((detUsed >= AliPIDResponse::kDetTOF) && (pidObj && pidObj->GetTOFsignal() < 99999)){
2257 
2258  // combined TPC + TOF nSigma cuts
2259  if(dP > 0.5) // && < 4 GeV TODO once TPC dEdx parametrisation is available
2260  {
2261  Double_t dNumSigmaCombined[5] = {-99,-99,-99,-99,-99};
2262 
2263  // discard candidates if no TOF is available if cut is on
2264  if(fCutPIDnSigmaCombinedNoTOFrejection && !bIsTOFok) ParticleFlag[0] = 1;// return kUnknown;
2265 
2266  // calculating combined nSigmas
2267  for(Short_t i(0); i < 5; i++)
2268  {
2269  if(bIsTOFok) { dNumSigmaCombined[i] = TMath::Sqrt(dNumSigmaTPC[i]*dNumSigmaTPC[i] + dNumSigmaTOF[i]*dNumSigmaTOF[i]); }
2270  else { dNumSigmaCombined[i] = dNumSigmaTPC[i]; }
2271  }
2272 
2273  if(fPIDnsigma){
2274  Double_t dMinSigmasCombined = TMath::MinElement(5,dNumSigmaCombined);
2275 
2276  // electron rejection
2277  if(dMinSigmasCombined == dNumSigmaCombined[0] && TMath::Abs(dNumSigmaCombined[0]) <= fCutPIDnSigmaPionMax) ParticleFlag[0] = 1;// return kUnknown;
2278 
2279  switch (fPIDnsigmaCombination) {
2280  case 1:
2281  //combination 1
2282  if(dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 3.) ParticleFlag[1] = 1;//return kPion;
2283  if(dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.5) ParticleFlag[2] = 1; //return kKaon;
2284  if(dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 3.) ParticleFlag[3] = 1;//return kProton;
2285  break;
2286  case 2:
2287  //combination 2
2288  if(dP < 2.5 && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 3.) ParticleFlag[1] = 1; //return kPion;
2289  if(dP > 2.5 && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 2.) ParticleFlag[1] = 1; //return kPion;
2290 
2291  if(dP < 2. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.5) ParticleFlag[2] = 1; //return kKaon;
2292  if(dP > 2. && dP < 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.) ParticleFlag[2] = 1; //return kKaon;
2293  if(dP > 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 1.5) ParticleFlag[2] = 1; //return kKaon;
2294 
2295  if(dP < 3. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 3.) ParticleFlag[3] = 1; //return kProton;
2296  if(dP > 3. && dP < 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 2.) ParticleFlag[3] = 1; //return kProton;
2297  if(dP > 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 1.5) ParticleFlag[3] = 1; //return kProton;
2298 
2299  break;
2300  case 3:
2301  //combination 3
2302  if(dP < 2.5 && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 3.) ParticleFlag[1] = 1; //return kPion;
2303  if(dP > 2.5 && dP < 4. && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 1.5) ParticleFlag[1] = 1; //return kPion;
2304  if(dP > 4. && dMinSigmasCombined == dNumSigmaCombined[2] && TMath::Abs(dNumSigmaCombined[2]) <= 1.) ParticleFlag[1] = 1; //return kPion;
2305 
2306  if(dP < 2. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 2.5) ParticleFlag[2] = 1; //return kKaon;
2307  if(dP > 2. && dP < 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 1.5)ParticleFlag[2] = 1; //return kKaon;
2308  if(dP > 3. && dMinSigmasCombined == dNumSigmaCombined[3] && TMath::Abs(dNumSigmaCombined[3]) <= 1.) ParticleFlag[2] = 1; //return kKaon;
2309 
2310  if(dP < 3. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 3.) ParticleFlag[3] = 1; //return kProton;
2311  if(dP > 3. && dP < 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 2.) ParticleFlag[3] = 1; //return kProton;
2312  if(dP > 5. && dMinSigmasCombined == dNumSigmaCombined[4] && TMath::Abs(dNumSigmaCombined[4]) <= 1.) ParticleFlag[3] = 1; //return kProton;
2313  break;
2314  default:
2315  break;
2316  }
2317  }
2318  if(fPIDbayesian){
2319  ProbBayes[0] = probTPCTOF[0];
2320  ProbBayes[1] = probTPCTOF[1];
2321  ProbBayes[2] = probTPCTOF[2];
2322  ProbBayes[3] = probTPCTOF[3];
2323  ProbBayes[4] = probTPCTOF[4];
2324  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2325  if(dMaxBayesianProb > fParticleProbability){
2326  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaCombined[0]) <= fCutPIDnSigmaPionMax) ParticleFlag[0] = 1; //return kUnknown;
2327  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaCombined[2]) <= fCutPIDnSigmaPionMax) ParticleFlag[1] = 1; //return kPion;
2328  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaCombined[3]) <= fCutPIDnSigmaKaonMax) ParticleFlag[2] = 1; //return kKaon;
2329  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaCombined[4]) <= fCutPIDnSigmaProtonMax) ParticleFlag[3] = 1; //return kProton;
2330  }else{ParticleFlag[0] = 1;}
2331  /*
2332  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
2333  probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
2334  ProbBayes[0] = probabilities[0];
2335  ProbBayes[1] = probabilities[1];
2336  ProbBayes[2] = probabilities[2];
2337  ProbBayes[3] = probabilities[3];
2338  ProbBayes[4] = probabilities[4];
2339 
2340  mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
2341 
2342  Double_t dMaxBayesianProb = TMath::MaxElement(5,ProbBayes);
2343  if(dMaxBayesianProb > fParticleProbability && mismProb < 0.5){
2344  // electron rejection
2345  if(dMaxBayesianProb == ProbBayes[0] && TMath::Abs(dNumSigmaCombined[0]) <= fCutPIDnSigmaPionMax) ParticleFlag[0] = 1; //return kUnknown;
2346  if(dMaxBayesianProb == ProbBayes[2] && TMath::Abs(dNumSigmaCombined[2]) <= fCutPIDnSigmaPionMax) ParticleFlag[1] = 1; //return kPion;
2347  if(dMaxBayesianProb == ProbBayes[3] && TMath::Abs(dNumSigmaCombined[3]) <= fCutPIDnSigmaKaonMax) ParticleFlag[2] = 1; //return kKaon;
2348  if(dMaxBayesianProb == ProbBayes[4] && TMath::Abs(dNumSigmaCombined[4]) <= fCutPIDnSigmaProtonMax) ParticleFlag[3] = 1; //return kProton;
2349  }*/
2350  }
2351  }//if(dP > 0.5)
2352  }//if ((detUsed >= AliPIDResponse::kDetTOF) && (pidObj && pidObj->GetTOFsignal() < 99999))
2353  }//TPC or TOF if (detUsed)
2354  }//TPC if (detUsed)
2355 
2356  PartSpecies species;
2357  if(!ParticleFlag[0] && ParticleFlag[1] && !ParticleFlag[2] && !ParticleFlag[3]){
2358  species = kPion;
2359  }else if(!ParticleFlag[0] && !ParticleFlag[1] && ParticleFlag[2] && !ParticleFlag[3]){
2360  species = kKaon;
2361  }else if(!ParticleFlag[0] && !ParticleFlag[1] && !ParticleFlag[2] && ParticleFlag[3]){
2362  species = kProton;
2363  }else{species = kUnknown;}
2364 
2365  if(fFillQA) FillPIDQA(1,track,species); // filling QA for tracks AFTER selection
2366 
2367  return species;
2368 }
2369 //_____________________________________________________________________________
2370 void AliAnalysisTaskFlowModes::FillPIDQA(const Short_t iQAindex, const AliAODTrack* track, const PartSpecies species)
2371 {
2372  // Filling various QA plots related to PID (pi,K,p) track selection
2373  // *************************************************************
2374  if(!track) return;
2375  if(!fPIDResponse || !fPIDCombined)
2376  {
2377  ::Error("FillPIDQA","AliPIDResponse or AliPIDCombined object not found!");
2378  return;
2379  }
2380 
2381  // TPC & TOF statuses & measures
2382  AliPIDResponse::EDetPidStatus pidStatusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, track);
2383  AliPIDResponse::EDetPidStatus pidStatusTOF = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
2384 
2385  fhQAPIDTOFstatus[iQAindex]->Fill((Int_t) pidStatusTOF );
2386  fhQAPIDTPCstatus[iQAindex]->Fill((Int_t) pidStatusTPC );
2387 
2388  Bool_t bIsTPCok = (pidStatusTPC == AliPIDResponse::kDetPidOk);
2389  //Bool_t bIsTOFok = ((pidStatusTOF == AliPIDResponse::kDetPidOk) && (track->GetStatus()& AliVTrack::kTOFout) && (track->GetStatus()& AliVTrack::kTIME));
2390 
2391  Bool_t bIsTOFok = ((pidStatusTOF == AliPIDResponse::kDetPidOk) && (track->GetStatus()& AliVTrack::kTOFout) && (track->GetStatus()& AliVTrack::kTIME) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000)); // checking TOF
2392 
2393  Double_t dNumSigmaTPC[5] = {-11}; // array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2394  Double_t dNumSigmaTOF[5] = {-11}; // array: 0: electron / 1: muon / 2: pion / 3: kaon / 4: proton
2395 
2396  Double_t dTPCdEdx = -5; // TPC dEdx for selected particle
2397  Double_t dTOFbeta = -0.05; //TOF beta for selected particle
2398 
2399  Double_t dP = track->P();
2400  Double_t dPt = track->Pt();
2401 
2402  // detector status dependent
2403  if(bIsTPCok)
2404  {
2405  dNumSigmaTPC[0] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2406  dNumSigmaTPC[1] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kMuon);
2407  dNumSigmaTPC[2] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
2408  dNumSigmaTPC[3] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
2409  dNumSigmaTPC[4] = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2410 
2411  dTPCdEdx = track->GetTPCsignal();
2412  fhQAPIDTPCdEdx[iQAindex]->Fill(track->P(), dTPCdEdx);
2413  }
2414  else // TPC status not OK
2415  {
2416  dNumSigmaTPC[0] = -11.;
2417  dNumSigmaTPC[1] = -11.;
2418  dNumSigmaTPC[2] = -11.;
2419  dNumSigmaTPC[3] = -11.;
2420  dNumSigmaTPC[4] = -11.;
2421 
2422  fhQAPIDTPCdEdx[iQAindex]->Fill(track->P(), -5.);
2423  }
2424 
2425  if(bIsTOFok)
2426  {
2427  dNumSigmaTOF[0] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
2428  dNumSigmaTOF[1] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kMuon);
2429  dNumSigmaTOF[2] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
2430  dNumSigmaTOF[3] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
2431  dNumSigmaTOF[4] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
2432 
2433  Double_t dTOF[5];
2434  track->GetIntegratedTimes(dTOF);
2435  dTOFbeta = dTOF[0] / track->GetTOFsignal();
2436  fhQAPIDTOFbeta[iQAindex]->Fill(dP,dTOFbeta);
2437  }
2438  else // TOF status not OK
2439  {
2440  dNumSigmaTOF[0] = -11.;
2441  dNumSigmaTOF[1] = -11.;
2442  dNumSigmaTOF[2] = -11.;
2443  dNumSigmaTOF[3] = -11.;
2444  dNumSigmaTOF[4] = -11.;
2445 
2446  fhQAPIDTOFbeta[iQAindex]->Fill(track->P(),-0.05);
2447  }
2448 
2449 
2450  // species dependent QA
2451  switch (species)
2452  {
2453  case kPion:
2454  fh2PIDPionPt->Fill(fIndexCentrality,track->Pt());
2455  fh2PIDPionPhi->Fill(fIndexCentrality,track->Phi());
2456  fh2PIDPionEta->Fill(fIndexCentrality,track->Eta());
2457  fhPIDPionCharge->Fill(track->Charge());
2458  fh2PIDPionTPCdEdx->Fill(dPt,dTPCdEdx);
2459  fh2PIDPionTOFbeta->Fill(dPt,dTOFbeta);
2460  fh2PIDPionTPCnSigmaPion->Fill(dPt,dNumSigmaTPC[2]);
2461  fh2PIDPionTOFnSigmaPion->Fill(dPt,dNumSigmaTOF[2]);
2462  fh2PIDPionTPCnSigmaKaon->Fill(dPt,dNumSigmaTPC[3]);
2463  fh2PIDPionTOFnSigmaKaon->Fill(dPt,dNumSigmaTOF[3]);
2464  fh2PIDPionTPCnSigmaProton->Fill(dPt,dNumSigmaTPC[4]);
2465  fh2PIDPionTOFnSigmaProton->Fill(dPt,dNumSigmaTOF[4]);
2466 
2467  fh3PIDPionTPCTOFnSigmaPion[iQAindex]->Fill(dNumSigmaTPC[2],dNumSigmaTOF[2],dPt);
2468  fh3PIDPionTPCTOFnSigmaKaon[iQAindex]->Fill(dNumSigmaTPC[3],dNumSigmaTOF[3],dPt);
2469  fh3PIDPionTPCTOFnSigmaProton[iQAindex]->Fill(dNumSigmaTPC[4],dNumSigmaTOF[4],dPt);
2470 
2471  break;
2472 
2473  case kKaon:
2474  fh2PIDKaonPt->Fill(fIndexCentrality,track->Pt());
2475  fh2PIDKaonPhi->Fill(fIndexCentrality,track->Phi());
2476  fh2PIDKaonEta->Fill(fIndexCentrality,track->Eta());
2477  fhPIDKaonCharge->Fill(track->Charge());
2478  fh2PIDKaonTPCdEdx->Fill(dP,dTPCdEdx);
2479  fh2PIDKaonTOFbeta->Fill(dP,dTOFbeta);
2480  fh2PIDKaonTPCnSigmaPion->Fill(dPt,dNumSigmaTPC[2]);
2481  fh2PIDKaonTOFnSigmaPion->Fill(dPt,dNumSigmaTOF[2]);
2482  fh2PIDKaonTPCnSigmaKaon->Fill(dPt,dNumSigmaTPC[3]);
2483  fh2PIDKaonTOFnSigmaKaon->Fill(dPt,dNumSigmaTOF[3]);
2484  fh2PIDKaonTPCnSigmaProton->Fill(dPt,dNumSigmaTPC[4]);
2485  fh2PIDKaonTOFnSigmaProton->Fill(dPt,dNumSigmaTOF[4]);
2486 
2487  fh3PIDKaonTPCTOFnSigmaPion[iQAindex]->Fill(dNumSigmaTPC[2],dNumSigmaTOF[2],dPt);
2488  fh3PIDKaonTPCTOFnSigmaKaon[iQAindex]->Fill(dNumSigmaTPC[3],dNumSigmaTOF[3],dPt);
2489  fh3PIDKaonTPCTOFnSigmaProton[iQAindex]->Fill(dNumSigmaTPC[4],dNumSigmaTOF[4],dPt);
2490 
2491  break;
2492 
2493  case kProton:
2494  fh2PIDProtonPt->Fill(fIndexCentrality,track->Pt());
2495  fh2PIDProtonPhi->Fill(fIndexCentrality,track->Phi());
2496  fh2PIDProtonEta->Fill(fIndexCentrality,track->Eta());
2497  fhPIDProtonCharge->Fill(track->Charge());
2498  fh2PIDProtonTPCdEdx->Fill(dP,dTPCdEdx);
2499  fh2PIDProtonTOFbeta->Fill(dP,dTOFbeta);
2500  fh2PIDProtonTPCnSigmaPion->Fill(dPt,dNumSigmaTPC[2]);
2501  fh2PIDProtonTOFnSigmaPion->Fill(dPt,dNumSigmaTOF[2]);
2502  fh2PIDProtonTPCnSigmaKaon->Fill(dPt,dNumSigmaTPC[3]);
2503  fh2PIDProtonTOFnSigmaKaon->Fill(dPt,dNumSigmaTOF[3]);
2504  fh2PIDProtonTPCnSigmaProton->Fill(dPt,dNumSigmaTPC[4]);
2505  fh2PIDProtonTOFnSigmaProton->Fill(dPt,dNumSigmaTOF[4]);
2506 
2507  fh3PIDProtonTPCTOFnSigmaPion[iQAindex]->Fill(dNumSigmaTPC[2],dNumSigmaTOF[2],dPt);
2508  fh3PIDProtonTPCTOFnSigmaKaon[iQAindex]->Fill(dNumSigmaTPC[3],dNumSigmaTOF[3],dPt);
2509  fh3PIDProtonTPCTOFnSigmaProton[iQAindex]->Fill(dNumSigmaTPC[4],dNumSigmaTOF[4],dPt);
2510 
2511  break;
2512 
2513  default:
2514  break;
2515  }
2516 
2517 
2518  //
2519 
2520 
2521  return;
2522 }
2523 //_____________________________________________________________________________
2525 {
2526 
2527  // main method for processing of (selected) event:
2528  // - Filtering of tracks / particles for flow calculations
2529  // - Phi,eta,pt weights for generic framework are calculated if specified
2530  // - Flow calculations
2531  // returns kTRUE if succesfull
2532  // *************************************************************
2533 
2534  // printf("======= EVENT ================\n");
2535 
2536 
2537  // checking the run number for applying weights & loading TList with weights
2538  if(fRunNumber < 0 || fRunNumber != fEventAOD->GetRunNumber() )
2539  {
2540  fRunNumber = fEventAOD->GetRunNumber();
2541 
2543  {
2544  TDirectory* dirFlowNUAWeights = (TDirectory*) fFlowNUAWeightsFile->Get(Form("000%d",fRunNumber));
2545  if(!dirFlowNUAWeights) {::Error("ProcessEvent","TList from flow weights not found."); return kFALSE; }
2546  fh3NUAWeightRefsPlus = (TH3D*) dirFlowNUAWeights->Get("ChargedPlus"); if(!fh3NUAWeightRefsPlus) { ::Error("ProcessEvent","Positive Refs NUA weights not found"); return kFALSE; }
2547  fh3NUAWeightRefsMinus = (TH3D*) dirFlowNUAWeights->Get("ChargedMinus"); if(!fh3NUAWeightRefsMinus) { ::Error("ProcessEvent","Negative Refs NUA weights not found"); return kFALSE; }
2548 
2549  fh3NUAWeightChargedPlus = (TH3D*) dirFlowNUAWeights->Get("ChargedPlus"); if(!fh3NUAWeightChargedPlus) { ::Error("ProcessEvent","Positive Charged NUA weights not found"); return kFALSE; }
2550  fh3NUAWeightChargedMinus = (TH3D*) dirFlowNUAWeights->Get("ChargedMinus"); if(!fh3NUAWeightChargedMinus) { ::Error("ProcessEvent","Nagative Charged NUA weights not found"); return kFALSE; }
2551 
2552  fh3NUAWeightPionPlus = (TH3D*) dirFlowNUAWeights->Get("PionPlus"); if(!fh3NUAWeightPionPlus) { ::Error("ProcessEvent","Positive Pion NUA weights not found"); return kFALSE; }
2553  fh3NUAWeightPionMinus = (TH3D*) dirFlowNUAWeights->Get("PionMinus"); if(!fh3NUAWeightPionMinus) { ::Error("ProcessEvent","Negative Pion NUA weights not found"); return kFALSE; }
2554 
2555 
2556  fh3NUAWeightKaonPlus = (TH3D*) dirFlowNUAWeights->Get("KaonPlus"); if(!fh3NUAWeightKaonPlus) { ::Error("ProcessEvent","Positive Kaon NUA weights not found"); return kFALSE; }
2557  fh3NUAWeightKaonMinus = (TH3D*) dirFlowNUAWeights->Get("KaonMinus"); if(!fh3NUAWeightKaonMinus) { ::Error("ProcessEvent","Negative Kaon NUA weights not found"); return kFALSE; }
2558 
2559 
2560  fh3NUAWeightProtonPlus = (TH3D*) dirFlowNUAWeights->Get("ProtonPlus"); if(!fh3NUAWeightProtonPlus) { ::Error("ProcessEvent","Positive Proton NUA weights not found"); return kFALSE; }
2561  fh3NUAWeightProtonMinus = (TH3D*) dirFlowNUAWeights->Get("ProtonMinus"); if(!fh3NUAWeightProtonMinus) { ::Error("ProcessEvent","Negative Proton NUA weights not found"); return kFALSE; }
2562 
2563 
2564  }
2565  }
2566 
2567  if(fColSystem == kPbPb){
2569  if(fIndexCentrality < 0) { return kFALSE;}
2570  }
2571 
2573  {
2574  TDirectory* dirFlowNUEWeights[6] = {0x0};
2575 
2576  const char* gCentrality[] = {"0-5cc","5-10cc","10-20cc","20-30cc","30-40cc","40-50cc"};
2577  for(int iCentrality = 0; iCentrality<6;iCentrality++){
2578  dirFlowNUEWeights[iCentrality] = (TDirectory*) fFlowNUEWeightsFile->Get(gCentrality[iCentrality]);
2579  if(!dirFlowNUEWeights[iCentrality]) {::Error("ProcessEvent","TDirectoy from NUE weights not found."); return kFALSE; }
2580  fhNUEWeightRefsPlus[iCentrality] = dynamic_cast<TH1F*>( dirFlowNUEWeights[iCentrality]->FindObject("ChargedPlus")); if(!fhNUEWeightRefsPlus[iCentrality]) { ::Error("ProcessEvent","Positive Refs NUE weights not found"); return kFALSE; }
2581  fhNUEWeightRefsMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ChargedMinus")); if(!fhNUEWeightRefsMinus[iCentrality]) { ::Error("ProcessEvent","Negative Refs NUE weights not found"); return kFALSE; }
2582  fhNUEWeightChargedPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ChargedPlus")); if(!fhNUEWeightChargedPlus[iCentrality]) { ::Error("ProcessEvent","Positive Charged NUE weights not found"); return kFALSE; }
2583  fhNUEWeightChargedMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ChargedMinus")); if(!fhNUEWeightChargedMinus[iCentrality]) { ::Error("ProcessEvent","Negative Charged NUE weights not found"); return kFALSE; }
2584  fhNUEWeightPionPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("PionPlus")); if(!fhNUEWeightPionPlus[iCentrality]) { ::Error("ProcessEvent","Positive Pion NUE weights not found"); return kFALSE; }
2585  fhNUEWeightPionMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("PionMinus")); if(!fhNUEWeightPionMinus[iCentrality]) { ::Error("ProcessEvent","Negative Pion NUE weights not found"); return kFALSE; }
2586  fhNUEWeightKaonPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("KaonPlus")); if(!fhNUEWeightKaonPlus[iCentrality]) { ::Error("ProcessEvent","Positive Kaon NUE weights not found"); return kFALSE; }
2587  fhNUEWeightKaonMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("KaonMinus")); if(!fhNUEWeightKaonMinus[iCentrality]) { ::Error("ProcessEvent","Negative Kaon NUE weights not found"); return kFALSE; }
2588  fhNUEWeightProtonPlus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ProtonPlus")); if(!fhNUEWeightProtonPlus[iCentrality]) { ::Error("ProcessEvent","Positive Proton NUE weights not found"); return kFALSE; }
2589  fhNUEWeightProtonMinus[iCentrality] = dynamic_cast<TH1F*>(dirFlowNUEWeights[iCentrality]->FindObject("ProtonMinus")); if(!fhNUEWeightProtonMinus[iCentrality]) { ::Error("ProcessEvent","Negative Proton NUE weights not found"); return kFALSE; }
2590 
2591  }
2592  }
2593 
2594  // filtering particles
2595  if(!Filtering()) return kFALSE;
2596  // at this point, centrality index (percentile) should be properly estimated, if not, skip event
2597 
2598  // if running in kFillWeights mode, skip the remaining part
2599  if(fRunMode == kFillWeights) { fEventCounter++; return kTRUE; }
2600 
2601  // checking if there is at least one charged track selected;
2602  // if not, event is skipped: unable to compute Reference flow (and thus any differential flow)
2603  if(fVectorCharged->size() < 1)
2604  return kFALSE;
2605 
2606  // at this point, all particles fullfiling relevant POIs (and REFs) criteria are filled in TClonesArrays
2607 
2608  // >>>> flow starts here <<<<
2609  // >>>> Flow a la General Framework <<<<
2610  for(Short_t iGap(0); iGap < fNumEtaGap; iGap++)
2611  {
2612 
2613  // Reference (pT integrated) flow
2614  DoFlowRefs(iGap);
2615 
2616  // pT differential
2617  if(fProcessCharged)
2618  {
2619  // charged track flow
2620  DoFlowCharged(iGap);
2621  }
2622 
2623  if(fProcessPID)
2624  {
2625  const Int_t iSizePion = fVectorPion->size();
2626  const Int_t iSizeKaon = fVectorKaon->size();
2627  const Int_t iSizeProton = fVectorProton->size();
2628 
2629  // pi,K,p flow
2630  if(iSizePion > 0) DoFlowPID(iGap,kPion);
2631  if(iSizeKaon > 0) DoFlowPID(iGap,kKaon);
2632  if(iSizeProton > 0) DoFlowPID(iGap,kProton);
2633  }
2634  } // endfor {iGap} eta gaps
2635 
2636  fEventCounter++; // counter of processed events
2637  //printf("event %d\n",fEventCounter);
2638 
2639  return kTRUE;
2640 }
2641 //_____________________________________________________________________________
2643 {
2644 
2645  // Estimate <2>, <4> and <6> for reference flow for all harmonics based on relevant flow vectors
2646  // *************************************************************
2647 
2648  //Float_t dEtaGap = fEtaGap[iEtaGapIndex];
2649  Short_t iHarmonics = 0;
2650  Double_t Cn2 = 0;
2651  Double_t Dn4GapP = 0;
2652  Double_t Dn6GapP = 0;
2653  TComplex vector = TComplex(0,0,kFALSE);
2654  Double_t dValue = 999;
2655 
2656  FillRefsVectors(iEtaGapIndex); // filling RFPs (Q) flow vectors
2658  // estimating <2>
2659  Cn2 = TwoGap(0,0).Re();
2660  if(Cn2 != 0)
2661  {
2662  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2663  {
2664  iHarmonics = fHarmonics[iHarm];
2665  vector = TwoGap(iHarmonics,-iHarmonics);
2666  dValue = vector.Re()/Cn2;
2667  // 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);
2668  if( TMath::Abs(dValue < 1) )
2669  fpRefsCor2[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexCentrality, dValue, Cn2);
2670 
2671  }
2672  }
2673  }
2674  //for mixed harmonics
2676  //estimating <4>
2677  Dn4GapP = FourGapPos(0,0,0,0).Re();
2678  if(Dn4GapP != 0)
2679  {
2680  // (2,2 | 2,2)_gap , referece flow for v4/psi2
2681  TComplex Four_2222_GapP = FourGapPos(2, 2, -2, -2);
2682  double c4_2222_GapP = Four_2222_GapP.Re()/Dn4GapP;
2683  fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][0]->Fill(fIndexCentrality, c4_2222_GapP, Dn4GapP);
2684  // (3,3 | 3,3)_gap, referece flow for v6/psi3
2685  TComplex Four_3333_GapP = FourGapPos(3, 3, -3, -3);
2686  double c4_3333_GapP = Four_3333_GapP.Re()/Dn4GapP;
2687  fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][1]->Fill(fIndexCentrality, c4_3333_GapP, Dn4GapP);
2688  // (3,2 | 3,2)_gap, reference flow for v5/psi23
2689  TComplex Four_3232_GapP = FourGapPos(3, 2, -3, -2);
2690  double c4_3232_GapP = Four_3232_GapP.Re()/Dn4GapP;
2691  fpMixedRefsCor4[fIndexSampling][iEtaGapIndex][2]->Fill(fIndexCentrality, c4_3232_GapP, Dn4GapP);
2692  }
2693  //estimating <6>
2694  Dn6GapP = SixGapPos(0,0,0,0,0,0).Re();
2695  if(Dn6GapP != 0)
2696  {
2697  // (2,2,2 | 2,2,2)_gap, reference flow for v6/psi2
2698  TComplex Six_222222_GapP = SixGapPos(2, 2, 2, -2, -2, -2);
2699  double c6_222222_GapP = Six_222222_GapP.Re()/Dn6GapP;
2700  fpMixedRefsCor6[fIndexSampling][iEtaGapIndex]->Fill(fIndexCentrality, c6_222222_GapP, Dn6GapP);
2701  }
2702  }
2703  return;
2704 }
2705 //_____________________________________________________________________________
2707 {
2708  // Estimate <2'>, <3'> and <4'> for pT diff flow of charged tracks for all harmonics based on relevant flow vectors
2709  // *************************************************************
2710 
2711  FillPOIsVectors(iEtaGapIndex,kCharged); // filling POIs (P,S) flow vectors
2712 
2713  const Double_t dPtBinWidth = (fFlowPOIsPtMax - fFlowPOIsPtMin) / fFlowPOIsPtNumBins;
2714 
2715  //Float_t dEtaGap = fEtaGap[iEtaGapIndex];
2716  Short_t iHarmonics = 0;
2717  Double_t Dn2 = 0;
2718  Double_t DDn3GapP=0;
2719  Double_t DDn4GapP=0;
2720  TComplex vector = TComplex(0,0,kFALSE);
2721  Double_t dValue = 999;
2722 
2723 
2724  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
2725  {
2727  // estimating <2'>
2728  // POIs in positive eta
2729  Dn2 = TwoDiffGapPos(0,0,iPt).Re();
2730  if(Dn2 != 0)
2731  {
2732  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2733  {
2734  iHarmonics = fHarmonics[iHarm];
2735  vector = TwoDiffGapPos(iHarmonics,-iHarmonics,iPt);
2736  dValue = vector.Re()/Dn2;
2737  if( TMath::Abs(dValue < 1) )
2738  fp2ChargedCor2Pos[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
2739  }
2740  }
2741 
2742  // POIs in negative eta
2743  Dn2 = TwoDiffGapNeg(0,0,iPt).Re();
2744  if(Dn2 != 0)
2745  {
2746  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2747  {
2748  iHarmonics = fHarmonics[iHarm];
2749  vector = TwoDiffGapNeg(iHarmonics,-iHarmonics,iPt);
2750  dValue = vector.Re()/Dn2;
2751  if( TMath::Abs(dValue < 1) )
2752  fp2ChargedCor2Neg[fIndexSampling][iEtaGapIndex][iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
2753  }
2754  }
2755  }
2757  // estimating <3'>
2758  // POIs in positive eta
2759  DDn3GapP = ThreeDiffGapPos(0,0,0,iPt).Re();
2760  if(DDn3GapP!=0)
2761  {
2762  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
2763  {
2764  if(iMixedHarm==0){ vector = ThreeDiffGapPos(4,-2,-2,iPt);}
2765  if(iMixedHarm==1){ vector = ThreeDiffGapPos(6,-3,-3,iPt);}
2766  if(iMixedHarm==2){ vector = ThreeDiffGapPos(5,-3,-2,iPt);}
2767  dValue = vector.Re()/DDn3GapP;
2768  if( TMath::Abs(dValue < 1) ){
2769  fpMixedChargedCor3Pos[fIndexSampling][iEtaGapIndex][iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
2770  }
2771  }
2772  }
2773  // POIs in negative eta
2774  DDn3GapP = ThreeDiffGapNeg(0,0,0,iPt).Re();
2775  if(DDn3GapP!=0)
2776  {
2777  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
2778  {
2779  if(iMixedHarm==0){ vector = ThreeDiffGapNeg(4,-2,-2,iPt);}
2780  if(iMixedHarm==1){ vector = ThreeDiffGapNeg(6,-3,-3,iPt);}
2781  if(iMixedHarm==2){ vector = ThreeDiffGapNeg(5,-2,-3,iPt);}
2782  dValue = vector.Re()/DDn3GapP;
2783  if( TMath::Abs(dValue < 1) ){
2784  fpMixedChargedCor3Neg[fIndexSampling][iEtaGapIndex][iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
2785  }
2786  }
2787  }
2788  // estimating <4'>
2789  // POIs in positive eta
2790  DDn4GapP = Four13DiffGapPos(0,0,0,0,iPt).Re();
2791  if(DDn4GapP!=0)
2792  {
2793  vector = Four13DiffGapPos(6,-2,-2,-2,iPt);
2794  dValue = vector.Re()/DDn4GapP;
2795  if( TMath::Abs(dValue < 1) ){
2796  fpMixedChargedCor4Pos[fIndexSampling][iEtaGapIndex]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
2797  }
2798  }
2799  // POIs in negative eta
2800  DDn4GapP = Four13DiffGapNeg(0,0,0,0,iPt).Re();
2801  if(DDn4GapP!=0)
2802  {
2803  vector = Four13DiffGapNeg(6,-2,-2,-2,iPt);
2804  dValue = vector.Re()/DDn4GapP;
2805  if( TMath::Abs(dValue < 1) ){
2806  fpMixedChargedCor4Neg[fIndexSampling][iEtaGapIndex]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
2807  }
2808  }
2809  }
2810  } // endfor {iPt}
2811  return;
2812 }
2813 //_____________________________________________________________________________
2814 void AliAnalysisTaskFlowModes::DoFlowPID(const Short_t iEtaGapIndex, const PartSpecies species)
2815 {
2816  // Estimate <2'>, <3'> and <4'> for pT diff flow of pi/K/p tracks for all harmonics based on relevant flow vectors
2817  // *************************************************************
2818 
2819  TProfile2D** profile2Pos = 0x0;
2820  TProfile2D** profile2Neg = 0x0;
2821  //TProfile2D** profile4 = 0x0;
2822  TProfile2D** profile3Pos = 0x0;
2823  TProfile2D** profile3Neg = 0x0;
2824  TProfile2D* profile4Pos = 0x0;
2825  TProfile2D* profile4Neg = 0x0;
2826 
2827 
2829 
2830  switch (species)
2831  {
2832  case kPion:
2833  profile2Pos = fp2PionCor2Pos[fIndexSampling][iEtaGapIndex];
2834  profile2Neg = fp2PionCor2Neg[fIndexSampling][iEtaGapIndex];
2835  //profile4 = fp2PionCor4;
2836  break;
2837 
2838  case kKaon:
2839  profile2Pos = fp2KaonCor2Pos[fIndexSampling][iEtaGapIndex];
2840  profile2Neg = fp2KaonCor2Neg[fIndexSampling][iEtaGapIndex];
2841  //profile4 = fp2KaonCor4;
2842  break;
2843 
2844  case kProton:
2845  profile2Pos = fp2ProtonCor2Pos[fIndexSampling][iEtaGapIndex];
2846  profile2Neg = fp2ProtonCor2Neg[fIndexSampling][iEtaGapIndex];
2847  //profile4 = fp2ProtonCor4;
2848  break;
2849 
2850  default:
2851  ::Error("DoFlowPID","Unexpected species! Terminating!");
2852  return;
2853  }
2854  }
2855 
2857 
2858  switch (species)
2859  {
2860  case kPion:
2861  profile3Pos = fpMixedPionCor3Pos[fIndexSampling][iEtaGapIndex];
2862  profile3Neg = fpMixedPionCor3Neg[fIndexSampling][iEtaGapIndex];
2863 
2864  profile4Pos = fpMixedPionCor4Pos[fIndexSampling][iEtaGapIndex];
2865  profile4Neg = fpMixedPionCor4Neg[fIndexSampling][iEtaGapIndex];
2866  break;
2867 
2868  case kKaon:
2869  profile3Pos = fpMixedKaonCor3Pos[fIndexSampling][iEtaGapIndex];
2870  profile3Neg = fpMixedKaonCor3Neg[fIndexSampling][iEtaGapIndex];
2871 
2872  profile4Pos = fpMixedKaonCor4Pos[fIndexSampling][iEtaGapIndex];
2873  profile4Neg = fpMixedKaonCor4Neg[fIndexSampling][iEtaGapIndex];
2874  break;
2875 
2876  case kProton:
2877  profile3Pos = fpMixedProtonCor3Pos[fIndexSampling][iEtaGapIndex];
2878  profile3Neg = fpMixedProtonCor3Neg[fIndexSampling][iEtaGapIndex];
2879 
2880  profile4Pos = fpMixedProtonCor4Pos[fIndexSampling][iEtaGapIndex];
2881  profile4Neg = fpMixedProtonCor4Neg[fIndexSampling][iEtaGapIndex];
2882  break;
2883 
2884  default:
2885  ::Error("DoFlowPID with mixed harmonics","Unexpected species! Terminating!");
2886  return;
2887  }
2888  }
2889 
2890  FillPOIsVectors(iEtaGapIndex,species); // Filling POIs vectors
2891 
2892  const Double_t dPtBinWidth = (fFlowPOIsPtMax - fFlowPOIsPtMin) / fFlowPOIsPtNumBins;
2893 
2894  //Float_t dEtaGap = fEtaGap[iEtaGapIndex];
2895  Short_t iHarmonics = 0;
2896  Double_t Dn2 = 0;
2897  TComplex vector = TComplex(0,0,kFALSE);
2898  Double_t dValue = 999;
2899  Double_t DDn3GapP=0;
2900  Double_t DDn4GapP=0;
2901 
2902 
2903  // filling POIs (P,S) flow vectors
2904 
2905  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
2906  {
2908  // estimating <2'>
2909  // POIs in positive eta
2910  Dn2 = TwoDiffGapPos(0,0,iPt).Re();
2911  if(Dn2 != 0)
2912  {
2913  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2914  {
2915  iHarmonics = fHarmonics[iHarm];
2916  vector = TwoDiffGapPos(iHarmonics,-iHarmonics,iPt);
2917  dValue = vector.Re()/Dn2;
2918  if( TMath::Abs(dValue < 1) )
2919  profile2Pos[iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
2920  }
2921  }
2922  // POIs in negative eta
2923  Dn2 = TwoDiffGapNeg(0,0,iPt).Re();
2924  if(Dn2 != 0)
2925  {
2926  for(Short_t iHarm(0); iHarm < fNumHarmonics; iHarm++)
2927  {
2928  iHarmonics = fHarmonics[iHarm];
2929  vector = TwoDiffGapNeg(iHarmonics,-iHarmonics,iPt);
2930  dValue = vector.Re()/Dn2;
2931  if( TMath::Abs(dValue < 1) )
2932  profile2Neg[iHarm]->Fill(fIndexCentrality, iPt*dPtBinWidth, dValue, Dn2);
2933  }
2934  }
2935  }
2937  // estimating <3'>
2938  // POIs in positive eta
2939  DDn3GapP = ThreeDiffGapPos(0,0,0,iPt).Re();
2940  if(DDn3GapP!=0)
2941  {
2942  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
2943  {
2944  if(iMixedHarm==0){ vector = ThreeDiffGapPos(4,-2,-2,iPt);}
2945  if(iMixedHarm==1){ vector = ThreeDiffGapPos(6,-3,-3,iPt);}
2946  if(iMixedHarm==2){ vector = ThreeDiffGapPos(5,-3,-2,iPt);}
2947  dValue = vector.Re()/DDn3GapP;
2948  if( TMath::Abs(dValue < 1) ){
2949  profile3Pos[iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
2950  }
2951  }
2952  }
2953  // POIs in negative eta
2954  DDn3GapP = ThreeDiffGapNeg(0,0,0,iPt).Re();
2955  if(DDn3GapP!=0)
2956  {
2957  for(Short_t iMixedHarm(0); iMixedHarm < fNumMixedHarmonics-1; iMixedHarm++)
2958  {
2959  if(iMixedHarm==0){ vector = ThreeDiffGapNeg(4,-2,-2,iPt);}
2960  if(iMixedHarm==1){ vector = ThreeDiffGapNeg(6,-3,-3,iPt);}
2961  if(iMixedHarm==2){ vector = ThreeDiffGapNeg(5,-2,-3,iPt);}
2962  dValue = vector.Re()/DDn3GapP;
2963  if( TMath::Abs(dValue < 1) ){
2964  profile3Neg[iMixedHarm]->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn3GapP);
2965  }
2966  }
2967  }
2968 
2969  //estimating <4'>
2970  //POIs in positive eta
2971  DDn4GapP = Four13DiffGapPos(0,0,0,0,iPt).Re();
2972  if(DDn4GapP!=0)
2973  {
2974  vector = Four13DiffGapPos(6,-2,-2,-2,iPt);
2975  dValue = vector.Re()/DDn4GapP;
2976  if( TMath::Abs(dValue < 1) ){
2977  profile4Pos->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
2978  }
2979  }
2980  //POIs in negative eta
2981  DDn4GapP = Four13DiffGapNeg(0,0,0,0,iPt).Re();
2982  if(DDn4GapP!=0)
2983  {
2984  vector = Four13DiffGapNeg(6,-2,-2,-2,iPt);
2985  dValue = vector.Re()/DDn4GapP;
2986  if( TMath::Abs(dValue < 1) ){
2987  profile4Neg->Fill(fIndexCentrality,iPt*dPtBinWidth, dValue, DDn4GapP);
2988  }
2989  }
2990  }
2991  } // endfor {iPt}
2992  return;
2993 }
2994 //_____________________________________________________________________________
2996 {
2997  // Filling Q flow vector with RFPs
2998  // return kTRUE if succesfull (i.e. no error occurs), kFALSE otherwise
2999  // *************************************************************
3000  const Float_t dEtaGap = fEtaGap[iEtaGapIndex];
3001  TH3D* h3NUAWeights = 0x0;
3002  TH1F* hNUEWeights = 0x0;
3003  Double_t dNUAWeight = 1.;
3004  Double_t dNUEWeight = 1.;
3005 
3006  Int_t c = -1;
3007  if(fIndexCentrality>=0 && fIndexCentrality<5) c = 0;
3008  if(fIndexCentrality>=5 && fIndexCentrality<10) c = 1;
3009  if(fIndexCentrality>=10 && fIndexCentrality<20) c = 2;
3010  if(fIndexCentrality>=20 && fIndexCentrality<30) c = 3;
3011  if(fIndexCentrality>=30 && fIndexCentrality<40) c = 4;
3012  if(fIndexCentrality>=40 && fIndexCentrality<50) c = 5;
3013  // clearing output (global) flow vectors
3016 
3017  Double_t dQcosPos, dQcosNeg, dQsinPos, dQsinNeg;
3018 
3019  // Double_t dQcosPos[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax] = {0};
3020  // Double_t dQcosNeg[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax] = {0};
3021  // Double_t dQsinPos[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax] = {0};
3022  // Double_t dQsinNeg[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax] = {0};
3023 
3024  for (auto part = fVectorCharged->begin(); part != fVectorCharged->end(); part++)
3025  {
3026  // checking species of used particles (just for double checking purpose)
3027  if( part->species != kCharged)
3028  {
3029  ::Warning("FillRefsVectors","Unexpected part. species (%d) in selected sample (expected %d)",part->species,kCharged);
3030  continue;
3031  }
3032  if(fFlowUseNUAWeights)
3033  {
3034  if(part->charge >0) h3NUAWeights = fh3NUAWeightRefsPlus;
3035  if(part->charge <0) h3NUAWeights = fh3NUAWeightRefsMinus;
3036 
3037  if(!h3NUAWeights) { ::Error("FillRefsVectors","Histogram with NUA weights not found."); continue; }
3038  }
3039 
3040  if(fFlowUseNUEWeights)
3041  {
3042  if(part->charge >0 && c>=0 && c<6) hNUEWeights = fhNUEWeightRefsPlus[c];
3043  if(part->charge <0 && c>=0 && c<6) hNUEWeights = fhNUEWeightRefsMinus[c];
3044 
3045  if(!hNUEWeights) { ::Error("FillRefsVectors","Histogram with NUE weights not found."); return; }
3046  }
3047  // RFPs pT check
3048  if(fCutFlowRFPsPtMin > 0. && part->pt < fCutFlowRFPsPtMin)
3049  continue;
3050 
3051  if(fCutFlowRFPsPtMax > 0. && part->pt > fCutFlowRFPsPtMax)
3052  continue;
3053 
3054  // 0-ing variables
3055  dQcosPos = 0;
3056  dQcosNeg = 0;
3057  dQsinPos = 0;
3058  dQsinNeg = 0;
3059 
3060  // loading weights if needed
3061  if(fFlowUseNUAWeights && h3NUAWeights)
3062  {
3063  dNUAWeight = h3NUAWeights->GetBinContent(h3NUAWeights->FindBin(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ()));
3064  if(dNUAWeight <= 0) dNUAWeight = 1.;
3065  if(iEtaGapIndex == 0) fh3AfterNUAWeightsRefs->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3066  }
3067  if(fFlowUseNUEWeights && hNUEWeights)
3068  {
3069  dNUEWeight = hNUEWeights->GetBinContent(hNUEWeights->FindBin(part->pt));
3070  if(dNUEWeight <= 0) dNUEWeight = 1.;
3071  if(iEtaGapIndex == 0) fhAfterNUEWeightsRefs->Fill(part->pt, dNUEWeight);
3072  }
3073 
3074  // RPF candidate passing all criteria: start filling flow vectors
3075 
3076  if(part->eta > dEtaGap / 2)
3077  {
3078  // RFP in positive eta acceptance
3079  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3080  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3081  {
3082  dQcosPos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3083  dQsinPos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3084  fFlowVecQpos[iHarm][iPower] += TComplex(dQcosPos,dQsinPos,kFALSE);
3085  }
3086  }
3087  }
3088  if(part->eta < -dEtaGap / 2 )
3089  {
3090  // RFP in negative eta acceptance
3091  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3092  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3093  {
3094  dQcosNeg = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3095  dQsinNeg = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3096  fFlowVecQneg[iHarm][iPower] += TComplex(dQcosNeg,dQsinNeg,kFALSE);
3097  }
3098  }
3099  }
3100  } // endfor {tracks} particle loop
3101 
3102  // // filling local flow vectors to global flow vector arrays
3103  // for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3104  // for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3105  // {
3106  // fFlowVecQpos[iHarm][iPower] = TComplex(dQcosPos[iHarm][iPower],dQsinPos[iHarm][iPower],kFALSE);
3107  // if(dEtaGap > -1)
3108  // fFlowVecQneg[iHarm][iPower] = TComplex(dQcosNeg[iHarm][iPower],dQsinNeg[iHarm][iPower],kFALSE);
3109  // }
3110 
3111  // printf("RFPs EtaGap %g : number %g (pos) %g (neg) \n", dEtaGap,fFlowVecQpos[0][0].Re(),fFlowVecQneg[0][0].Re());
3112  return;
3113 }
3114 //_____________________________________________________________________________
3115 void AliAnalysisTaskFlowModes::FillPOIsVectors(const Short_t iEtaGapIndex, const PartSpecies species, const Short_t iMassIndex)
3116 {
3117  // Filling p,q and s flow vectors with POIs (given by species) for differential flow calculation
3118  // *************************************************************
3119 
3120  if(species == kUnknown) return;
3121 
3122  Double_t dNUAWeight = 1.; // for generic framework != 1
3123  Double_t dNUEWeight = 1.; // for generic framework != 1
3124 
3125  // clearing output (global) flow vectors
3129 
3130  std::vector<FlowPart>* vector = 0x0;
3131  //TH3D* hist = 0x0;
3132  //Double_t dMassLow = 0, dMassHigh = 0;
3133  TH3D* h3NUAWeights = 0x0;
3134  TH1F* hNUEWeights = 0x0;
3135 
3136  Int_t c = -1;
3137  if(fIndexCentrality>=0 && fIndexCentrality<5) c = 0;
3138  if(fIndexCentrality>=5 && fIndexCentrality<10) c = 1;
3139  if(fIndexCentrality>=10 && fIndexCentrality<20) c = 2;
3140  if(fIndexCentrality>=20 && fIndexCentrality<30) c = 3;
3141  if(fIndexCentrality>=30 && fIndexCentrality<40) c = 4;
3142  if(fIndexCentrality>=40 && fIndexCentrality<50) c = 5;
3143 
3144  // swich based on species
3145  switch (species)
3146  {
3147  case kCharged:
3148  vector = fVectorCharged;
3149  break;
3150 
3151  case kPion:
3152  vector = fVectorPion;
3153  break;
3154 
3155  case kKaon:
3156  vector = fVectorKaon;
3157  break;
3158 
3159  case kProton:
3160  vector = fVectorProton;
3161  break;
3162 
3163  default:
3164  ::Error("FillPOIsVectors","Selected species unknown.");
3165  return;
3166  }//switch(species)
3167 
3168  const Double_t dEtaGap = fEtaGap[iEtaGapIndex];
3169  //const Double_t dMass = (dMassLow+dMassHigh)/2;
3170 
3171  Short_t iPtBin = 0;
3172  Double_t dCos = 0, dSin = 0;
3173 
3174  for (auto part = vector->begin(); part != vector->end(); part++)
3175  {
3176 
3177  // swich based on species
3178  switch (species)
3179  {
3180  case kCharged:
3181  if(fFlowUseNUAWeights) {
3182  if(part->charge >0) h3NUAWeights = fh3NUAWeightChargedPlus;
3183  if(part->charge <0) h3NUAWeights = fh3NUAWeightChargedMinus;
3184  }
3185  if(fFlowUseNUEWeights) {
3186  if(part->charge >0 && c>=0 && c<6) hNUEWeights = fhNUEWeightChargedPlus[c];
3187  if(part->charge <0 && c>=0 && c<6) hNUEWeights = fhNUEWeightChargedMinus[c];
3188  }
3189  break;
3190 
3191  case kPion:
3192  if(fFlowUseNUAWeights) {
3193  if(part->charge >0) h3NUAWeights = fh3NUAWeightPionPlus;
3194  if(part->charge <0) h3NUAWeights = fh3NUAWeightPionMinus;
3195  }
3196  if(fFlowUseNUEWeights) {
3197  if(part->charge >0 && c>=0 && c<6) hNUEWeights = fhNUEWeightPionPlus[c];
3198  if(part->charge <0 && c>=0 && c<6) hNUEWeights = fhNUEWeightPionMinus[c];
3199  }
3200  break;
3201 
3202  case kKaon:
3203  if(fFlowUseNUAWeights) {
3204  if(part->charge >0) h3NUAWeights = fh3NUAWeightKaonPlus;
3205  if(part->charge <0) h3NUAWeights = fh3NUAWeightKaonMinus;
3206  }
3207  if(fFlowUseNUEWeights) {
3208  if(part->charge >0 && c>=0 && c<6) hNUEWeights = fhNUEWeightKaonPlus[c];
3209  if(part->charge <0 && c>=0 && c<6) hNUEWeights = fhNUEWeightKaonMinus[c];
3210  }
3211  break;
3212 
3213  case kProton:
3214  if(fFlowUseNUAWeights) {
3215  if(part->charge >0) h3NUAWeights = fh3NUAWeightProtonPlus;
3216  if(part->charge <0) h3NUAWeights = fh3NUAWeightProtonMinus;
3217  }
3218  if(fFlowUseNUEWeights) {
3219  if(part->charge >0 && c>=0 && c<6) hNUEWeights = fhNUEWeightProtonPlus[c];
3220  if(part->charge <0 && c>=0 && c<6) hNUEWeights = fhNUEWeightProtonMinus[c];
3221  }
3222  break;
3223 
3224  default:
3225  ::Error("FillPOIsVectors","Selected species unknown.");
3226  return;
3227  }//switch(species)
3228 
3229  if(fFlowUseNUAWeights && !h3NUAWeights) { ::Error("FillPOIsVectors","Histogram with NUA weights not found."); continue; }
3230  if(fFlowUseNUEWeights && !hNUEWeights) { ::Error("FillPOIsVectors","Histogram with NUE weights not found."); return; }
3231 
3232  // checking species of used particles (just for double checking purpose)
3233  if( part->species != species)
3234  {
3235  ::Warning("FillPOIsVectors","Unexpected part. species (%d) in selected sample (expected %d)",part->species,species);
3236  continue;
3237  }//endif{part->species != species}
3238 
3239  // assign iPtBin based on particle momenta
3240  iPtBin = GetPOIsPtBinIndex(part->pt);
3241 
3242  // 0-ing variables
3243  dCos = 0;
3244  dSin = 0;
3245 
3246  // POIs candidate passing all criteria: start filling flow vectors
3247 
3248  // loading weights if needed
3249  if(fFlowUseNUAWeights && h3NUAWeights)
3250  {
3251  dNUAWeight = h3NUAWeights->GetBinContent(h3NUAWeights->FindBin(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ()));
3252  if(dNUAWeight <= 0){ dNUAWeight = 1.;}
3253  if(iEtaGapIndex == 0){
3254  switch (species)
3255  {
3256  case kCharged:
3257  fh3AfterNUAWeightsCharged->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3258  break;
3259  case kPion:
3260  fh3AfterNUAWeightsPion->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3261  break;
3262  case kKaon:
3263  fh3AfterNUAWeightsKaon->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3264  break;
3265  case kProton:
3266  fh3AfterNUAWeightsProton->Fill(part->phi,part->eta,fEventAOD->GetPrimaryVertex()->GetZ(), dNUAWeight);
3267  break;
3268  default:
3269  ::Error("Fill fh3AfterNUAWeights","Selected species unknown.");
3270  return;
3271  }
3272  }
3273  }//endif{fFlowUseNUAWeights && h3NUAWeights}
3274 
3275  if(fFlowUseNUEWeights && hNUEWeights)
3276  {
3277  dNUEWeight = hNUEWeights->GetBinContent(hNUEWeights->FindBin(part->pt));
3278  if(dNUEWeight <= 0) dNUEWeight = 1.;
3279  if(iEtaGapIndex == 0){
3280  switch (species)
3281  {
3282  case kCharged:
3283  fhAfterNUEWeightsCharged->Fill(part->pt, dNUEWeight);
3284  break;
3285  case kPion:
3286  fhAfterNUEWeightsPion->Fill(part->pt, dNUEWeight);
3287  break;
3288  case kKaon:
3289  fhAfterNUEWeightsKaon->Fill(part->pt, dNUEWeight);
3290  break;
3291  case kProton:
3292  fhAfterNUEWeightsProton->Fill(part->pt, dNUEWeight);
3293  break;
3294  default:
3295  ::Error("Fill fhAfterNUEWeights","Selected species unknown.");
3296  return;
3297  }
3298  }
3299  }//endif{fFlowUseNUEWeights && hNUEWeights}
3300 
3301  if(part->eta > dEtaGap / 2 )
3302  {
3303  // particle in positive eta acceptance
3304  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3305  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3306  {
3307  dCos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3308  dSin = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3309  fFlowVecPpos[iHarm][iPower][iPtBin] += TComplex(dCos,dSin,kFALSE);
3310  }//endfor{Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++}
3311  }//endfor{Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++}
3312  }//endif{part->eta > dEtaGap / 2 }
3313  if(part->eta < -dEtaGap / 2 ){
3314  // particle in negative eta acceptance
3315  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++){
3316  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3317  {
3318  dCos = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Cos(iHarm * part->phi);
3319  dSin = TMath::Power(dNUAWeight*dNUEWeight,iPower) * TMath::Sin(iHarm * part->phi);
3320  fFlowVecPneg[iHarm][iPower][iPtBin] += TComplex(dCos,dSin,kFALSE);
3321  }//endfor{Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++}
3322  }//endfor{Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++}
3323  }//endif{part->eta < -dEtaGap / 2 }
3324  } // endfor {tracks}
3325  return;
3326 }
3327 //_____________________________________________________________________________
3329 {
3330  // Return POIs pT bin index based on pT value
3331  // *************************************************************
3332  const Double_t dPtBinWidth = (fFlowPOIsPtMax - fFlowPOIsPtMin) / fFlowPOIsPtNumBins;
3333  // printf("Pt %g | index %d\n",pt,(Short_t) (pt / dPtBinWidth) );
3334  return (Short_t) (pt / dPtBinWidth);
3335 }
3336 //_____________________________________________________________________________
3338 {
3339  // Reset RFPs (Q) array values to TComplex(0,0,kFALSE) for given array
3340  // *************************************************************
3341  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3342  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3343  array[iHarm][iPower] = TComplex(0,0,kFALSE);
3344  return;
3345 }
3346 //_____________________________________________________________________________
3348 {
3349  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3350  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3351  for(Short_t iPt(0); iPt < fFlowPOIsPtNumBins; iPt++)
3352  array[iHarm][iPower][iPt] = TComplex(0,0,kFALSE);
3353  return;
3354 }
3355 //_____________________________________________________________________________
3357 {
3358  // List all values of given flow vector TComplex array
3359  // *************************************************************
3360  printf(" ### Listing (TComplex) flow vector array ###########################\n");
3361  for(Short_t iHarm(0); iHarm < fFlowNumHarmonicsMax; iHarm++)
3362  {
3363  printf("Harm %d (power):",iHarm);
3364  for(Short_t iPower(0); iPower < fFlowNumWeightPowersMax; iPower++)
3365  {
3366  printf("|(%d) %g+%g(i)",iPower, array[iHarm][iPower].Re(), array[iHarm][iPower].Im());
3367  }
3368  printf("\n");
3369  }
3370  return;
3371 }
3372 //_____________________________________________________________________________
3374 {
3375  // Assessing sampling index based on generated random number
3376  // returns centrality index
3377  // *************************************************************
3378 
3379  Short_t index = 0x0;
3380 
3381  if(fSampling && fNumSamples > 1)
3382  {
3383  TRandom3 rr(0);
3384  Double_t ranNum = rr.Rndm(); // getting random number in (0,1)
3385  Double_t generated = ranNum * fNumSamples; // getting random number in range (0, fNumSamples)
3386  // finding right index for sampling based on generated number and total number of samples
3387  for(Short_t i(0); i < fNumSamples; i++)
3388  {
3389  if(generated < (i+1) )
3390  {
3391  index = i;
3392  break;
3393  }
3394  }
3395  }
3396 
3397  return index;
3398 }
3399 
3400 //_____________________________________________________________________________
3402 {
3403  // Estimating centrality percentile based on selected estimator.
3404  // (Default) If no multiplicity estimator is specified (fMultEstimator == '' || Charged), percentile is estimated as number of selected / filtered charged tracks.
3405  // If a valid multiplicity estimator is specified, centrality percentile is estimated via AliMultSelection
3406  // otherwise -1 is returned (and event is skipped)
3407  // *************************************************************
3408  fMultEstimator.ToUpper();
3409  //cout<<"++++++++++++++++++GetCentralityIndex+++++++++++++++++"<<endl;
3410  if(
3411  fMultEstimator.EqualTo("V0A") || fMultEstimator.EqualTo("V0C") || fMultEstimator.EqualTo("V0M") ||
3412  fMultEstimator.EqualTo("CL0") || fMultEstimator.EqualTo("CL1") || fMultEstimator.EqualTo("ZNA") ||
3413  fMultEstimator.EqualTo("ZNC") || fMultEstimator.EqualTo("TRK") ){
3414  //cout<<"++++++++++++++++++CentralityIndex+++++++++++++++++"<<endl;
3415 
3416  // some of supported AliMultSelection estimators (listed above)
3417  Float_t dPercentile = 300;
3418 
3419  // checking AliMultSelection
3420  AliMultSelection* multSelection = 0x0;
3421  multSelection = (AliMultSelection*) fEventAOD->FindListObject("MultSelection");
3422  if(!multSelection) { AliError("AliMultSelection object not found! Returning -1"); return -1;}
3423 
3424  dPercentile = multSelection->GetMultiplicityPercentile(fMultEstimator.Data());
3425  //cout<<"fMultEstimator.Data()"<<fMultEstimator.Data()<<endl;
3426  //cout<<"dPercentile = "<<dPercentile<<endl;
3428  if(dPercentile > 100 || dPercentile < 0) { AliWarning("Centrality percentile estimated not within 0-100 range. Returning -1"); return -1;}
3429  else { return dPercentile;}
3430  }
3431  if(!fFullCentralityRange){
3432  if(dPercentile > 100 || dPercentile <50) { AliWarning("Centrality percentile estimated not within 50-100 range. Returning -1"); return -1;}
3433  else { return dPercentile;}
3434  }
3435  }
3436  else if(fMultEstimator.EqualTo("") || fMultEstimator.EqualTo("CHARGED"))
3437  {
3438  // assigning centrality based on number of selected charged tracks
3439  return fVectorCharged->size();
3440  }
3441  else
3442  {
3443  AliWarning(Form("Multiplicity estimator '%s' not supported. Returning -1\n",fMultEstimator.Data()));
3444  return -1;
3445  }
3446 
3447 
3448  return -1;
3449 }
3450 //____________________________
3451 Double_t AliAnalysisTaskFlowModes::GetWDist(const AliAODVertex* v0, const AliAODVertex* v1) {
3452 // calculate sqrt of weighted distance to other vertex
3453  if (!v0 || !v1) {
3454  printf("One of vertices is not valid\n");
3455  return kFALSE;
3456  }
3457  static TMatrixDSym vVb(3);
3458  double dist = -1;
3459  double dx = v0->GetX()-v1->GetX();
3460  double dy = v0->GetY()-v1->GetY();
3461  double dz = v0->GetZ()-v1->GetZ();
3462  double cov0[6],cov1[6];
3463  v0->GetCovarianceMatrix(cov0);
3464  v1->GetCovarianceMatrix(cov1);
3465  vVb(0,0) = cov0[0]+cov1[0];
3466  vVb(1,1) = cov0[2]+cov1[2];
3467  vVb(2,2) = cov0[5]+cov1[5];
3468  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
3469  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
3470  vVb.InvertFast();
3471  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
3472  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;
3473  return dist>0 ? TMath::Sqrt(dist) : -1;
3474 }
3475 //_________________________________________________
3477 {
3478  // called on end of task, after all events are processed
3479  // *************************************************************
3480 
3481  return;
3482 }
3483 //_____________________________________________________________________________
3484 // Set of methods returning given complex flow vector based on flow harmonics (n) and weight power indexes (p)
3485 // a la General Framework implementation.
3486 // Q: flow vector of RFPs (with/out eta gap)
3487 // P: flow vector of POIs (with/out eta gap) (in usual notation p)
3488 // S: flow vector of overlaping RFPs and POIs (in usual notation q)
3489 
3490 TComplex AliAnalysisTaskFlowModes::Q(const Short_t n, const Short_t p)
3491 {
3492  if (n < 0) return TComplex::Conjugate(fFlowVecQpos[-n][p]);
3493  else return fFlowVecQpos[n][p];
3494 }
3495 //____________________________________________________________________
3497 {
3498  if (n < 0) return TComplex::Conjugate(fFlowVecQpos[-n][p]);
3499  else return fFlowVecQpos[n][p];
3500 }
3501 //____________________________________________________________________
3503 {
3504  if(n < 0) return TComplex::Conjugate(fFlowVecQneg[-n][p]);
3505  else return fFlowVecQneg[n][p];
3506 }
3507 //____________________________________________________________________
3508 TComplex AliAnalysisTaskFlowModes::P(const Short_t n, const Short_t p, const Short_t pt)
3509 {
3510  if(n < 0) return TComplex::Conjugate(fFlowVecPpos[-n][p][pt]);
3511  else return fFlowVecPpos[n][p][pt];
3512 }
3513 //____________________________________________________________________
3514 TComplex AliAnalysisTaskFlowModes::PGapPos(const Short_t n, const Short_t p, const Short_t pt)
3515 {
3516  if(n < 0) return TComplex::Conjugate(fFlowVecPpos[-n][p][pt]);
3517  else return fFlowVecPpos[n][p][pt];
3518 }
3519 //____________________________________________________________________
3520 TComplex AliAnalysisTaskFlowModes::PGapNeg(const Short_t n, const Short_t p, const Short_t pt)
3521 {
3522  if(n < 0) return TComplex::Conjugate(fFlowVecPneg[-n][p][pt]);
3523  else return fFlowVecPneg[n][p][pt];
3524 }
3525 //____________________________________________________________________
3526 TComplex AliAnalysisTaskFlowModes::S(const Short_t n, const Short_t p, const Short_t pt)
3527 {
3528  if(n < 0) return TComplex::Conjugate(fFlowVecS[-n][p][pt]);
3529  else return fFlowVecS[n][p][pt];
3530 }
3531 //____________________________________________________________________
3532 
3533 // Set of flow calculation methods for cumulants of different orders with/out eta gap
3534 
3535 TComplex AliAnalysisTaskFlowModes::Two(const Short_t n1, const Short_t n2)
3536 {
3537  TComplex formula = Q(n1,1)*Q(n2,1) - Q(n1+n2,2);
3538  return formula;
3539 }
3540 //____________________________________________________________________
3542 {
3543  TComplex formula = QGapPos(n1,1)*QGapNeg(n2,1);
3544  return formula;
3545 }
3546 //____________________________________________________________________
3547 TComplex AliAnalysisTaskFlowModes::TwoDiff(const Short_t n1, const Short_t n2, const Short_t pt)
3548 {
3549  TComplex formula = P(n1,1,pt)*Q(n2,1) - S(n1+n2,1,pt);
3550  return formula;
3551 }
3552 //____________________________________________________________________
3553 TComplex AliAnalysisTaskFlowModes::TwoDiffGapPos(const Short_t n1, const Short_t n2, const Short_t pt)
3554 {
3555  TComplex formula = PGapPos(n1,1,pt)*QGapNeg(n2,1);
3556  return formula;
3557 }
3558 //____________________________________________________________________
3559 TComplex AliAnalysisTaskFlowModes::TwoDiffGapNeg(const Short_t n1, const Short_t n2, const Short_t pt)
3560 {
3561  TComplex formula = PGapNeg(n1,1,pt)*QGapPos(n2,1);
3562  return formula;
3563 }
3564 //____________________________________________________________________
3565 TComplex AliAnalysisTaskFlowModes::ThreeDiffGapPos(const Short_t n1, const Short_t n2, const Short_t n3,const Short_t pt)
3566 {
3567  TComplex formula = PGapPos(n1,1,pt)*QGapNeg(n2,1)*QGapNeg(n3,1)- PGapPos(n1,1,pt)*QGapNeg(n2+n3,2);
3568  return formula;
3569 }
3570 //____________________________________________________________________
3571 TComplex AliAnalysisTaskFlowModes::ThreeDiffGapNeg(const Short_t n1, const Short_t n2, const Short_t n3,const Short_t pt)
3572 {
3573  TComplex formula = PGapNeg(n1,1,pt)*QGapPos(n2,1)*QGapPos(n3,1)- PGapNeg(n1,1,pt)*QGapPos(n2+n3,2);
3574  return formula;
3575 }
3576 //____________________________________________________________________
3577 TComplex AliAnalysisTaskFlowModes::Four(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4)
3578 {
3579  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)
3580  - 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)
3581  + 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)
3582  + 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)
3583  + 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);
3584  return formula;
3585 }
3586 // //____________________________________________________________________
3587 TComplex AliAnalysisTaskFlowModes::FourDiff(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t pt)
3588 {
3589  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)
3590  - 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)
3591  + 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)
3592  + 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)
3593  + 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);
3594  return formula;
3595 }
3596 //____________________________________________________________________
3597 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
3598 {
3599  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);
3600  //TComplex *out = (TComplex*) &formula;
3601  return formula;
3602 }
3603 //____________________________________________________________________
3604 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
3605 {
3606  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);
3607  //TComplex *out = (TComplex*) &formula;
3608  return formula;
3609 }
3610 //____________________________________________________________________
3611 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
3612 {
3613  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);
3614  return formula;
3615 }
3616 //____________________________________________________________________
3617 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
3618 {
3619  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);
3620  return formula;
3621 }
3622 //____________________________________________________________________
3623 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
3624 {
3625  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);
3626  return formula;
3627 }
3628 //____________________________________________________________________
3629 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
3630 {
3631 
3632  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);
3633  return formula;
3634 }
3635 
3636 //-----------------------------------------------------------------------
3638  //set priors for the bayesian pid selection
3639  fCurrCentr = centrCur;
3640 
3641  fBinLimitPID[0] = 0.300000;
3642  fBinLimitPID[1] = 0.400000;
3643  fBinLimitPID[2] = 0.500000;
3644  fBinLimitPID[3] = 0.600000;
3645  fBinLimitPID[4] = 0.700000;
3646  fBinLimitPID[5] = 0.800000;
3647  fBinLimitPID[6] = 0.900000;
3648  fBinLimitPID[7] = 1.000000;
3649  fBinLimitPID[8] = 1.200000;
3650  fBinLimitPID[9] = 1.400000;
3651  fBinLimitPID[10] = 1.600000;
3652  fBinLimitPID[11] = 1.800000;
3653  fBinLimitPID[12] = 2.000000;
3654  fBinLimitPID[13] = 2.200000;
3655  fBinLimitPID[14] = 2.400000;
3656  fBinLimitPID[15] = 2.600000;
3657  fBinLimitPID[16] = 2.800000;
3658  fBinLimitPID[17] = 3.000000;
3659 
3660  for(Int_t i=18;i<fgkPIDptBin;i++){
3661  fBinLimitPID[i] = 3.0 + 0.2 * (i-17);
3662  }
3663 
3664  // 0-10%
3665  if(centrCur < 10){
3666  fC[0][0] = 0.005;
3667  fC[0][1] = 0.005;
3668  fC[0][2] = 1.0000;
3669  fC[0][3] = 0.0010;
3670  fC[0][4] = 0.0010;
3671 
3672  fC[1][0] = 0.005;
3673  fC[1][1] = 0.005;
3674  fC[1][2] = 1.0000;
3675  fC[1][3] = 0.0168;
3676  fC[1][4] = 0.0122;
3677 
3678  fC[2][0] = 0.005;
3679  fC[2][1] = 0.005;
3680  fC[2][2] = 1.0000;
3681  fC[2][3] = 0.0272;
3682  fC[2][4] = 0.0070;
3683 
3684  fC[3][0] = 0.005;
3685  fC[3][1] = 0.005;
3686  fC[3][2] = 1.0000;
3687  fC[3][3] = 0.0562;
3688  fC[3][4] = 0.0258;
3689 
3690  fC[4][0] = 0.005;
3691  fC[4][1] = 0.005;
3692  fC[4][2] = 1.0000;
3693  fC[4][3] = 0.0861;
3694  fC[4][4] = 0.0496;
3695 
3696  fC[5][0] = 0.005;
3697  fC[5][1] = 0.005;
3698  fC[5][2] = 1.0000;
3699  fC[5][3] = 0.1168;
3700  fC[5][4] = 0.0740;
3701 
3702  fC[6][0] = 0.005;
3703  fC[6][1] = 0.005;
3704  fC[6][2] = 1.0000;
3705  fC[6][3] = 0.1476;
3706  fC[6][4] = 0.0998;
3707 
3708  fC[7][0] = 0.005;
3709  fC[7][1] = 0.005;
3710  fC[7][2] = 1.0000;
3711  fC[7][3] = 0.1810;
3712  fC[7][4] = 0.1296;
3713 
3714  fC[8][0] = 0.005;
3715  fC[8][1] = 0.005;
3716  fC[8][2] = 1.0000;
3717  fC[8][3] = 0.2240;
3718  fC[8][4] = 0.1827;
3719 
3720  fC[9][0] = 0.005;
3721  fC[9][1] = 0.005;
3722  fC[9][2] = 1.0000;
3723  fC[9][3] = 0.2812;
3724  fC[9][4] = 0.2699;
3725 
3726  fC[10][0] = 0.005;
3727  fC[10][1] = 0.005;
3728  fC[10][2] = 1.0000;
3729  fC[10][3] = 0.3328;
3730  fC[10][4] = 0.3714;
3731 
3732  fC[11][0] = 0.005;
3733  fC[11][1] = 0.005;
3734  fC[11][2] = 1.0000;
3735  fC[11][3] = 0.3780;
3736  fC[11][4] = 0.4810;
3737 
3738  fC[12][0] = 0.005;
3739  fC[12][1] = 0.005;
3740  fC[12][2] = 1.0000;
3741  fC[12][3] = 0.4125;
3742  fC[12][4] = 0.5771;
3743 
3744  fC[13][0] = 0.005;
3745  fC[13][1] = 0.005;
3746  fC[13][2] = 1.0000;
3747  fC[13][3] = 0.4486;
3748  fC[13][4] = 0.6799;
3749 
3750  fC[14][0] = 0.005;
3751  fC[14][1] = 0.005;
3752  fC[14][2] = 1.0000;
3753  fC[14][3] = 0.4840;
3754  fC[14][4] = 0.7668;
3755 
3756  fC[15][0] = 0.005;
3757  fC[15][1] = 0.005;
3758  fC[15][2] = 1.0000;
3759  fC[15][3] = 0.4971;
3760  fC[15][4] = 0.8288;
3761 
3762  fC[16][0] = 0.005;
3763  fC[16][1] = 0.005;
3764  fC[16][2] = 1.0000;
3765  fC[16][3] = 0.4956;
3766  fC[16][4] = 0.8653;
3767 
3768  fC[17][0] = 0.005;
3769  fC[17][1] = 0.005;
3770  fC[17][2] = 1.0000;
3771  fC[17][3] = 0.5173;
3772  fC[17][4] = 0.9059;
3773 
3774  for(Int_t i=18;i<fgkPIDptBin;i++){
3775  fC[i][0] = fC[17][0];
3776  fC[i][1] = fC[17][1];
3777  fC[i][2] = fC[17][2];
3778  fC[i][3] = fC[17][3];
3779  fC[i][4] = fC[17][4];
3780  }
3781  }
3782  // 10-20%
3783  else if(centrCur < 20){
3784  fC[0][0] = 0.005;
3785  fC[0][1] = 0.005;
3786  fC[0][2] = 1.0000;
3787  fC[0][3] = 0.0010;
3788  fC[0][4] = 0.0010;
3789 
3790  fC[1][0] = 0.005;
3791  fC[1][1] = 0.005;
3792  fC[1][2] = 1.0000;
3793  fC[1][3] = 0.0132;
3794  fC[1][4] = 0.0088;
3795 
3796  fC[2][0] = 0.005;
3797  fC[2][1] = 0.005;
3798  fC[2][2] = 1.0000;
3799  fC[2][3] = 0.0283;
3800  fC[2][4] = 0.0068;
3801 
3802  fC[3][0] = 0.005;
3803  fC[3][1] = 0.005;
3804  fC[3][2] = 1.0000;
3805  fC[3][3] = 0.0577;
3806  fC[3][4] = 0.0279;
3807 
3808  fC[4][0] = 0.005;
3809  fC[4][1] = 0.005;
3810  fC[4][2] = 1.0000;
3811  fC[4][3] = 0.0884;
3812  fC[4][4] = 0.0534;
3813 
3814  fC[5][0] = 0.005;
3815  fC[5][1] = 0.005;
3816  fC[5][2] = 1.0000;
3817  fC[5][3] = 0.1179;
3818  fC[5][4] = 0.0794;
3819 
3820  fC[6][0] = 0.005;
3821  fC[6][1] = 0.005;
3822  fC[6][2] = 1.0000;
3823  fC[6][3] = 0.1480;
3824  fC[6][4] = 0.1058;
3825 
3826  fC[7][0] = 0.005;
3827  fC[7][1] = 0.005;
3828  fC[7][2] = 1.0000;
3829  fC[7][3] = 0.1807;
3830  fC[7][4] = 0.1366;
3831 
3832  fC[8][0] = 0.005;
3833  fC[8][1] = 0.005;
3834  fC[8][2] = 1.0000;
3835  fC[8][3] = 0.2219;
3836  fC[8][4] = 0.1891;
3837 
3838  fC[9][0] = 0.005;
3839  fC[9][1] = 0.005;
3840  fC[9][2] = 1.0000;
3841  fC[9][3] = 0.2804;
3842  fC[9][4] = 0.2730;
3843 
3844  fC[10][0] = 0.005;
3845  fC[10][1] = 0.005;
3846  fC[10][2] = 1.0000;
3847  fC[10][3] = 0.3283;
3848  fC[10][4] = 0.3660;
3849 
3850  fC[11][0] = 0.005;
3851  fC[11][1] = 0.005;
3852  fC[11][2] = 1.0000;
3853  fC[11][3] = 0.3710;
3854  fC[11][4] = 0.4647;
3855 
3856  fC[12][0] = 0.005;
3857  fC[12][1] = 0.005;
3858  fC[12][2] = 1.0000;
3859  fC[12][3] = 0.4093;
3860  fC[12][4] = 0.5566;
3861 
3862  fC[13][0] = 0.005;
3863  fC[13][1] = 0.005;
3864  fC[13][2] = 1.0000;
3865  fC[13][3] = 0.4302;
3866  fC[13][4] = 0.6410;
3867 
3868  fC[14][0] = 0.005;
3869  fC[14][1] = 0.005;
3870  fC[14][2] = 1.0000;
3871  fC[14][3] = 0.4649;
3872  fC[14][4] = 0.7055;
3873 
3874  fC[15][0] = 0.005;
3875  fC[15][1] = 0.005;
3876  fC[15][2] = 1.0000;
3877  fC[15][3] = 0.4523;
3878  fC[15][4] = 0.7440;
3879 
3880  fC[16][0] = 0.005;
3881  fC[16][1] = 0.005;
3882  fC[16][2] = 1.0000;
3883  fC[16][3] = 0.4591;
3884  fC[16][4] = 0.7799;
3885 
3886  fC[17][0] = 0.005;
3887  fC[17][1] = 0.005;
3888  fC[17][2] = 1.0000;
3889  fC[17][3] = 0.4804;
3890  fC[17][4] = 0.8218;
3891 
3892  for(Int_t i=18;i<fgkPIDptBin;i++){
3893  fC[i][0] = fC[17][0];
3894  fC[i][1] = fC[17][1];
3895  fC[i][2] = fC[17][2];
3896  fC[i][3] = fC[17][3];
3897  fC[i][4] = fC[17][4];
3898  }
3899  }
3900  // 20-30%
3901  else if(centrCur < 30){
3902  fC[0][0] = 0.005;
3903  fC[0][1] = 0.005;
3904  fC[0][2] = 1.0000;
3905  fC[0][3] = 0.0010;
3906  fC[0][4] = 0.0010;
3907 
3908  fC[1][0] = 0.005;
3909  fC[1][1] = 0.005;
3910  fC[1][2] = 1.0000;
3911  fC[1][3] = 0.0102;
3912  fC[1][4] = 0.0064;
3913 
3914  fC[2][0] = 0.005;
3915  fC[2][1] = 0.005;
3916  fC[2][2] = 1.0000;
3917  fC[2][3] = 0.0292;
3918  fC[2][4] = 0.0066;
3919 
3920  fC[3][0] = 0.005;
3921  fC[3][1] = 0.005;
3922  fC[3][2] = 1.0000;
3923  fC[3][3] = 0.0597;
3924  fC[3][4] = 0.0296;
3925 
3926  fC[4][0] = 0.005;
3927  fC[4][1] = 0.005;
3928  fC[4][2] = 1.0000;
3929  fC[4][3] = 0.0900;
3930  fC[4][4] = 0.0589;
3931 
3932  fC[5][0] = 0.005;
3933  fC[5][1] = 0.005;
3934  fC[5][2] = 1.0000;
3935  fC[5][3] = 0.1199;
3936  fC[5][4] = 0.0859;
3937 
3938  fC[6][0] = 0.005;
3939  fC[6][1] = 0.005;
3940  fC[6][2] = 1.0000;
3941  fC[6][3] = 0.1505;
3942  fC[6][4] = 0.1141;
3943 
3944  fC[7][0] = 0.005;
3945  fC[7][1] = 0.005;
3946  fC[7][2] = 1.0000;
3947  fC[7][3] = 0.1805;
3948  fC[7][4] = 0.1454;
3949 
3950  fC[8][0] = 0.005;
3951  fC[8][1] = 0.005;
3952  fC[8][2] = 1.0000;
3953  fC[8][3] = 0.2221;
3954  fC[8][4] = 0.2004;
3955 
3956  fC[9][0] = 0.005;
3957  fC[9][1] = 0.005;
3958  fC[9][2] = 1.0000;
3959  fC[9][3] = 0.2796;
3960  fC[9][4] = 0.2838;
3961 
3962  fC[10][0] = 0.005;
3963  fC[10][1] = 0.005;
3964  fC[10][2] = 1.0000;
3965  fC[10][3] = 0.3271;
3966  fC[10][4] = 0.3682;
3967 
3968  fC[11][0] = 0.005;
3969  fC[11][1] = 0.005;
3970  fC[11][2] = 1.0000;
3971  fC[11][3] = 0.3648;
3972  fC[11][4] = 0.4509;
3973 
3974  fC[12][0] = 0.005;
3975  fC[12][1] = 0.005;
3976  fC[12][2] = 1.0000;
3977  fC[12][3] = 0.3988;
3978  fC[12][4] = 0.5339;
3979 
3980  fC[13][0] = 0.005;
3981  fC[13][1] = 0.005;
3982  fC[13][2] = 1.0000;
3983  fC[13][3] = 0.4315;
3984  fC[13][4] = 0.5995;
3985 
3986  fC[14][0] = 0.005;
3987  fC[14][1] = 0.005;
3988  fC[14][2] = 1.0000;
3989  fC[14][3] = 0.4548;
3990  fC[14][4] = 0.6612;
3991 
3992  fC[15][0] = 0.005;
3993  fC[15][1] = 0.005;
3994  fC[15][2] = 1.0000;
3995  fC[15][3] = 0.4744;
3996  fC[15][4] = 0.7060;
3997 
3998  fC[16][0] = 0.005;
3999  fC[16][1] = 0.005;
4000  fC[16][2] = 1.0000;
4001  fC[16][3] = 0.4899;
4002  fC[16][4] = 0.7388;
4003 
4004  fC[17][0] = 0.005;
4005  fC[17][1] = 0.005;
4006  fC[17][2] = 1.0000;
4007  fC[17][3] = 0.4411;
4008  fC[17][4] = 0.7293;
4009 
4010  for(Int_t i=18;i<fgkPIDptBin;i++){
4011  fC[i][0] = fC[17][0];
4012  fC[i][1] = fC[17][1];
4013  fC[i][2] = fC[17][2];
4014  fC[i][3] = fC[17][3];
4015  fC[i][4] = fC[17][4];
4016  }
4017  }
4018  // 30-40%
4019  else if(centrCur < 40){
4020  fC[0][0] = 0.005;
4021  fC[0][1] = 0.005;
4022  fC[0][2] = 1.0000;
4023  fC[0][3] = 0.0010;
4024  fC[0][4] = 0.0010;
4025 
4026  fC[1][0] = 0.005;
4027  fC[1][1] = 0.005;
4028  fC[1][2] = 1.0000;
4029  fC[1][3] = 0.0102;
4030  fC[1][4] = 0.0048;
4031 
4032  fC[2][0] = 0.005;
4033  fC[2][1] = 0.005;
4034  fC[2][2] = 1.0000;
4035  fC[2][3] = 0.0306;
4036  fC[2][4] = 0.0079;
4037 
4038  fC[3][0] = 0.005;
4039  fC[3][1] = 0.005;
4040  fC[3][2] = 1.0000;
4041  fC[3][3] = 0.0617;
4042  fC[3][4] = 0.0338;
4043 
4044  fC[4][0] = 0.005;
4045  fC[4][1] = 0.005;
4046  fC[4][2] = 1.0000;
4047  fC[4][3] = 0.0920;
4048  fC[4][4] = 0.0652;
4049 
4050  fC[5][0] = 0.005;
4051  fC[5][1] = 0.005;
4052  fC[5][2] = 1.0000;
4053  fC[5][3] = 0.1211;
4054  fC[5][4] = 0.0955;
4055 
4056  fC[6][0] = 0.005;
4057  fC[6][1] = 0.005;
4058  fC[6][2] = 1.0000;
4059  fC[6][3] = 0.1496;
4060  fC[6][4] = 0.1242;
4061 
4062  fC[7][0] = 0.005;
4063  fC[7][1] = 0.005;
4064  fC[7][2] = 1.0000;
4065  fC[7][3] = 0.1807;
4066  fC[7][4] = 0.1576;
4067 
4068  fC[8][0] = 0.005;
4069  fC[8][1] = 0.005;
4070  fC[8][2] = 1.0000;
4071  fC[8][3] = 0.2195;
4072  fC[8][4] = 0.2097;
4073 
4074  fC[9][0] = 0.005;
4075  fC[9][1] = 0.005;
4076  fC[9][2] = 1.0000;
4077  fC[9][3] = 0.2732;
4078  fC[9][4] = 0.2884;
4079 
4080  fC[10][0] = 0.005;
4081  fC[10][1] = 0.005;
4082  fC[10][2] = 1.0000;
4083  fC[10][3] = 0.3204;
4084  fC[10][4] = 0.3679;
4085 
4086  fC[11][0] = 0.005;
4087  fC[11][1] = 0.005;
4088  fC[11][2] = 1.0000;
4089  fC[11][3] = 0.3564;
4090  fC[11][4] = 0.4449;
4091 
4092  fC[12][0] = 0.005;
4093  fC[12][1] = 0.005;
4094  fC[12][2] = 1.0000;
4095  fC[12][3] = 0.3791;
4096  fC[12][4] = 0.5052;
4097 
4098  fC[13][0] = 0.005;
4099  fC[13][1] = 0.005;
4100  fC[13][2] = 1.0000;
4101  fC[13][3] = 0.4062;
4102  fC[13][4] = 0.5647;
4103 
4104  fC[14][0] = 0.005;
4105  fC[14][1] = 0.005;
4106  fC[14][2] = 1.0000;
4107  fC[14][3] = 0.4234;
4108  fC[14][4] = 0.6203;
4109 
4110  fC[15][0] = 0.005;
4111  fC[15][1] = 0.005;
4112  fC[15][2] = 1.0000;
4113  fC[15][3] = 0.4441;
4114  fC[15][4] = 0.6381;
4115 
4116  fC[16][0] = 0.005;
4117  fC[16][1] = 0.005;
4118  fC[16][2] = 1.0000;
4119  fC[16][3] = 0.4629;
4120  fC[16][4] = 0.6496;
4121 
4122  fC[17][0] = 0.005;
4123  fC[17][1] = 0.005;
4124  fC[17][2] = 1.0000;
4125  fC[17][3] = 0.4293;
4126  fC[17][4] = 0.6491;
4127 
4128  for(Int_t i=18;i<fgkPIDptBin;i++){
4129  fC[i][0] = fC[17][0];
4130  fC[i][1] = fC[17][1];
4131  fC[i][2] = fC[17][2];
4132  fC[i][3] = fC[17][3];
4133  fC[i][4] = fC[17][4];
4134  }
4135  }
4136  // 40-50%
4137  else if(centrCur < 50){
4138  fC[0][0] = 0.005;
4139  fC[0][1] = 0.005;
4140  fC[0][2] = 1.0000;
4141  fC[0][3] = 0.0010;
4142  fC[0][4] = 0.0010;
4143 
4144  fC[1][0] = 0.005;
4145  fC[1][1] = 0.005;
4146  fC[1][2] = 1.0000;
4147  fC[1][3] = 0.0093;
4148  fC[1][4] = 0.0057;
4149 
4150  fC[2][0] = 0.005;
4151  fC[2][1] = 0.005;
4152  fC[2][2] = 1.0000;
4153  fC[2][3] = 0.0319;
4154  fC[2][4] = 0.0075;
4155 
4156  fC[3][0] = 0.005;
4157  fC[3][1] = 0.005;
4158  fC[3][2] = 1.0000;
4159  fC[3][3] = 0.0639;
4160  fC[3][4] = 0.0371;
4161 
4162  fC[4][0] = 0.005;
4163  fC[4][1] = 0.005;
4164  fC[4][2] = 1.0000;
4165  fC[4][3] = 0.0939;
4166  fC[4][4] = 0.0725;
4167 
4168  fC[5][0] = 0.005;
4169  fC[5][1] = 0.005;
4170  fC[5][2] = 1.0000;
4171  fC[5][3] = 0.1224;
4172  fC[5][4] = 0.1045;
4173 
4174  fC[6][0] = 0.005;
4175  fC[6][1] = 0.005;
4176  fC[6][2] = 1.0000;
4177  fC[6][3] = 0.1520;
4178  fC[6][4] = 0.1387;
4179 
4180  fC[7][0] = 0.005;
4181  fC[7][1] = 0.005;
4182  fC[7][2] = 1.0000;
4183  fC[7][3] = 0.1783;
4184  fC[7][4] = 0.1711;
4185 
4186  fC[8][0] = 0.005;
4187  fC[8][1] = 0.005;
4188  fC[8][2] = 1.0000;
4189  fC[8][3] = 0.2202;
4190  fC[8][4] = 0.2269;
4191 
4192  fC[9][0] = 0.005;
4193  fC[9][1] = 0.005;
4194  fC[9][2] = 1.0000;
4195  fC[9][3] = 0.2672;
4196  fC[9][4] = 0.2955;
4197 
4198  fC[10][0] = 0.005;
4199  fC[10][1] = 0.005;
4200  fC[10][2] = 1.0000;
4201  fC[10][3] = 0.3191;
4202  fC[10][4] = 0.3676;
4203 
4204  fC[11][0] = 0.005;
4205  fC[11][1] = 0.005;
4206  fC[11][2] = 1.0000;
4207  fC[11][3] = 0.3434;
4208  fC[11][4] = 0.4321;
4209 
4210  fC[12][0] = 0.005;
4211  fC[12][1] = 0.005;
4212  fC[12][2] = 1.0000;
4213  fC[12][3] = 0.3692;
4214  fC[12][4] = 0.4879;
4215 
4216  fC[13][0] = 0.005;
4217  fC[13][1] = 0.005;
4218  fC[13][2] = 1.0000;
4219  fC[13][3] = 0.3993;
4220  fC[13][4] = 0.5377;
4221 
4222  fC[14][0] = 0.005;
4223  fC[14][1] = 0.005;
4224  fC[14][2] = 1.0000;
4225  fC[14][3] = 0.3818;
4226  fC[14][4] = 0.5547;
4227 
4228  fC[15][0] = 0.005;
4229  fC[15][1] = 0.005;
4230  fC[15][2] = 1.0000;
4231  fC[15][3] = 0.4003;
4232  fC[15][4] = 0.5484;
4233 
4234  fC[16][0] = 0.005;
4235  fC[16][1] = 0.005;
4236  fC[16][2] = 1.0000;
4237  fC[16][3] = 0.4281;
4238  fC[16][4] = 0.5383;
4239 
4240  fC[17][0] = 0.005;
4241  fC[17][1] = 0.005;
4242  fC[17][2] = 1.0000;
4243  fC[17][3] = 0.3960;
4244  fC[17][4] = 0.5374;
4245 
4246  for(Int_t i=18;i<fgkPIDptBin;i++){
4247  fC[i][0] = fC[17][0];
4248  fC[i][1] = fC[17][1];
4249  fC[i][2] = fC[17][2];
4250  fC[i][3] = fC[17][3];
4251  fC[i][4] = fC[17][4];
4252  }
4253  }
4254  // 50-60%
4255  else if(centrCur < 60){
4256  fC[0][0] = 0.005;
4257  fC[0][1] = 0.005;
4258  fC[0][2] = 1.0000;
4259  fC[0][3] = 0.0010;
4260  fC[0][4] = 0.0010;
4261 
4262  fC[1][0] = 0.005;
4263  fC[1][1] = 0.005;
4264  fC[1][2] = 1.0000;
4265  fC[1][3] = 0.0076;
4266  fC[1][4] = 0.0032;
4267 
4268  fC[2][0] = 0.005;
4269  fC[2][1] = 0.005;
4270  fC[2][2] = 1.0000;
4271  fC[2][3] = 0.0329;
4272  fC[2][4] = 0.0085;
4273 
4274  fC[3][0] = 0.005;
4275  fC[3][1] = 0.005;
4276  fC[3][2] = 1.0000;
4277  fC[3][3] = 0.0653;
4278  fC[3][4] = 0.0423;
4279 
4280  fC[4][0] = 0.005;
4281  fC[4][1] = 0.005;
4282  fC[4][2] = 1.0000;
4283  fC[4][3] = 0.0923;
4284  fC[4][4] = 0.0813;
4285 
4286  fC[5][0] = 0.005;
4287  fC[5][1] = 0.005;
4288  fC[5][2] = 1.0000;
4289  fC[5][3] = 0.1219;
4290  fC[5][4] = 0.1161;
4291 
4292  fC[6][0] = 0.005;
4293  fC[6][1] = 0.005;
4294  fC[6][2] = 1.0000;
4295  fC[6][3] = 0.1519;
4296  fC[6][4] = 0.1520;
4297 
4298  fC[7][0] = 0.005;
4299  fC[7][1] = 0.005;
4300  fC[7][2] = 1.0000;
4301  fC[7][3] = 0.1763;
4302  fC[7][4] = 0.1858;
4303 
4304  fC[8][0] = 0.005;
4305  fC[8][1] = 0.005;
4306  fC[8][2] = 1.0000;
4307  fC[8][3] = 0.2178;
4308  fC[8][4] = 0.2385;
4309 
4310  fC[9][0] = 0.005;
4311  fC[9][1] = 0.005;
4312  fC[9][2] = 1.0000;
4313  fC[9][3] = 0.2618;
4314  fC[9][4] = 0.3070;
4315 
4316  fC[10][0] = 0.005;
4317  fC[10][1] = 0.005;
4318  fC[10][2] = 1.0000;
4319  fC[10][3] = 0.3067;
4320  fC[10][4] = 0.3625;
4321 
4322  fC[11][0] = 0.005;
4323  fC[11][1] = 0.005;
4324  fC[11][2] = 1.0000;
4325  fC[11][3] = 0.3336;
4326  fC[11][4] = 0.4188;
4327 
4328  fC[12][0] = 0.005;
4329  fC[12][1] = 0.005;
4330  fC[12][2] = 1.0000;
4331  fC[12][3] = 0.3706;
4332  fC[12][4] = 0.4511;
4333 
4334  fC[13][0] = 0.005;
4335  fC[13][1] = 0.005;
4336  fC[13][2] = 1.0000;
4337  fC[13][3] = 0.3765;
4338  fC[13][4] = 0.4729;
4339 
4340  fC[14][0] = 0.005;
4341  fC[14][1] = 0.005;
4342  fC[14][2] = 1.0000;
4343  fC[14][3] = 0.3942;
4344  fC[14][4] = 0.4855;
4345 
4346  fC[15][0] = 0.005;
4347  fC[15][1] = 0.005;
4348  fC[15][2] = 1.0000;
4349  fC[15][3] = 0.4051;
4350  fC[15][4] = 0.4762;
4351 
4352  fC[16][0] = 0.005;
4353  fC[16][1] = 0.005;
4354  fC[16][2] = 1.0000;
4355  fC[16][3] = 0.3843;
4356  fC[16][4] = 0.4763;
4357 
4358  fC[17][0] = 0.005;
4359  fC[17][1] = 0.005;
4360  fC[17][2] = 1.0000;
4361  fC[17][3] = 0.4237;
4362  fC[17][4] = 0.4773;
4363 
4364  for(Int_t i=18;i<fgkPIDptBin;i++){
4365  fC[i][0] = fC[17][0];
4366  fC[i][1] = fC[17][1];
4367  fC[i][2] = fC[17][2];
4368  fC[i][3] = fC[17][3];
4369  fC[i][4] = fC[17][4];
4370  }
4371  }
4372  // 60-70%
4373  else if(centrCur < 70){
4374  fC[0][0] = 0.005;
4375  fC[0][1] = 0.005;
4376  fC[0][2] = 1.0000;
4377  fC[0][3] = 0.0010;
4378  fC[0][4] = 0.0010;
4379 
4380  fC[1][0] = 0.005;
4381  fC[1][1] = 0.005;
4382  fC[1][2] = 1.0000;
4383  fC[1][3] = 0.0071;
4384  fC[1][4] = 0.0012;
4385 
4386  fC[2][0] = 0.005;
4387  fC[2][1] = 0.005;
4388  fC[2][2] = 1.0000;
4389  fC[2][3] = 0.0336;
4390  fC[2][4] = 0.0097;
4391 
4392  fC[3][0] = 0.005;
4393  fC[3][1] = 0.005;
4394  fC[3][2] = 1.0000;
4395  fC[3][3] = 0.0662;
4396  fC[3][4] = 0.0460;
4397 
4398  fC[4][0] = 0.005;
4399  fC[4][1] = 0.005;
4400  fC[4][2] = 1.0000;
4401  fC[4][3] = 0.0954;
4402  fC[4][4] = 0.0902;
4403 
4404  fC[5][0] = 0.005;
4405  fC[5][1] = 0.005;
4406  fC[5][2] = 1.0000;
4407  fC[5][3] = 0.1181;
4408  fC[5][4] = 0.1306;
4409 
4410  fC[6][0] = 0.005;
4411  fC[6][1] = 0.005;
4412  fC[6][2] = 1.0000;
4413  fC[6][3] = 0.1481;
4414  fC[6][4] = 0.1662;
4415 
4416  fC[7][0] = 0.005;
4417  fC[7][1] = 0.005;
4418  fC[7][2] = 1.0000;
4419  fC[7][3] = 0.1765;
4420  fC[7][4] = 0.1963;
4421 
4422  fC[8][0] = 0.005;
4423  fC[8][1] = 0.005;
4424  fC[8][2] = 1.0000;
4425  fC[8][3] = 0.2155;
4426  fC[8][4] = 0.2433;
4427 
4428  fC[9][0] = 0.005;
4429  fC[9][1] = 0.005;
4430  fC[9][2] = 1.0000;
4431  fC[9][3] = 0.2580;
4432  fC[9][4] = 0.3022;
4433 
4434  fC[10][0] = 0.005;
4435  fC[10][1] = 0.005;
4436  fC[10][2] = 1.0000;
4437  fC[10][3] = 0.2872;
4438  fC[10][4] = 0.3481;
4439 
4440  fC[11][0] = 0.005;
4441  fC[11][1] = 0.005;
4442  fC[11][2] = 1.0000;
4443  fC[11][3] = 0.3170;
4444  fC[11][4] = 0.3847;
4445 
4446  fC[12][0] = 0.005;
4447  fC[12][1] = 0.005;
4448  fC[12][2] = 1.0000;
4449  fC[12][3] = 0.3454;
4450  fC[12][4] = 0.4258;
4451 
4452  fC[13][0] = 0.005;
4453  fC[13][1] = 0.005;
4454  fC[13][2] = 1.0000;
4455  fC[13][3] = 0.3580;
4456  fC[13][4] = 0.4299;
4457 
4458  fC[14][0] = 0.005;
4459  fC[14][1] = 0.005;
4460  fC[14][2] = 1.0000;
4461  fC[14][3] = 0.3903;
4462  fC[14][4] = 0.4326;
4463 
4464  fC[15][0] = 0.005;
4465  fC[15][1] = 0.005;
4466  fC[15][2] = 1.0000;
4467  fC[15][3] = 0.3690;
4468  fC[15][4] = 0.4491;
4469 
4470  fC[16][0] = 0.005;
4471