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