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