PWG/FLOW/Base/AliFlowAnalysisWithMixedHarmonics.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 /* $Id$ */
  17
  18 /**********************************************************
  19  * In this class azimuthal correlators in mixed harmonics *
  20  * are implemented in terms of Q-vectors. This approach   *
  21  * doesn't require evaluation of nested loops. This class *
  22  * can be used to:                                        *
  23  *                                                        *
  24  *  a) Extract subdominant harmonics (like v1 and v4);    *
  25  *  b) Study flow of two-particle resonances;             *
  26  *  c) Study strong parity violation.                     *
  27  *                                                        *
  28  * Author: Ante Bilandzic (abilandzic@gmail.com)          *
  29  *********************************************************/
  30
  31 #define AliFlowAnalysisWithMixedHarmonics_cxx
  32
  33 #include "Riostream.h"
  34 #include "AliFlowCommonConstants.h"
  35 #include "AliFlowCommonHist.h"
  36 #include "AliFlowCommonHistResults.h"
  37
  38 #include "TMath.h"
  39 #include "TFile.h"
  40 #include "TList.h"
  41 #include "TProfile.h"
  42 #include "TProfile2D.h"
  43
  44 #include "AliFlowEventSimple.h"
  45 #include "AliFlowTrackSimple.h"
  46 #include "AliFlowAnalysisWithMixedHarmonics.h"
  47
  48 class TH1;
  49 class TList;
  50
  51 using std::endl;
  52 using std::cout;
  53 ClassImp(AliFlowAnalysisWithMixedHarmonics)
  54
  55 //================================================================================================================
  56 AliFlowAnalysisWithMixedHarmonics::AliFlowAnalysisWithMixedHarmonics():
  57 fHistList(NULL),
  58 fHistListName(NULL),
  59 fHarmonic(1),
  60 fAnalysisLabel(NULL),
  61 fAnalysisSettings(NULL),
  62 fNoOfMultipicityBins(100),
  63 fMultipicityBinWidth(1),
  64 fMinMultiplicity(3),
  65 fOppositeChargesPOI(kFALSE),
  66 fEvaluateDifferential3pCorrelator(kFALSE),
  67 fCorrectForDetectorEffects(kFALSE),
  68 fPrintOnTheScreen(kTRUE),
  69 fCalculateVsM(kFALSE),
  70 fShowBinLabelsVsM(kFALSE),
  71 fCommonHists(NULL),
  72 fnBinsPhi(0),
  73 fPhiMin(0),
  74 fPhiMax(0),
  75 fPhiBinWidth(0),
  76 fnBinsPt(0),
  77 fPtMin(0),
  78 fPtMax(0),
  79 fPtBinWidth(0),
  80 fnBinsEta(0),
  81 fEtaMin(0),
  82 fEtaMax(0),
  83 fEtaBinWidth(0),
  84 fCommonConstants(NULL),
  85 fWeightsList(NULL),
  86 fUsePhiWeights(kFALSE),
  87 fUsePtWeights(kFALSE),
  88 fUseEtaWeights(kFALSE),
  89 fUseParticleWeights(NULL),
  90 fPhiWeights(NULL),
  91 fPtWeights(NULL),
  92 fEtaWeights(NULL),
  93 fReQnk(NULL),
  94 fImQnk(NULL),
  95 fSpk(NULL),
  96 fProfileList(NULL),
  97 f3pCorrelatorPro(NULL),
  98 f5pCorrelatorPro(NULL),
  99 fNonIsotropicTermsPro(NULL),
 100 f3pCorrelatorVsMPro(NULL),
 101 f3pPOICorrelatorVsM(NULL),
 102 fNonIsotropicTermsVsMPro(NULL),
 103 fNonIsotropicTermsList(NULL),
 104 f2pCorrelatorCosPsiDiffPtDiff(NULL),
 105 f2pCorrelatorCosPsiSumPtDiff(NULL),
 106 f2pCorrelatorSinPsiDiffPtDiff(NULL),
 107 f2pCorrelatorSinPsiSumPtDiff(NULL),
 108 f2pCorrelatorCosPsiDiffPtSum(NULL),
 109 f2pCorrelatorCosPsiSumPtSum(NULL),
 110 f2pCorrelatorSinPsiDiffPtSum(NULL),
 111 f2pCorrelatorSinPsiSumPtSum(NULL),
 112 f2pCorrelatorCosPsiDiffEtaDiff(NULL),
 113 f2pCorrelatorCosPsiSumEtaDiff(NULL),
 114 f2pCorrelatorSinPsiDiffEtaDiff(NULL),
 115 f2pCorrelatorSinPsiSumEtaDiff(NULL),
 116 f2pCorrelatorCosPsiDiffEtaSum(NULL),
 117 f2pCorrelatorCosPsiSumEtaSum(NULL),
 118 f2pCorrelatorSinPsiDiffEtaSum(NULL),
 119 f2pCorrelatorSinPsiSumEtaSum(NULL),
 120 fResultsList(NULL),
 121 f3pCorrelatorHist(NULL),
 122 fDetectorBiasHist(NULL),
 123 f3pCorrelatorVsMHist(NULL),
 124 fDetectorBiasVsMHist(NULL)
 125 {
 126  // Constructor.
 127
 128  // Base list to hold all output objects:
 129  fHistList = new TList();
 130  fHistListName = new TString("cobjMH");
 131  fHistList->SetName(fHistListName->Data());
 132  fHistList->SetOwner(kTRUE);
 133
 134  // List to hold histograms with phi, pt and eta weights:
 135  fWeightsList = new TList();
 136
 137  // List to hold all all-event profiles:
 138  fProfileList = new TList();
 139
 140  // List to hold profiles with all non-isotropic terms for diff. correlators:
 141  fNonIsotropicTermsList = new TList();
 142
 143  // List to hold objects with final results:
 144  fResultsList = new TList();
 145
 146  // Initialize all arrays:
 147  this->InitializeArrays();
 148
 149 } // AliFlowAnalysisWithMixedHarmonics::AliFlowAnalysisWithMixedHarmonics()
 150
 151 //================================================================================================================
 152
 153 AliFlowAnalysisWithMixedHarmonics::~AliFlowAnalysisWithMixedHarmonics()
 154 {
 155  // Destructor.
 156
 157  delete fHistList;
 158
 159 } // end of AliFlowAnalysisWithMixedHarmonics::~AliFlowAnalysisWithMixedHarmonics()
 160
 161 //================================================================================================================
 162
 163 void AliFlowAnalysisWithMixedHarmonics::Init()
 164 {
 165  // Initialize and book all objects.
 166
 167  // a) Cross check if the user settings make sense before starting;
 168  // b) Access all common constants;
 169  // c) Book and nest all lists in the base list fHistList;
 170  // d) Book common control histograms;
 171  // e) Book all event-by-event quantities;
 172  // f) Book all all-event quantities;
 173  // g) Book and fill histograms to hold phi, pt and eta weights;
 174  // h) Store harmonic n used in cos[n*(phi1+phi2-2phi3)] and cos[n*(psi1+psi2-2phi3)].
 175
 176  //save old value and prevent histograms from being added to directory
 177  //to avoid name clashes in case multiple analaysis objects are used
 178  //in an analysis
 179  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 180  TH1::AddDirectory(kFALSE);
 181
 182  TH1::SetDefaultSumw2();
 183
 184  this->CrossCheckSettings();
 185  this->AccessConstants("Init");
 186  this->BookAndNestAllLists();
 187  this->BookProfileHoldingSettings();
 188  this->BookCommonHistograms();
 189  this->BookAllEventByEventQuantities();
 190  this->BookAllAllEventQuantities();
 191  this->BookAndFillWeightsHistograms();
 192  this->StoreHarmonic();
 193
 194  TH1::AddDirectory(oldHistAddStatus);
 195
 196 } // end of void AliFlowAnalysisWithMixedHarmonics::Init()
 197
 198 //================================================================================================================
 199
 200 void AliFlowAnalysisWithMixedHarmonics::Make(AliFlowEventSimple* anEvent)
 201 {
 202  // Running over data only in this method.
 203
 204  // a) Check all pointers used in this method;
 205  // b) Define local variables;
 206  // c) Fill common control histograms;
 207  // d) Loop over data and calculate e-b-e quantities Q_{n,k} and S_{p,k};
 208  // e) Calculate 3-p azimuthal correlator cos[n(phi1+phi2-2*phi3)] and non-isotropic terms in terms of Q_{n,k} and S_{p,k};
 209  // f) Calculate differential 3-p azimuthal correlator cos[n(psi1+psi2-2*phi3)] in terms of Q_{2n} and p_{n}:
 210  // g) Reset all event-by-event quantities.
 211
 212  // a) Check all pointers used in this method:
 213  this->CheckPointersUsedInMake();
 214
 215  // b) Define local variables:
 216  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 217  Double_t dPt  = 0.; // transverse momentum
 218  Double_t dEta = 0.; // pseudorapidity
 219  Double_t wPhi = 1.; // phi weight
 220  Double_t wPt  = 1.; // pt weight
 221  Double_t wEta = 1.; // eta weight
 222  AliFlowTrackSimple *aftsTrack = NULL; // simple track
 223
 224  // c) Fill common control histograms:
 225  fCommonHists->FillControlHistograms(anEvent);
 226
 227  // d) Loop over data and calculate e-b-e quantities:
 228  Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI
 229                                            // nRP   = # of particles used to determine the reaction plane ("Reference Particles");
 230                                            // nPOI  = # of particles of interest for a detailed flow analysis ("Particles of Interest");
 231
 232  Int_t nRefMult = anEvent->GetReferenceMultiplicity();
 233
 234  // Start loop over data:
 235  for(Int_t i=0;i<nPrim;i++)
 236  {
 237   aftsTrack=anEvent->GetTrack(i);
 238   if(aftsTrack)
 239   {
 240    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue; // consider only tracks which are either RPs or POIs
 241    Int_t n = fHarmonic;
 242    if(aftsTrack->InRPSelection()) // checking RP condition:
 243    {
 244     dPhi = aftsTrack->Phi();
 245     dPt  = aftsTrack->Pt();
 246     dEta = aftsTrack->Eta();
 247     if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi-weight for this particle:
 248     {
 249      wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 250     }
 251     if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt-weight for this particle:
 252     {
 253      wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 254     }
 255     if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta-weight for this particle:
 256     {
 257      wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 258     }
 259     // Calculate Re[Q_{m,k}] and Im[Q_{m,k}], (m = 1,2,3,4,5,6 and k = 0,1,2,3) for this event:
 260     for(Int_t m=0;m<6;m++)
 261     {
 262      for(Int_t k=0;k<4;k++) // to be improved (what is the maximum k that I need?)
 263      {
 264       (*fReQnk)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1)*n*dPhi);
 265       (*fImQnk)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1)*n*dPhi);
 266      }
 267     }
 268     // Calculate partially S_{p,k} for this event (final calculation of S_{p,k} follows after the loop over data bellow):
 269     for(Int_t p=0;p<4;p++) // to be improved (what is maximum p that I need?)
 270     {
 271      for(Int_t k=0;k<4;k++) // to be improved (what is maximum k that I need?)
 272      {
 273       (*fSpk)(p,k)+=pow(wPhi*wPt*wEta,k);
 274      }
 275     }
 276    } // end of if(aftsTrack->InRPSelection())
 277    // POIs:
 278    if(fEvaluateDifferential3pCorrelator)
 279    {
 280     if(aftsTrack->InPOISelection()) // 1st POI
 281     {
 282      Double_t dPsi1 = aftsTrack->Phi();
 283      Double_t dPt1 = aftsTrack->Pt();
 284      Double_t dEta1 = aftsTrack->Eta();
 285      Int_t iCharge1 = aftsTrack->Charge();
 286      Bool_t b1stPOIisAlsoRP = kFALSE;
 287      if(aftsTrack->InRPSelection()){b1stPOIisAlsoRP = kTRUE;}
 288      for(Int_t j=0;j<nPrim;j++)
 289      {
 290       if(j==i){continue;}
 291       aftsTrack=anEvent->GetTrack(j);
 292       if(aftsTrack->InPOISelection()) // 2nd POI
 293       {
 294        Double_t dPsi2 = aftsTrack->Phi();
 295        Double_t dPt2 = aftsTrack->Pt();
 296        Double_t dEta2 = aftsTrack->Eta();
 297        Int_t iCharge2 = aftsTrack->Charge();
 298        if(fOppositeChargesPOI && iCharge1 == iCharge2){continue;}
 299        Bool_t b2ndPOIisAlsoRP = kFALSE;
 300        if(aftsTrack->InRPSelection()){b2ndPOIisAlsoRP = kTRUE;}
 301
 302        // Fill:Pt
 303        fRePEBE[0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1+dPsi2)),1.);
 304        fImPEBE[0]->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi1+dPsi2)),1.);
 305        fRePEBE[1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1+dPsi2)),1.);
 306        fImPEBE[1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1+dPsi2)),1.);
 307
 308        // Fill:Eta
 309        fReEtaEBE[0]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1+dPsi2)),1.);
 310        fImEtaEBE[0]->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi1+dPsi2)),1.);
 311        fReEtaEBE[1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1+dPsi2)),1.);
 312        fImEtaEBE[1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi1+dPsi2)),1.);
 313
 314        //=========================================================//
 315        //2particle correlator <cos(n*(psi1 - ps12))> vs |Pt1-Pt2|
 316        f2pCorrelatorCosPsiDiffPtDiff->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1-dPsi2)));
 317        f2pCorrelatorCosPsiSumPtDiff->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1+dPsi2)));
 318        f2pCorrelatorSinPsiDiffPtDiff->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1-dPsi2)));
 319        f2pCorrelatorSinPsiSumPtDiff->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1+dPsi2)));
 320        //_______________________________________________________//
 321        //2particle correlator <cos(n*(psi1 - ps12))> vs (Pt1+Pt2)/2
 322        f2pCorrelatorCosPsiDiffPtSum->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1-dPsi2)));
 323        f2pCorrelatorCosPsiSumPtSum->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1+dPsi2)));
 324        f2pCorrelatorSinPsiDiffPtSum->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi1-dPsi2)));
 325        f2pCorrelatorSinPsiSumPtSum->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi1+dPsi2)));
 326        //_______________________________________________________//
 327        //2particle correlator <cos(n*(psi1 - ps12))> vs |eta1-eta2|
 328        f2pCorrelatorCosPsiDiffEtaDiff->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1-dPsi2)));
 329        f2pCorrelatorCosPsiSumEtaDiff->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1+dPsi2)));
 330        f2pCorrelatorSinPsiDiffEtaDiff->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi1-dPsi2)));
 331        f2pCorrelatorSinPsiSumEtaDiff->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi1+dPsi2)));
 332        //_______________________________________________________//
 333        //2particle correlator <cos(n*(psi1 - ps12))> vs (Pt1+Pt2)/2
 334        f2pCorrelatorCosPsiDiffEtaSum->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1-dPsi2)));
 335        f2pCorrelatorCosPsiSumEtaSum->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1+dPsi2)));
 336        f2pCorrelatorSinPsiDiffEtaSum->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi1-dPsi2)));
 337        f2pCorrelatorSinPsiSumEtaSum->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi1+dPsi2)));
 338        //=========================================================//
 339
 340        // non-isotropic terms, 1st POI:
 341        fReNITEBE[0][0][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1)),1.);
 342        fReNITEBE[0][0][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1)),1.);
 343        fReNITEBE[0][0][2]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1)),1.);
 344        fReNITEBE[0][0][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1)),1.);
 345        fImNITEBE[0][0][0]->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi1)),1.);
 346        fImNITEBE[0][0][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1)),1.);
 347        fImNITEBE[0][0][2]->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi1)),1.);
 348        fImNITEBE[0][0][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi1)),1.);
 349        // non-isotropic terms, 2nd POI:
 350        fReNITEBE[1][0][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi2)),1.);
 351        fReNITEBE[1][0][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi2)),1.);
 352        fReNITEBE[1][0][2]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi2)),1.);
 353        fReNITEBE[1][0][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi2)),1.);
 354        fImNITEBE[1][0][0]->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi2)),1.);
 355        fImNITEBE[1][0][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi2)),1.);
 356        fImNITEBE[1][0][2]->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi2)),1.);
 357        fImNITEBE[1][0][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi2)),1.);
 358
 359        if(b1stPOIisAlsoRP)
 360        {
 361         fOverlapEBE[0][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 362         fOverlapEBE[0][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 363         fOverlapEBE2[0][0]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 364         fOverlapEBE2[0][1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 365         // non-isotropic terms, 1st POI:
 366         fReNITEBE[0][1][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1)),1.);
 367         fReNITEBE[0][1][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1)),1.);
 368         fReNITEBE[0][1][2]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1)),1.);
 369         fReNITEBE[0][1][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1)),1.);
 370         fImNITEBE[0][1][0]->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi1)),1.);
 371         fImNITEBE[0][1][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1)),1.);
 372         fImNITEBE[0][1][2]->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi1)),1.);
 373         fImNITEBE[0][1][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi1)),1.);
 374        }
 375        if(b2ndPOIisAlsoRP)
 376        {
 377         fOverlapEBE[1][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 378         fOverlapEBE[1][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 379         fOverlapEBE2[1][0]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 380         fOverlapEBE2[1][1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 381         // non-isotropic terms, 2nd POI:
 382         fReNITEBE[1][1][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi2)),1.);
 383         fReNITEBE[1][1][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi2)),1.);
 384         fReNITEBE[1][1][2]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi2)),1.);
 385         fReNITEBE[1][1][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi2)),1.);
 386         fImNITEBE[1][1][0]->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi2)),1.);
 387         fImNITEBE[1][1][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi2)),1.);
 388         fImNITEBE[1][1][2]->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi2)),1.);
 389         fImNITEBE[1][1][3]->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi2)),1.);
 390        }
 391       } // end of if(aftsTrack->InPOISelection()) // 2nd POI
 392      } // end of for(Int_t j=i+1;j<nPrim;j++)
 393     } // end of if(aftsTrack->InPOISelection()) // 1st POI
 394    } // end of if(fEvaluateDifferential3pCorrelator)
 395   } else // to if(aftsTrack)
 396     {
 397      cout<<endl;
 398      cout<<" WARNING (MH): No particle! (i.e. aftsTrack is a NULL pointer in Make().)"<<endl;
 399      cout<<endl;
 400     }
 401  } // end of for(Int_t i=0;i<nPrim;i++)
 402
 403  // Calculate the final expressions for S_{p,k}:
 404  for(Int_t p=0;p<4;p++) // to be improved (what is maximum p that I need?)
 405  {
 406   for(Int_t k=0;k<4;k++) // to be improved (what is maximum k that I need?)
 407   {
 408    (*fSpk)(p,k)=pow((*fSpk)(p,k),p+1);
 409   }
 410  }
 411
 412  // e) Calculate 3-p correlator cos[n(phi1+phi2-2*phi3)] in terms of Q_{n,k} and S_{p,k}:
 413  if(anEvent->GetEventNSelTracksRP() >= 3)
 414  {
 415   this->Calculate3pCorrelator();
 416   this->CalculateNonIsotropicTerms();
 417   if(anEvent->GetEventNSelTracksRP() >= 5)
 418   {
 419    this->Calculate5pCorrelator();
 420   } // end of if(anEvent->GetEventNSelTracksRP() >= 5)
 421  } // end of if(anEvent->GetEventNSelTracksRP() >= 3)
 422
 423  // f) Calculate differential 3-p azimuthal correlator cos[n(psi1+psi2-2*phi3)] in terms of Q_{2n} and p_{n}:
 424  if(fEvaluateDifferential3pCorrelator && anEvent->GetEventNSelTracksRP() >= 1)
 425  {
 426    Double_t gIntegrated3pCorrelator = 0.;
 427    this->CalculateDifferential3pCorrelator(gIntegrated3pCorrelator); // to be improved - add relevant if statements for the min # POIs as well
 428
 429   //3particle correlator vs ref. mult
 430    if(fCalculateVsM)
 431      f3pPOICorrelatorVsM->Fill(nRefMult,gIntegrated3pCorrelator);
 432  }
 433
 434  // g) Reset all event-by-event quantities:
 435  this->ResetEventByEventQuantities();
 436
 437 } // end of AliFlowAnalysisWithMixedHarmonics::Make(AliFlowEventSimple* anEvent)
 438
 439 //================================================================================================================
 440
 441 void AliFlowAnalysisWithMixedHarmonics::Finish()
 442 {
 443  // Calculate the final results.
 444
 445  // a) Check all pointers used in this method;
 446  // b) Access common constants;
 447  // c) Access settings for analysis with mixed harmonics;
 448  // d) Correct for detector effects;
 449  // e) Print on the screen the final results.
 450
 451  this->CheckPointersUsedInFinish();
 452  this->AccessConstants("Finish");
 453  this->AccessSettings();
 454  this->CorrectForDetectorEffects();
 455  if(fCalculateVsM){this->CorrectForDetectorEffectsVsM();}
 456  if(fPrintOnTheScreen){this->PrintOnTheScreen();}
 457
 458 } // end of AliFlowAnalysisWithMixedHarmonics::Finish()
 459
 460 //================================================================================================================
 461
 462 void AliFlowAnalysisWithMixedHarmonics::GetOutputHistograms(TList *outputListHistos)
 463 {
 464  // Get pointers to all objects saved in the output file.
 465
 466  // a) Get pointers for common control histograms.
 467  if(outputListHistos)
 468  {
 469   this->SetHistList(outputListHistos);
 470   if(!fHistList)
 471   {
 472    cout<<endl;
 473    cout<<" WARNING (MH): fHistList is NULL in GetOutputHistograms() !!!!"<<endl;
 474    cout<<endl;
 475    exit(0);
 476   }
 477   this->GetPointersForBaseHistograms();
 478   this->GetPointersForCommonHistograms();
 479   this->GetPointersForAllEventProfiles();
 480   this->GetPointersForResultsHistograms();
 481  } else
 482    {
 483     cout<<endl;
 484     cout<<" WARNING (MH): outputListHistos is NULL in GetOutputHistograms() !!!!"<<endl;
 485     cout<<endl;
 486     exit(0);
 487    }
 488
 489 } // end of void AliFlowAnalysisWithMixedHarmonics::GetOutputHistograms(TList *outputListHistos)
 490
 491 //================================================================================================================
 492
 493 void AliFlowAnalysisWithMixedHarmonics::GetPointersForBaseHistograms()
 494 {
 495  // Get pointers to base histograms.
 496
 497  TString analysisSettingsName = "fAnalysisSettings";
 498  TProfile *analysisSettings = dynamic_cast<TProfile*>(fHistList->FindObject(analysisSettingsName.Data()));
 499  if(analysisSettings)
 500  {
 501   this->SetAnalysisSettings(analysisSettings);
 502  } else
 503    {
 504     cout<<endl;
 505     cout<<" WARNING (MH): analysisSettings is NULL in GetPointersForBaseHistograms() !!!!"<<endl;
 506     cout<<endl;
 507     exit(0);
 508    }
 509
 510  TString sCommonConstantsName = "fCommonConstants";
 511  fCommonConstants = dynamic_cast<TProfile*>(fHistList->FindObject(sCommonConstantsName.Data()));
 512  if(!fCommonConstants)
 513  {
 514   printf("\n WARNING (MH): fCommonConstants is NULL in GetPointersForBaseHistograms() !!!!\n\n");
 515   exit(0);
 516  }
 517
 518 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForBaseHistograms()
 519
 520 //================================================================================================================
 521
 522 void AliFlowAnalysisWithMixedHarmonics::GetPointersForCommonHistograms()
 523 {
 524  // Get pointers to common control histograms.
 525
 526  TString commonHistsName = "AliFlowCommonHistMH";
 527  AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(fHistList->FindObject(commonHistsName.Data()));
 528  if(commonHist)
 529  {
 530   this->SetCommonHists(commonHist);
 531  } else
 532    {
 533     cout<<endl;
 534     cout<<" WARNING (MH): commonHist is NULL in GetPointersForCommonHistograms() !!!!"<<endl;
 535     cout<<endl;
 536     exit(0);
 537    }
 538
 539 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForCommonHistograms()
 540
 541 //================================================================================================================
 542
 543 void AliFlowAnalysisWithMixedHarmonics::GetPointersForAllEventProfiles()
 544 {
 545  // Get pointers to profiles holding final results.
 546
 547  TList *profileList = NULL;
 548  profileList = dynamic_cast<TList*>(fHistList->FindObject("Profiles"));
 549  if(!profileList)
 550  {
 551   cout<<endl;
 552   cout<<" WARNING (MH): profileList is NULL in GetPointersForAllEventProfiles() !!!!"<<endl;
 553   cout<<endl;
 554   exit(0);
 555  }
 556  TList *nonIsotropicTermsList = NULL;
 557  nonIsotropicTermsList = dynamic_cast<TList*>(profileList->FindObject("Nonisotropic Terms"));
 558  if(!nonIsotropicTermsList && fEvaluateDifferential3pCorrelator)
 559  {
 560   cout<<endl;
 561   cout<<" WARNING (MH): nonIsotropicTerms is NULL in GetPointersForAllEventProfiles() !!!!"<<endl;
 562   cout<<endl;
 563   exit(0);
 564  }
 565
 566  TString s3pCorrelatorProName = "f3pCorrelatorPro";
 567  TProfile *p3pCorrelatorPro = dynamic_cast<TProfile*>(profileList->FindObject(s3pCorrelatorProName.Data()));
 568  if(p3pCorrelatorPro)
 569  {
 570   this->Set3pCorrelatorPro(p3pCorrelatorPro);
 571  }
 572  TString s3pCorrelatorVsMProName = "f3pCorrelatorVsMPro";
 573  TProfile *p3pCorrelatorVsMPro = dynamic_cast<TProfile*>(profileList->FindObject(s3pCorrelatorVsMProName.Data()));
 574  if(p3pCorrelatorVsMPro)
 575  {
 576   this->Set3pCorrelatorVsMPro(p3pCorrelatorVsMPro);
 577  }
 578  TString s3pPOICorrelatorVsMName = "f3pPOICorrelatorVsM";
 579  TProfile *p3pPOICorrelatorVsM = dynamic_cast<TProfile*>(profileList->FindObject(s3pPOICorrelatorVsMName.Data()));
 580  if(p3pPOICorrelatorVsM)
 581  {
 582   this->Set3pPOICorrelatorVsM(p3pPOICorrelatorVsM);
 583  }
 584  TString nonIsotropicTermsProName = "fNonIsotropicTermsPro";
 585  TProfile *nonIsotropicTermsPro = dynamic_cast<TProfile*>(profileList->FindObject(nonIsotropicTermsProName.Data()));
 586  if(nonIsotropicTermsPro)
 587  {
 588   this->SetNonIsotropicTermsPro(nonIsotropicTermsPro);
 589  }
 590  TString nonIsotropicTermsVsMProName = "fNonIsotropicTermsVsMPro";
 591  TProfile2D *nonIsotropicTermsVsMPro = dynamic_cast<TProfile2D*>(profileList->FindObject(nonIsotropicTermsVsMProName.Data()));
 592  if(nonIsotropicTermsVsMPro)
 593  {
 594   this->SetNonIsotropicTermsVsMPro(nonIsotropicTermsVsMPro);
 595  }
 596  TString psdFlag[2] = {"PtSum","PtDiff"};
 597  TString psdFlag2[2] = {"EtaSum","EtaDiff"};
 598  TString nonIsotropicTerm[10] = {"#LT#LTcos(#psi_{POI_1})#GT#GT","#LT#LTsin(#psi_{POI_1})#GT#GT",
 599                                  "#LT#LTcos(#psi_{POI_2})#GT#GT","#LT#LTsin(#psi_{POI_2})#GT#GT",
 600                                  "#LT#LTcos(#psi_{POI_1}-2#phi_{RP})#GT#GT","#LT#LTsin(#psi_{POI_1}-2#phi_{RP})#GT#GT",
 601                                  "#LT#LTcos(#psi_{POI_2}-2#phi_{RP})#GT#GT","#LT#LTsin(#psi_{POI_2}-2#phi_{RP})#GT#GT",
 602                                  "#LT#LTcos(#psi_{POI_1}+#psi_{POI_2})#GT#GT","#LT#LTsin(#psi_{POI_1}+#psi_{POI_2})#GT#GT"};
 603  for(Int_t sd=0;sd<2;sd++)
 604  {
 605   TProfile *p3pCorrelatorVsPtSumDiffPro = dynamic_cast<TProfile*>(profileList->FindObject(Form("f3pCorrelatorVs%sPro",psdFlag[sd].Data())));
 606   if(p3pCorrelatorVsPtSumDiffPro)
 607   {
 608    this->Set3pCorrelatorVsPtSumDiffPro(p3pCorrelatorVsPtSumDiffPro,sd);
 609   }
 610   TProfile *p3pCorrelatorVsEtaSumDiffPro = dynamic_cast<TProfile*>(profileList->FindObject(Form("f3pCorrelatorVs%sPro",psdFlag2[sd].Data())));
 611   if(p3pCorrelatorVsEtaSumDiffPro)
 612   {
 613    this->Set3pCorrelatorVsEtaSumDiffPro(p3pCorrelatorVsEtaSumDiffPro,sd);
 614   }
 615   if(nonIsotropicTermsList)
 616   {
 617    for(Int_t t=0;t<10;t++)
 618    {
 619     // Pt:
 620     TProfile *pNonIsotropicTermsVsPtSumDiffPro = dynamic_cast<TProfile*>
 621                       (nonIsotropicTermsList->FindObject(Form("fNonIsotropicTermsVs%sPro %s",psdFlag[sd].Data(),nonIsotropicTerm[t].Data())));
 622     if(pNonIsotropicTermsVsPtSumDiffPro)
 623     {
 624      this->SetNonIsotropicTermsVsPtSumDiffPro(pNonIsotropicTermsVsPtSumDiffPro,sd,t);
 625     }
 626     // Eta:
 627     TProfile *pNonIsotropicTermsVsEtaSumDiffPro = dynamic_cast<TProfile*>
 628                       (nonIsotropicTermsList->FindObject(Form("fNonIsotropicTermsVs%sPro %s",psdFlag2[sd].Data(),nonIsotropicTerm[t].Data())));
 629     if(pNonIsotropicTermsVsEtaSumDiffPro)
 630     {
 631      this->SetNonIsotropicTermsVsEtaSumDiffPro(pNonIsotropicTermsVsEtaSumDiffPro,sd,t);
 632     }
 633    } // end of for(Int_t t=0;t<10;t++)
 634   } // end of if(nonIsotropicTermsList)
 635  } // end of for(Int_t sd=0;sd<2;sd++)
 636
 637  //2p correlator vs |Pt1-Pt2|
 638  TProfile *g2pCorrelatorCosPsiDiffPtDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiDiffPtDiff"));
 639  if(g2pCorrelatorCosPsiDiffPtDiff)
 640    this->Set2pCorrelatorCosPsiDiffPtDiff(g2pCorrelatorCosPsiDiffPtDiff);
 641  TProfile *g2pCorrelatorCosPsiSumPtDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiSumPtDiff"));
 642  if(g2pCorrelatorCosPsiSumPtDiff)
 643    this->Set2pCorrelatorCosPsiSumPtDiff(g2pCorrelatorCosPsiSumPtDiff);
 644  TProfile *g2pCorrelatorSinPsiDiffPtDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiDiffPtDiff"));
 645  if(g2pCorrelatorSinPsiDiffPtDiff)
 646    this->Set2pCorrelatorSinPsiDiffPtDiff(g2pCorrelatorSinPsiDiffPtDiff);
 647  TProfile *g2pCorrelatorSinPsiSumPtDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiSumPtDiff"));
 648  if(g2pCorrelatorSinPsiSumPtDiff)
 649    this->Set2pCorrelatorSinPsiSumPtDiff(g2pCorrelatorSinPsiSumPtDiff);
 650
 651  //2p correlator vs (Pt1+Pt2)/2
 652  TProfile *g2pCorrelatorCosPsiDiffPtSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiDiffPtSum"));
 653  if(g2pCorrelatorCosPsiDiffPtSum)
 654    this->Set2pCorrelatorCosPsiDiffPtSum(g2pCorrelatorCosPsiDiffPtSum);
 655  TProfile *g2pCorrelatorCosPsiSumPtSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiSumPtSum"));
 656  if(g2pCorrelatorCosPsiSumPtSum)
 657    this->Set2pCorrelatorCosPsiSumPtSum(g2pCorrelatorCosPsiSumPtSum);
 658  TProfile *g2pCorrelatorSinPsiDiffPtSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiDiffPtSum"));
 659  if(g2pCorrelatorSinPsiDiffPtSum)
 660    this->Set2pCorrelatorSinPsiDiffPtSum(g2pCorrelatorSinPsiDiffPtSum);
 661  TProfile *g2pCorrelatorSinPsiSumPtSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiSumPtSum"));
 662  if(g2pCorrelatorSinPsiSumPtSum)
 663    this->Set2pCorrelatorSinPsiSumPtSum(g2pCorrelatorSinPsiSumPtSum);
 664
 665  //2p correlator vs |eta1-eta2|
 666  TProfile *g2pCorrelatorCosPsiDiffEtaDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiDiffEtaDiff"));
 667  if(g2pCorrelatorCosPsiDiffEtaDiff)
 668    this->Set2pCorrelatorCosPsiDiffEtaDiff(g2pCorrelatorCosPsiDiffEtaDiff);
 669  TProfile *g2pCorrelatorCosPsiSumEtaDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiSumEtaDiff"));
 670  if(g2pCorrelatorCosPsiSumEtaDiff)
 671    this->Set2pCorrelatorCosPsiSumEtaDiff(g2pCorrelatorCosPsiSumEtaDiff);
 672  TProfile *g2pCorrelatorSinPsiDiffEtaDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiDiffEtaDiff"));
 673  if(g2pCorrelatorSinPsiDiffEtaDiff)
 674    this->Set2pCorrelatorSinPsiDiffEtaDiff(g2pCorrelatorSinPsiDiffEtaDiff);
 675  TProfile *g2pCorrelatorSinPsiSumEtaDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiSumEtaDiff"));
 676  if(g2pCorrelatorSinPsiSumEtaDiff)
 677    this->Set2pCorrelatorSinPsiSumEtaDiff(g2pCorrelatorSinPsiSumEtaDiff);
 678
 679  //2p correlator vs (eta1+eta2)/2
 680  TProfile *g2pCorrelatorCosPsiDiffEtaSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiDiffEtaSum"));
 681  if(g2pCorrelatorCosPsiDiffEtaSum)
 682    this->Set2pCorrelatorCosPsiDiffEtaSum(g2pCorrelatorCosPsiDiffEtaSum);
 683  TProfile *g2pCorrelatorCosPsiSumEtaSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiSumEtaSum"));
 684  if(g2pCorrelatorCosPsiSumEtaSum)
 685    this->Set2pCorrelatorCosPsiSumEtaSum(g2pCorrelatorCosPsiSumEtaSum);
 686  TProfile *g2pCorrelatorSinPsiDiffEtaSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiDiffEtaSum"));
 687  if(g2pCorrelatorSinPsiDiffEtaSum)
 688    this->Set2pCorrelatorSinPsiDiffEtaSum(g2pCorrelatorSinPsiDiffEtaSum);
 689  TProfile *g2pCorrelatorSinPsiSumEtaSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiSumEtaSum"));
 690  if(g2pCorrelatorSinPsiSumEtaSum)
 691    this->Set2pCorrelatorSinPsiSumEtaSum(g2pCorrelatorSinPsiSumEtaSum);
 692
 693  TString s5pCorrelatorProName = "f5pCorrelatorPro";
 694  TProfile *p5pCorrelatorPro = dynamic_cast<TProfile*>(profileList->FindObject(s5pCorrelatorProName.Data()));
 695  if(p5pCorrelatorPro)
 696  {
 697   this->Set5pCorrelatorPro(p5pCorrelatorPro);
 698  }
 699
 700 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForAllEventProfiles()
 701
 702 //================================================================================================================
 703
 704 void AliFlowAnalysisWithMixedHarmonics::GetPointersForResultsHistograms()
 705 {
 706  // Get pointers to histograms holding final results.
 707
 708  TList *resultsList = NULL;
 709  resultsList = dynamic_cast<TList*>(fHistList->FindObject("Results"));
 710  if(!resultsList)
 711  {
 712   cout<<endl;
 713   cout<<" WARNING (MH): resultsList is NULL in GetPointersForResultsHistograms() !!!!"<<endl;
 714   cout<<endl;
 715   exit(0);
 716  }
 717  TString s3pCorrelatorHistName = "f3pCorrelatorHist";
 718  TH1D *h3pCorrelatorHist = dynamic_cast<TH1D*>(resultsList->FindObject(s3pCorrelatorHistName.Data()));
 719  if(h3pCorrelatorHist)
 720  {
 721   this->Set3pCorrelatorHist(h3pCorrelatorHist);
 722  }
 723  TString s3pCorrelatorVsMHistName = "f3pCorrelatorVsMHist";
 724  TH1D *h3pCorrelatorVsMHist = dynamic_cast<TH1D*>(resultsList->FindObject(s3pCorrelatorVsMHistName.Data()));
 725  if(h3pCorrelatorVsMHist)
 726  {
 727   this->Set3pCorrelatorVsMHist(h3pCorrelatorVsMHist);
 728  }
 729  TString detectorBiasHistName = "fDetectorBiasHist";
 730  TH1D *detectorBiasHist = dynamic_cast<TH1D*>(resultsList->FindObject(detectorBiasHistName.Data()));
 731  if(detectorBiasHist)
 732  {
 733   this->SetDetectorBiasHist(detectorBiasHist);
 734  }
 735  TString detectorBiasVsMHistName = "fDetectorBiasVsMHist";
 736  TH1D *detectorBiasVsMHist = dynamic_cast<TH1D*>(resultsList->FindObject(detectorBiasVsMHistName.Data()));
 737  if(detectorBiasVsMHist)
 738  {
 739   this->SetDetectorBiasVsMHist(detectorBiasVsMHist);
 740  }
 741
 742  TString psdFlag[2] = {"PtSum","PtDiff"};
 743  TString psdFlag2[2] = {"EtaSum","EtaDiff"};
 744  for(Int_t sd=0;sd<2;sd++)
 745  {
 746   TH1D *h3pCorrelatorVsPtSumDiffHist = dynamic_cast<TH1D*>(resultsList->FindObject(Form("f3pCorrelatorVs%sHist",psdFlag[sd].Data())));
 747   if(h3pCorrelatorVsPtSumDiffHist)
 748   {
 749    this->Set3pCorrelatorVsPtSumDiffHist(h3pCorrelatorVsPtSumDiffHist,sd);
 750   }
 751   TH1D *h3pCorrelatorVsEtaSumDiffHist = dynamic_cast<TH1D*>(resultsList->FindObject(Form("f3pCorrelatorVs%sHist",psdFlag2[sd].Data())));
 752   if(h3pCorrelatorVsEtaSumDiffHist)
 753   {
 754    this->Set3pCorrelatorVsEtaSumDiffHist(h3pCorrelatorVsEtaSumDiffHist,sd);
 755   }
 756  } // end of for(Int_t sd=0;sd<2;sd++)
 757
 758 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForResultsHistograms()
 759
 760 //================================================================================================================
 761
 762 void AliFlowAnalysisWithMixedHarmonics::WriteHistograms(TString outputFileName)
 763 {
 764  // Store the final results in output .root file.
 765  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 766  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
 767  delete output;
 768 }
 769
 770 //================================================================================================================
 771
 772 void AliFlowAnalysisWithMixedHarmonics::WriteHistograms(TDirectoryFile *outputFileName)
 773 {
 774  // Store the final results in output .root file.
 775  fHistList->SetName("cobjMH");
 776  fHistList->SetOwner(kTRUE);
 777  outputFileName->Add(fHistList);
 778  outputFileName->Write(outputFileName->GetName(),TObject::kSingleKey);
 779 }
 780
 781 //================================================================================================================
 782
 783 void AliFlowAnalysisWithMixedHarmonics::StoreHarmonic()
 784 {
 785  // Store harmonic n used in cos[n*(phi1+phi2-2phi3)] and cos[n*(psi1+psi2-2phi3)].
 786
 787  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
 788
 789 } // end of void AliFlowAnalysisWithMixedHarmonics::StoreHarmonic()
 790
 791 //================================================================================================================
 792
 793 void AliFlowAnalysisWithMixedHarmonics::InitializeArrays()
 794 {
 795  // Initialize arrays.
 796
 797  for(Int_t sd=0;sd<2;sd++)
 798  {
 799   fRePEBE[sd] = NULL;
 800   fImPEBE[sd] = NULL;
 801   fReEtaEBE[sd] = NULL;
 802   fImEtaEBE[sd] = NULL;
 803   f3pCorrelatorVsPtSumDiffPro[sd] = NULL;
 804   f3pCorrelatorVsEtaSumDiffPro[sd] = NULL;
 805   f3pCorrelatorVsPtSumDiffHist[sd] = NULL;
 806   f3pCorrelatorVsEtaSumDiffHist[sd] = NULL;
 807   for(Int_t t=0;t<10;t++) // non-isotropic terms for diff. correlators
 808   {
 809    fNonIsotropicTermsVsPtSumDiffPro[sd][t] = NULL;
 810    fNonIsotropicTermsVsEtaSumDiffPro[sd][t] = NULL;
 811   }
 812  } // end of for(Int_t sd=0;sd<2;sd++)
 813  for(Int_t fs=0;fs<2;fs++) // 1st/2nd POI which is also RP
 814  {
 815   for(Int_t sd=0;sd<2;sd++)
 816   {
 817    fOverlapEBE[fs][sd] = NULL;
 818    fOverlapEBE2[fs][sd] = NULL;
 819   }
 820  } // end of for(Int_t fs=0;fs<2;fs++) // 1st/2nd POI which is also RP
 821  for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
 822  {
 823   for(Int_t ao=0;ao<2;ao++) // all/overlap
 824   {
 825    for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
 826    {
 827     fReNITEBE[p12][ao][pe] = NULL;
 828     fImNITEBE[p12][ao][pe] = NULL;
 829    } // end of for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
 830   } // end of for(Int_t ao=0;ao<2;ao++) // all/overlap
 831  } // end of for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
 832
 833 } // end of AliFlowAnalysisWithMixedHarmonics::InitializeArrays()
 834
 835 //================================================================================================================
 836
 837 void AliFlowAnalysisWithMixedHarmonics::BookAndNestAllLists()
 838 {
 839  // Book and nest all list in base list fHistList.
 840
 841  // Weights:
 842  fWeightsList->SetName("Weights");
 843  fWeightsList->SetOwner(kTRUE);
 844  fHistList->Add(fWeightsList);
 845  // Profiles:
 846  fProfileList->SetName("Profiles");
 847  fProfileList->SetOwner(kTRUE);
 848  fHistList->Add(fProfileList);
 849  // Results:
 850  fResultsList->SetName("Results");
 851  fResultsList->SetOwner(kTRUE);
 852  fHistList->Add(fResultsList);
 853  // Profiles with non-isotropic terms for diff. correlators:
 854  fNonIsotropicTermsList->SetName("Nonisotropic Terms");
 855  fNonIsotropicTermsList->SetOwner(kTRUE);
 856  if(fEvaluateDifferential3pCorrelator){fProfileList->Add(fNonIsotropicTermsList);}
 857
 858 } // end of void AliFlowAnalysisWithMixedHarmonics::BookAndNestAllLists()
 859
 860 //================================================================================================================
 861
 862 void AliFlowAnalysisWithMixedHarmonics::BookProfileHoldingSettings()
 863 {
 864  // Book profile to hold all analysis settings.
 865
 866  TString analysisSettingsName = "fAnalysisSettings";
 867  fAnalysisSettings = new TProfile(analysisSettingsName.Data(),"Settings for analysis with mixed harmonics",10,0,10);
 868  fAnalysisSettings->SetStats(kFALSE);
 869  fAnalysisSettings->GetXaxis()->SetLabelSize(0.03);
 870  fAnalysisSettings->GetXaxis()->SetBinLabel(1,"Corr. for det. effects?");
 871  fAnalysisSettings->Fill(0.5,(Int_t)fCorrectForDetectorEffects);
 872  fAnalysisSettings->GetXaxis()->SetBinLabel(2,"# of mult. bins");
 873  fAnalysisSettings->Fill(1.5,fNoOfMultipicityBins);
 874  fAnalysisSettings->GetXaxis()->SetBinLabel(3,"Width of mult. bins");
 875  fAnalysisSettings->Fill(2.5,fMultipicityBinWidth);
 876  fAnalysisSettings->GetXaxis()->SetBinLabel(4,"Minimal mult.");
 877  fAnalysisSettings->Fill(3.5,fMinMultiplicity);
 878  fAnalysisSettings->GetXaxis()->SetBinLabel(5,"Print on the screen?");
 879  fAnalysisSettings->Fill(4.5,(Int_t)fPrintOnTheScreen);
 880  fAnalysisSettings->GetXaxis()->SetBinLabel(6,"fHarmonic");
 881  fAnalysisSettings->Fill(5.5,(Int_t)fHarmonic);
 882  fAnalysisSettings->GetXaxis()->SetBinLabel(7,"fOppositeChargesPOI");
 883  fAnalysisSettings->Fill(6.5,(Int_t)fOppositeChargesPOI);
 884  fAnalysisSettings->GetXaxis()->SetBinLabel(8,"fEvaluateDifferential3pCorrelator");
 885  fAnalysisSettings->Fill(7.5,(Int_t)fEvaluateDifferential3pCorrelator);
 886  fAnalysisSettings->GetXaxis()->SetBinLabel(9,"fCalculateVsM");
 887  fAnalysisSettings->Fill(8.5,(Int_t)fCalculateVsM);
 888  fAnalysisSettings->GetXaxis()->SetBinLabel(10,"fShowBinLabelsVsM");
 889  fAnalysisSettings->Fill(9.5,(Int_t)fShowBinLabelsVsM);
 890  fHistList->Add(fAnalysisSettings);
 891
 892 } // end of void AliFlowAnalysisWithMixedHarmonics::BookProfileHoldingSettings()
 893
 894 //================================================================================================================
 895
 896 void AliFlowAnalysisWithMixedHarmonics::BookCommonHistograms()
 897 {
 898  // Book common control histograms and common histograms for final results.
 899
 900  TString commonHistsName = "AliFlowCommonHistMH";
 901  fCommonHists = new AliFlowCommonHist(commonHistsName.Data());
 902  fHistList->Add(fCommonHists);
 903
 904 } // end of void AliFlowAnalysisWithMixedHarmonics::BookCommonHistograms()
 905
 906 //================================================================================================================
 907
 908 void AliFlowAnalysisWithMixedHarmonics::BookAllEventByEventQuantities()
 909 {
 910  // Book all event-by-event quantitites.
 911
 912   // Q_{n,k} and S{p,k}:
 913  fReQnk = new TMatrixD(6,9); // to be improved (check bound on k!)
 914  fImQnk = new TMatrixD(6,9); // to be improved (check bound on k!)
 915  fSpk = new TMatrixD(4,4); // to be improved (check bound on p and k!)
 916
 917  // p_n vs [(p1+p2)/2,|p1-p2|]
 918  if(!fEvaluateDifferential3pCorrelator){return;}
 919  TString psdFlag[2] = {"PtSum","PtDiff"};
 920  TString p2sdFlag[2] = {"PtSum","PtDiff"};
 921  TString fsFlag[2] = {"1st","2nd"};
 922  for(Int_t sd=0;sd<2;sd++)
 923  {
 924   fRePEBE[sd] = new TProfile(Form("fRePEBE%s",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
 925   fImPEBE[sd] = new TProfile(Form("fImPEBE%s",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
 926   fReEtaEBE[sd] = new TProfile(Form("fReEtaEBE%s",p2sdFlag[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
 927   fImEtaEBE[sd] = new TProfile(Form("fImEtaEBE%s",p2sdFlag[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
 928  }
 929  for(Int_t fs=0;fs<2;fs++)
 930  {
 931   for(Int_t sd=0;sd<2;sd++)
 932   {
 933    fOverlapEBE[fs][sd] = new TProfile(Form("%s POI, %s",fsFlag[sd].Data(),psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
 934    fOverlapEBE2[fs][sd] = new TProfile(Form("%s POI 2, %s",fsFlag[sd].Data(),p2sdFlag[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
 935   }
 936  }
 937  Int_t nBinsPtEta[4] = {fnBinsPt,fnBinsPt,fnBinsEta,fnBinsEta};
 938  Double_t dPtEtaMin[4] = {0.,0.,fEtaMin,fEtaMin};
 939  Double_t dPtEtaMax[4] = {fPtMax,fPtMax,fEtaMax,fEtaMax};
 940  for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
 941  {
 942   for(Int_t ao=0;ao<2;ao++) // all/overlap
 943   {
 944    for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
 945    {
 946     fReNITEBE[p12][ao][pe] = new TProfile(Form("fReNITEBE%d%d%d",p12,ao,pe),"",nBinsPtEta[pe],dPtEtaMin[pe],dPtEtaMax[pe]);
 947     fImNITEBE[p12][ao][pe] = new TProfile(Form("fImNITEBE%d%d%d",p12,ao,pe),"",nBinsPtEta[pe],dPtEtaMin[pe],dPtEtaMax[pe]);
 948    } // end of for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
 949   } // end of for(Int_t ao=0;ao<2;ao++) // all/overlap
 950  } // end of for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
 951
 952 } // end fo void AliFlowAnalysisWithMixedHarmonics::BookAllEventByEventQuantities()
 953
 954 //================================================================================================================
 955
 956 void AliFlowAnalysisWithMixedHarmonics::BookAllAllEventQuantities()
 957 {
 958  // Book all all-event quantitites.
 959
 960  // a) Book histos and profiles without any binning in multiplicity, pt or eta;
 961  // b) Book quantites with multiplicity binning;
 962  // c) Book quantites with binning in (p1+p2)/2 and |p1-p2|.
 963
 964  this->BookDefault();
 965  if(fCalculateVsM){this->BookVsM();}
 966  if(fEvaluateDifferential3pCorrelator){this->BookDifferential();}
 967
 968 } // end of void AliFlowAnalysisWithMixedHarmonics::BookAllAllEventQuantities()
 969
 970 //================================================================================================================
 971
 972 void AliFlowAnalysisWithMixedHarmonics::BookDefault()
 973 {
 974  // Book histos and profiles without any binning in multiplicity, pt or eta.
 975
 976  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects);
 977  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>>;
 978  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects;
 979  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>>;
 980  // e) 5-p correlator <<cos[n*(2phi1+2phi2+2phi3-3phi4-3phi5)]>> for all events (not corrected for detector effects - not supported yet).
 981
 982  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects);
 983  TString s3pCorrelatorProName = "f3pCorrelatorPro";
 984  f3pCorrelatorPro = new TProfile(s3pCorrelatorProName.Data(),"",1,0,1);
 985  f3pCorrelatorPro->SetStats(kFALSE);
 986  f3pCorrelatorPro->GetXaxis()->SetLabelOffset(0.01);
 987  f3pCorrelatorPro->GetXaxis()->SetLabelSize(0.05);
 988  if(fHarmonic == 1)
 989  {
 990   f3pCorrelatorPro->GetXaxis()->SetBinLabel(1,"#LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
 991  } else
 992    {
 993     f3pCorrelatorPro->GetXaxis()->SetBinLabel(1,Form("#LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
 994    }
 995  fProfileList->Add(f3pCorrelatorPro);
 996
 997  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>>:
 998  TString nonIsotropicTermsProName = "fNonIsotropicTermsPro";
 999  fNonIsotropicTermsPro = new TProfile(nonIsotropicTermsProName.Data(),"",8,0,8);
1000  fNonIsotropicTermsPro->SetStats(kFALSE);
1001  if(fHarmonic == 1)
1002  {
1003   fNonIsotropicTermsPro->SetTitle("Non-isotropic terms in decomposition of #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
1004   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(1,"cos(#phi_{1})");
1005   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(2,"sin(#phi_{1})");
1006   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(3,"cos(2#phi_{1})");
1007   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(4,"sin(2#phi_{1})");
1008   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(5,"cos(#phi_{1}+#phi_{2})");
1009   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(6,"sin(#phi_{1}+#phi_{2})");
1010   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(7,"cos(2#phi_{1}-#phi_{2})");
1011   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(8,"sin(2#phi_{1}-#phi_{2})");
1012   // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(9,"cos(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
1013   // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(10,"sin(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
1014  } else
1015    {
1016     fNonIsotropicTermsPro->SetTitle(Form("Non-isotropic terms in decomposition of #LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
1017     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(1,Form("cos(%d#phi_{1})",fHarmonic));
1018     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(2,Form("sin(%d#phi_{1})",fHarmonic));
1019     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(3,Form("cos(%d#phi_{1})",2*fHarmonic));
1020     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(4,Form("sin(%d#phi_{1})",2*fHarmonic));
1021     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(5,Form("cos[%d(#phi_{1}+#phi_{2})]",fHarmonic));
1022     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(6,Form("sin[%d(#phi_{1}+#phi_{2})]",fHarmonic));
1023     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(7,Form("cos[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
1024     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(8,Form("sin[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
1025     // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(9,Form("cos(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
1026     // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(10,Form("sin(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
1027    }
1028  fProfileList->Add(fNonIsotropicTermsPro);
1029
1030  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects:
1031  TString s3pCorrelatorHistName = "f3pCorrelatorHist";
1032  f3pCorrelatorHist = new TH1D(s3pCorrelatorHistName.Data(),"",1,0,1);
1033  f3pCorrelatorHist->SetStats(kFALSE);
1034  f3pCorrelatorHist->GetXaxis()->SetLabelOffset(0.01);
1035  f3pCorrelatorHist->GetXaxis()->SetLabelSize(0.05);
1036  if(fHarmonic == 1)
1037  {
1038   f3pCorrelatorHist->GetXaxis()->SetBinLabel(1,"#LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
1039  } else
1040    {
1041     f3pCorrelatorHist->GetXaxis()->SetBinLabel(1,Form("#LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
1042    }
1043  fResultsList->Add(f3pCorrelatorHist);
1044
1045  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>>:
1046  TString detectorBiasHistName = "fDetectorBiasHist";
1047  fDetectorBiasHist = new TH1D(detectorBiasHistName.Data(),"Bias coming from detector inefficiences",1,0,1);
1048  fDetectorBiasHist->SetStats(kFALSE);
1049  if(fHarmonic == 1)
1050  {
1051   fDetectorBiasHist->GetXaxis()->SetBinLabel(1,"#frac{corrected}{measured} #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
1052  } else
1053    {
1054     fDetectorBiasHist->GetXaxis()->SetBinLabel(1,Form("#frac{corrected}{measured} #LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
1055    }
1056  fResultsList->Add(fDetectorBiasHist);
1057
1058  // e) 5-p correlator <<cos[n*(2phi1+2phi2+2phi3-3phi4-3phi5)]>> for all events (not corrected for detector effects - not supported yet):
1059  TString s5pCorrelatorProName = "f5pCorrelatorPro";
1060  f5pCorrelatorPro = new TProfile(s5pCorrelatorProName.Data(),"",1,0,1);
1061  f5pCorrelatorPro->SetStats(kFALSE);
1062  f5pCorrelatorPro->GetXaxis()->SetLabelOffset(0.01);
1063  f5pCorrelatorPro->GetXaxis()->SetLabelSize(0.05);
1064  if(fHarmonic == 1)
1065  {
1066   f5pCorrelatorPro->GetXaxis()->SetBinLabel(1,"#LT#LTcos(2#phi_{1}+2#phi_{2}+2#phi_{3}-3#phi_{4}-3#phi_{5})#GT#GT");
1067  } else
1068    {
1069     f5pCorrelatorPro->GetXaxis()->SetBinLabel(1,Form("#LT#LTcos[%i(2#phi_{1}+2#phi_{2}+2#phi_{3}-3#phi_{4}-3#phi_{5})]#GT#GT",fHarmonic));
1070    }
1071  fProfileList->Add(f5pCorrelatorPro);
1072
1073 } // end of void AliFlowAnalysisWithMixedHarmonics::BookDefault()
1074
1075 //================================================================================================================
1076
1077 void AliFlowAnalysisWithMixedHarmonics::BookVsM()
1078 {
1079  // Book histos and profiles holding results vs multiplicity.
1080
1081  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects) vs M;
1082  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>> vs M;
1083  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects vs M;
1084  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> vs M.
1085
1086  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects) vs M:
1087  TString s3pCorrelatorVsMProName = "f3pCorrelatorVsMPro";
1088  f3pCorrelatorVsMPro = new TProfile(s3pCorrelatorVsMProName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
1089  f3pCorrelatorVsMPro->SetStats(kFALSE);
1090  if(fHarmonic == 1)
1091  {
1092   f3pCorrelatorVsMPro->SetTitle("#LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
1093  } else
1094    {
1095     f3pCorrelatorVsMPro->SetTitle(Form("#LT#LTcos[%d(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} M",fHarmonic));
1096    }
1097  if(fShowBinLabelsVsM)
1098  {
1099   f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
1100   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
1101   {
1102    f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
1103   }
1104   f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
1105  } else
1106    {
1107     f3pCorrelatorVsMPro->GetXaxis()->SetTitle("M");
1108    }
1109  fProfileList->Add(f3pCorrelatorVsMPro);
1110
1111  TString s3pPOICorrelatorVsMName = "f3pPOICorrelatorVsM";
1112  f3pPOICorrelatorVsM = new TProfile(s3pPOICorrelatorVsMName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
1113  f3pPOICorrelatorVsM->SetStats(kFALSE);
1114  if(fHarmonic == 1)
1115  {
1116   f3pPOICorrelatorVsM->SetTitle("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
1117  } else
1118    {
1119     f3pPOICorrelatorVsM->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} M",fHarmonic));
1120    }
1121  if(fShowBinLabelsVsM)
1122  {
1123   f3pPOICorrelatorVsM->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
1124   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
1125   {
1126    f3pPOICorrelatorVsM->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
1127   }
1128   f3pPOICorrelatorVsM->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
1129  } else
1130    {
1131     f3pPOICorrelatorVsM->GetXaxis()->SetTitle("M");
1132    }
1133  fProfileList->Add(f3pPOICorrelatorVsM);
1134
1135  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>> vs M:
1136  TString s3pCorrelatorVsMHistName = "f3pCorrelatorVsMHist";
1137  f3pCorrelatorVsMHist = new TH1D(s3pCorrelatorVsMHistName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
1138  f3pCorrelatorVsMHist->SetStats(kFALSE);
1139  if(fHarmonic == 1)
1140  {
1141   f3pCorrelatorVsMHist->SetTitle("cos(#phi_{1}+#phi_{2}-2#phi_{3}) #font[72]{vs} M");
1142  } else
1143    {
1144     f3pCorrelatorVsMHist->SetTitle(Form("cos[%d(#phi_{1}+#phi_{2}-2#phi_{3})] #font[72]{vs} M",fHarmonic));
1145    }
1146  if(fShowBinLabelsVsM)
1147  {
1148   f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
1149   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
1150   {
1151    f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
1152   }
1153   f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
1154  } else
1155    {
1156     f3pCorrelatorVsMHist->GetXaxis()->SetTitle("M");
1157    }
1158  fResultsList->Add(f3pCorrelatorVsMHist);
1159
1160  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects vs M:
1161  TString nonIsotropicTermsVsMProName = "fNonIsotropicTermsVsMPro";
1162  fNonIsotropicTermsVsMPro = new TProfile2D(nonIsotropicTermsVsMProName.Data(),"",8,0,8,fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
1163  fNonIsotropicTermsVsMPro->SetStats(kFALSE);
1164  if(fHarmonic == 1)
1165  {
1166   fNonIsotropicTermsVsMPro->SetTitle("Non-isotropic terms in decomposition of #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
1167   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(1,"cos(#phi_{1})");
1168   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(2,"sin(#phi_{1})");
1169   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(3,"cos(2#phi_{1})");
1170   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(4,"sin(2#phi_{1})");
1171   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(5,"cos(#phi_{1}+#phi_{2})");
1172   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(6,"sin(#phi_{1}+#phi_{2})");
1173   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(7,"cos(2#phi_{1}-#phi_{2})");
1174   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(8,"sin(2#phi_{1}-#phi_{2})");
1175   // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(9,"cos(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
1176   // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(10,"sin(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
1177  } else
1178    {
1179     fNonIsotropicTermsVsMPro->SetTitle(Form("Non-isotropic terms in decomposition of #LT#LTcos[%d(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
1180     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(1,Form("cos(%d#phi_{1})",fHarmonic));
1181     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(2,Form("sin(%d#phi_{1})",fHarmonic));
1182     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(3,Form("cos(%d#phi_{1})",2*fHarmonic));
1183     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(4,Form("sin(%d#phi_{1})",2*fHarmonic));
1184     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(5,Form("cos[%d(#phi_{1}+#phi_{2})]",fHarmonic));
1185     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(6,Form("sin[%d(#phi_{1}+#phi_{2})]",fHarmonic));
1186     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(7,Form("cos[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
1187     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(8,Form("sin[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
1188     // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(9,Form("cos(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
1189     // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(10,Form("sin(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
1190    }
1191  if(fShowBinLabelsVsM)
1192  {
1193   fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
1194   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
1195   {
1196    fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
1197   }
1198   fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
1199  } else
1200    {
1201     fNonIsotropicTermsVsMPro->GetYaxis()->SetTitle("M");
1202    }
1203  fProfileList->Add(fNonIsotropicTermsVsMPro);
1204
1205  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> vs M:
1206  TString detectorBiasVsMHistName = "fDetectorBiasVsMHist";
1207  fDetectorBiasVsMHist = new TH1D(detectorBiasVsMHistName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
1208  fDetectorBiasVsMHist->SetStats(kFALSE);
1209  if(fHarmonic == 1)
1210  {
1211   fDetectorBiasVsMHist->SetTitle("#frac{corrected}{measured} #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
1212  } else
1213    {
1214     fDetectorBiasVsMHist->SetTitle(Form("#frac{corrected}{measured} #LT#LTcos[%d(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} M",fHarmonic));
1215    }
1216  if(fShowBinLabelsVsM)
1217  {
1218   fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
1219   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
1220   {
1221    fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
1222   }
1223   fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
1224  } else
1225    {
1226     fDetectorBiasVsMHist->GetXaxis()->SetTitle("M");
1227    }
1228  fResultsList->Add(fDetectorBiasVsMHist);
1229
1230 } // end of void AliFlowAnalysisWithMixedHarmonics::BookVsM()
1231
1232 //================================================================================================================
1233
1234 void AliFlowAnalysisWithMixedHarmonics::BookDifferential()
1235 {
1236  // Book histos and profiles holding results vs (p1+p2)/2 and |p1-p2|.
1237
1238  TString psdFlag[2] = {"PtSum","PtDiff"};
1239  TString psdTitleFlag[2] = {"(p_{T,1}+ p_{T,2})/2","#left|p_{T,1}- p_{T,2}#right|"};
1240  TString psdFlag2[2] = {"EtaSum","EtaDiff"};
1241  TString psdTitleFlag2[2] = {"(#eta_{1}+ #eta_{2})/2","#left|#eta_{1}- #eta_{2}#right|"};
1242  //TString s3pCorrelatorVsPtSumDiffProName = "f3pCorrelatorVsPtSumDiffPro";
1243  for(Int_t sd=0;sd<2;sd++)
1244  {
1245   f3pCorrelatorVsPtSumDiffPro[sd] = new TProfile(Form("f3pCorrelatorVs%sPro",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
1246   f3pCorrelatorVsPtSumDiffPro[sd]->SetStats(kFALSE);
1247   f3pCorrelatorVsEtaSumDiffPro[sd] = new TProfile(Form("f3pCorrelatorVs%sPro",psdFlag2[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
1248   f3pCorrelatorVsEtaSumDiffPro[sd]->SetStats(kFALSE);
1249   //f3pCorrelatorVsPtSumDiffPro[sd]->SetLabelSize(0.05);
1250   //f3pCorrelatorVsPtSumDiffPro[sd]->SetMarkerStyle(25);
1251   if(fHarmonic == 1)
1252   {
1253    f3pCorrelatorVsPtSumDiffPro[sd]->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} %s",psdTitleFlag[sd].Data()));
1254    f3pCorrelatorVsEtaSumDiffPro[sd]->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} %s",psdTitleFlag2[sd].Data()));
1255   } else
1256     {
1257      f3pCorrelatorVsPtSumDiffPro[sd]->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[sd].Data()));
1258      f3pCorrelatorVsEtaSumDiffPro[sd]->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[sd].Data()));
1259     }
1260   f3pCorrelatorVsPtSumDiffPro[sd]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
1261   fProfileList->Add(f3pCorrelatorVsPtSumDiffPro[sd]);
1262   f3pCorrelatorVsEtaSumDiffPro[sd]->GetXaxis()->SetTitle(psdTitleFlag2[sd].Data());
1263   fProfileList->Add(f3pCorrelatorVsEtaSumDiffPro[sd]);
1264  }
1265
1266  // Corrected for detector effects:
1267  for(Int_t sd=0;sd<2;sd++)
1268  {
1269   f3pCorrelatorVsPtSumDiffHist[sd] = new TH1D(Form("f3pCorrelatorVs%sHist",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
1270   f3pCorrelatorVsPtSumDiffHist[sd]->SetStats(kFALSE);
1271   f3pCorrelatorVsEtaSumDiffHist[sd] = new TH1D(Form("f3pCorrelatorVs%sHist",psdFlag2[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
1272   f3pCorrelatorVsEtaSumDiffHist[sd]->SetStats(kFALSE);
1273   //f3pCorrelatorVsPtSumDiffHist[sd]->SetLabelSize(0.05);
1274   //f3pCorrelatorVsPtSumDiffHist[sd]->SetMarkerStyle(25);
1275   if(fHarmonic == 1)
1276   {
1277    f3pCorrelatorVsPtSumDiffHist[sd]->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} %s",psdTitleFlag[sd].Data()));
1278    f3pCorrelatorVsEtaSumDiffHist[sd]->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} %s",psdTitleFlag2[sd].Data()));
1279   } else
1280     {
1281      f3pCorrelatorVsPtSumDiffHist[sd]->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[sd].Data()));
1282      f3pCorrelatorVsEtaSumDiffHist[sd]->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[sd].Data()));
1283     }
1284   f3pCorrelatorVsPtSumDiffHist[sd]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
1285   fResultsList->Add(f3pCorrelatorVsPtSumDiffHist[sd]);
1286   f3pCorrelatorVsEtaSumDiffHist[sd]->GetXaxis()->SetTitle(psdTitleFlag2[sd].Data());
1287   fResultsList->Add(f3pCorrelatorVsEtaSumDiffHist[sd]);
1288  }
1289
1290  //TString psdFlag[2] = {"PtSum","PtDiff"};
1291  //TString psdTitleFlag[2] = {"(p_{T,1}+ p_{T,2})/2","#left|p_{T,1}- p_{T,2}#right|"};
1292  //TString psdFlag2[2] = {"EtaSum","EtaDiff"};
1293  //TString psdTitleFlag2[2] = {"(#eta_{1}+ #eta_{2})/2","#left|#eta_{1}- #eta_{2}#right|"};
1294  TString nonIsotropicTerm[10] = {"#LT#LTcos(#psi_{POI_1})#GT#GT","#LT#LTsin(#psi_{POI_1})#GT#GT",
1295                                  "#LT#LTcos(#psi_{POI_2})#GT#GT","#LT#LTsin(#psi_{POI_2})#GT#GT",
1296                                  "#LT#LTcos(#psi_{POI_1}-2#phi_{RP})#GT#GT","#LT#LTsin(#psi_{POI_1}-2#phi_{RP})#GT#GT",
1297                                  "#LT#LTcos(#psi_{POI_2}-2#phi_{RP})#GT#GT","#LT#LTsin(#psi_{POI_2}-2#phi_{RP})#GT#GT",
1298                                  "#LT#LTcos(#psi_{POI_1}+#psi_{POI_2})#GT#GT","#LT#LTsin(#psi_{POI_1}+#psi_{POI_2})#GT#GT"};
1299  for(Int_t sd=0;sd<2;sd++)
1300  {
1301   for(Int_t t=0;t<10;t++)
1302   {
1303    // Pt:
1304    fNonIsotropicTermsVsPtSumDiffPro[sd][t] = new TProfile(Form("fNonIsotropicTermsVs%sPro %s",psdFlag[sd].Data(),nonIsotropicTerm[t].Data()),"",fnBinsPt,0.,fPtMax);
1305    fNonIsotropicTermsVsPtSumDiffPro[sd][t]->SetTitle(Form("%s vs %s",nonIsotropicTerm[t].Data(),psdTitleFlag[sd].Data()));
1306    fNonIsotropicTermsVsPtSumDiffPro[sd][t]->SetStats(kFALSE);
1307    fNonIsotropicTermsVsPtSumDiffPro[sd][t]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
1308    fNonIsotropicTermsList->Add(fNonIsotropicTermsVsPtSumDiffPro[sd][t]);
1309    // Eta:
1310    fNonIsotropicTermsVsEtaSumDiffPro[sd][t] = new TProfile(Form("fNonIsotropicTermsVs%sPro %s",psdFlag2[sd].Data(),nonIsotropicTerm[t].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
1311    fNonIsotropicTermsVsEtaSumDiffPro[sd][t]->SetTitle(Form("%s vs %s",nonIsotropicTerm[t].Data(),psdTitleFlag2[sd].Data()));
1312    fNonIsotropicTermsVsEtaSumDiffPro[sd][t]->SetStats(kFALSE);
1313    fNonIsotropicTermsVsEtaSumDiffPro[sd][t]->GetXaxis()->SetTitle(psdTitleFlag2[sd].Data());
1314    fNonIsotropicTermsList->Add(fNonIsotropicTermsVsEtaSumDiffPro[sd][t]);
1315   } // end of for(Int_t t=0;t<10;t++)
1316  } // end of for(Int_t sd=0;sd<2;sd++)
1317
1318  //2p correlator vs |Pt1-Pt2|
1319  f2pCorrelatorCosPsiDiffPtDiff = new TProfile("f2pCorrelatorCosPsiDiffPtDiff",";|p_{T,1}-p_{T,2}| (GeV/c);#LT #LT cos(#psi_{1} - #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1320  f2pCorrelatorCosPsiDiffPtDiff->SetStats(kFALSE);
1321  f2pCorrelatorCosPsiSumPtDiff = new TProfile("f2pCorrelatorCosPsiSumPtDiff",";|p_{T,1}-p_{T,2}| (GeV/c);#LT #LT cos(#psi_{1} + #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1322  f2pCorrelatorCosPsiSumPtDiff->SetStats(kFALSE);
1323  f2pCorrelatorSinPsiDiffPtDiff = new TProfile("f2pCorrelatorSinPsiDiffPtDiff",";|p_{T,1}-p_{T,2}| (GeV/c);#LT #LT sin(#psi_{1} - #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1324  f2pCorrelatorSinPsiDiffPtDiff->SetStats(kFALSE);
1325  f2pCorrelatorSinPsiSumPtDiff = new TProfile("f2pCorrelatorSinPsiSumPtDiff",";|p_{T,1}-p_{T,2}| (GeV/c);#LT #LT sin(#psi_{1} + #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1326  f2pCorrelatorSinPsiSumPtDiff->SetStats(kFALSE);
1327  if(fHarmonic == 1) {
1328    f2pCorrelatorCosPsiDiffPtDiff->SetTitle(Form("#LT#LTcos(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1329    f2pCorrelatorCosPsiSumPtDiff->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1330    f2pCorrelatorSinPsiDiffPtDiff->SetTitle(Form("#LT#LTsin(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1331    f2pCorrelatorSinPsiSumPtDiff->SetTitle(Form("#LT#LTsin(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1332  }
1333  else {
1334    f2pCorrelatorCosPsiDiffPtDiff->SetTitle(Form("#LT#LTcos[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[1].Data()));
1335    f2pCorrelatorCosPsiSumPtDiff->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[1].Data()));
1336    f2pCorrelatorSinPsiDiffPtDiff->SetTitle(Form("#LT#LTsin[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[1].Data()));
1337    f2pCorrelatorSinPsiSumPtDiff->SetTitle(Form("#LT#LTsin[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[1].Data()));
1338  }
1339  fProfileList->Add(f2pCorrelatorCosPsiDiffPtDiff);
1340  fProfileList->Add(f2pCorrelatorCosPsiSumPtDiff);
1341  fProfileList->Add(f2pCorrelatorSinPsiDiffPtDiff);
1342  fProfileList->Add(f2pCorrelatorSinPsiSumPtDiff);
1343
1344  //2p correlator vs (Pt1+Pt2)/2
1345  f2pCorrelatorCosPsiDiffPtSum = new TProfile("f2pCorrelatorCosPsiDiffPtSum",";(p_{T,1}+p_{T,2})/2 (GeV/c);#LT #LT cos(#psi_{1} - #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1346  f2pCorrelatorCosPsiDiffPtSum->SetStats(kFALSE);
1347  f2pCorrelatorCosPsiSumPtSum = new TProfile("f2pCorrelatorCosPsiSumPtSum",";(p_{T,1}+p_{T,2})/2 (GeV/c);#LT #LT cos(#psi_{1} + #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1348  f2pCorrelatorCosPsiSumPtSum->SetStats(kFALSE);
1349  f2pCorrelatorSinPsiDiffPtSum = new TProfile("f2pCorrelatorSinPsiDiffPtSum",";(p_{T,1}+p_{T,2})/2 (GeV/c);#LT #LT sin(#psi_{1} - #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1350  f2pCorrelatorSinPsiDiffPtSum->SetStats(kFALSE);
1351  f2pCorrelatorSinPsiSumPtSum = new TProfile("f2pCorrelatorSinPsiSumPtSum",";(p_{T,1}+p_{T,2})/2 (GeV/c);#LT #LT sin(#psi_{1} + #psi_{2}) #GT #GT",fnBinsPt,0.,fPtMax);
1352  f2pCorrelatorSinPsiSumPtSum->SetStats(kFALSE);
1353  if(fHarmonic == 1) {
1354    f2pCorrelatorCosPsiDiffPtSum->SetTitle(Form("#LT#LTcos(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[0].Data()));
1355    f2pCorrelatorCosPsiSumPtSum->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[0].Data()));
1356    f2pCorrelatorSinPsiDiffPtSum->SetTitle(Form("#LT#LTsin(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[0].Data()));
1357    f2pCorrelatorSinPsiSumPtSum->SetTitle(Form("#LT#LTsin(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[0].Data()));
1358  }
1359  else {
1360    f2pCorrelatorCosPsiDiffPtSum->SetTitle(Form("#LT#LTcos[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[0].Data()));
1361    f2pCorrelatorCosPsiSumPtSum->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[0].Data()));
1362    f2pCorrelatorSinPsiDiffPtSum->SetTitle(Form("#LT#LTsin[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[0].Data()));
1363    f2pCorrelatorSinPsiSumPtSum->SetTitle(Form("#LT#LTsin[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[0].Data()));
1364  }
1365  fProfileList->Add(f2pCorrelatorCosPsiDiffPtSum);
1366  fProfileList->Add(f2pCorrelatorCosPsiSumPtSum);
1367  fProfileList->Add(f2pCorrelatorSinPsiDiffPtSum);
1368  fProfileList->Add(f2pCorrelatorSinPsiSumPtSum);
1369
1370  //2p correlator vs |eta1-eta2|
1371  f2pCorrelatorCosPsiDiffEtaDiff = new TProfile("f2pCorrelatorCosPsiDiffEtaDiff",";|#eta_{1}-#eta_{2}|;#LT #LT cos(#psi_{1} - #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1372  f2pCorrelatorCosPsiDiffEtaDiff->SetStats(kFALSE);
1373  f2pCorrelatorCosPsiSumEtaDiff = new TProfile("f2pCorrelatorCosPsiSumEtaDiff",";|#eta_{1}-#eta_{2}|;#LT #LT cos(#psi_{1} + #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1374  f2pCorrelatorCosPsiSumEtaDiff->SetStats(kFALSE);
1375  f2pCorrelatorSinPsiDiffEtaDiff = new TProfile("f2pCorrelatorSinPsiDiffEtaDiff",";|#eta_{1}-#eta_{2}|;#LT #LT sin(#psi_{1} - #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1376  f2pCorrelatorSinPsiDiffEtaDiff->SetStats(kFALSE);
1377  f2pCorrelatorSinPsiSumEtaDiff = new TProfile("f2pCorrelatorSinPsiSumEtaDiff",";|#eta_{1}-#eta_{2}|;#LT #LT sin(#psi_{1} + #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1378  f2pCorrelatorSinPsiSumEtaDiff->SetStats(kFALSE);
1379  if(fHarmonic == 1) {
1380    f2pCorrelatorCosPsiDiffEtaDiff->SetTitle(Form("#LT#LTcos(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[1].Data()));
1381    f2pCorrelatorCosPsiSumEtaDiff->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[1].Data()));
1382    f2pCorrelatorSinPsiDiffEtaDiff->SetTitle(Form("#LT#LTsin(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[1].Data()));
1383    f2pCorrelatorSinPsiSumEtaDiff->SetTitle(Form("#LT#LTsin(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[1].Data()));
1384  }
1385  else {
1386    f2pCorrelatorCosPsiDiffEtaDiff->SetTitle(Form("#LT#LTcos[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1387    f2pCorrelatorCosPsiSumEtaDiff->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1388    f2pCorrelatorSinPsiDiffEtaDiff->SetTitle(Form("#LT#LTsin[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1389    f2pCorrelatorSinPsiSumEtaDiff->SetTitle(Form("#LT#LTsin[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1390  }
1391  fProfileList->Add(f2pCorrelatorCosPsiDiffEtaDiff);
1392  fProfileList->Add(f2pCorrelatorCosPsiSumEtaDiff);
1393  fProfileList->Add(f2pCorrelatorSinPsiDiffEtaDiff);
1394  fProfileList->Add(f2pCorrelatorSinPsiSumEtaDiff);
1395
1396  //2p correlator vs (eta1+eta2)/2
1397  f2pCorrelatorCosPsiDiffEtaSum = new TProfile("f2pCorrelatorCosPsiDiffEtaSum",";(#eta_{1}+#eta_{,2})/2;#LT #LT cos(#psi_{1} - #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1398  f2pCorrelatorCosPsiDiffEtaSum->SetStats(kFALSE);
1399  f2pCorrelatorCosPsiSumEtaSum = new TProfile("f2pCorrelatorCosPsiSumEtaSum",";(#eta_{1}+#eta_{,2})/2;#LT #LT cos(#psi_{1} + #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1400  f2pCorrelatorCosPsiSumEtaSum->SetStats(kFALSE);
1401  f2pCorrelatorSinPsiDiffEtaSum = new TProfile("f2pCorrelatorSinPsiDiffEtaSum",";(#eta_{1}+#eta_{,2})/2;#LT #LT sin(#psi_{1} - #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1402  f2pCorrelatorSinPsiDiffEtaSum->SetStats(kFALSE);
1403  f2pCorrelatorSinPsiSumEtaSum = new TProfile("f2pCorrelatorSinPsiSumEtaSum",";(#eta_{1}+#eta_{,2})/2;#LT #LT sin(#psi_{1} + #psi_{2}) #GT #GT",fnBinsEta,fEtaMin,fEtaMax);
1404  f2pCorrelatorSinPsiSumEtaSum->SetStats(kFALSE);
1405  if(fHarmonic == 1) {
1406    f2pCorrelatorCosPsiDiffEtaSum->SetTitle(Form("#LT#LTcos(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[0].Data()));
1407    f2pCorrelatorCosPsiSumEtaSum->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[0].Data()));
1408    f2pCorrelatorSinPsiDiffEtaSum->SetTitle(Form("#LT#LTsin(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[0].Data()));
1409    f2pCorrelatorSinPsiSumEtaSum->SetTitle(Form("#LT#LTsin(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag2[0].Data()));
1410  }
1411  else {
1412    f2pCorrelatorCosPsiDiffEtaSum->SetTitle(Form("#LT#LTcos[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[0].Data()));
1413    f2pCorrelatorCosPsiSumEtaSum->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[0].Data()));
1414    f2pCorrelatorSinPsiDiffEtaSum->SetTitle(Form("#LT#LTsin[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[0].Data()));
1415    f2pCorrelatorSinPsiSumEtaSum->SetTitle(Form("#LT#LTsin[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[0].Data()));
1416  }
1417  fProfileList->Add(f2pCorrelatorCosPsiDiffEtaSum);
1418  fProfileList->Add(f2pCorrelatorCosPsiSumEtaSum);
1419  fProfileList->Add(f2pCorrelatorSinPsiDiffEtaSum);
1420  fProfileList->Add(f2pCorrelatorSinPsiSumEtaSum);
1421
1422 } // end of void AliFlowAnalysisWithMixedHarmonics::BookDifferential()
1423
1424 //================================================================================================================
1425
1426 void AliFlowAnalysisWithMixedHarmonics::AccessConstants(TString method)
1427 {
1428  // Access and store common constants.
1429
1430  // a) If this method was called in Init() access common constants from AliFlowCommonConstants;
1431  // b) If this method was called in Init() book and fill TProfile to hold constants accessed in a);
1432  // c) If this method was called in Finish() access common constants from TProfile booked and filled in b).
1433
1434  if(method == "Init")
1435  {
1436   // a) If this method was called in Init() access common constants from AliFlowCommonConstants:
1437   fnBinsPhi = AliFlowCommonConstants::GetMaster()->GetNbinsPhi();
1438   fPhiMin = AliFlowCommonConstants::GetMaster()->GetPhiMin();
1439   fPhiMax = AliFlowCommonConstants::GetMaster()->GetPhiMax();
1440   if(fnBinsPhi){fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;}
1441   fnBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
1442   fPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
1443   fPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
1444   if(fnBinsPt){fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;}
1445   fnBinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
1446   fEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
1447   fEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
1448   if(fnBinsEta){fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;}
1449
1450   // b) If this method was called in Init() book and fill TProfile to hold constants accessed in a):
1451   TString fCommonConstantsName = "fCommonConstants";
1452   fCommonConstants = new TProfile(fCommonConstantsName.Data(),"Common constants",9,0.,9.);
1453   fCommonConstants->SetLabelSize(0.05);
1454   fCommonConstants->GetXaxis()->SetBinLabel(1,"nBins (#phi)");
1455   fCommonConstants->Fill(0.5,fnBinsPhi);
1456   fCommonConstants->GetXaxis()->SetBinLabel(2,"#phi_{min}");
1457   fCommonConstants->Fill(1.5,fPhiMin);
1458   fCommonConstants->GetXaxis()->SetBinLabel(3,"#phi_{max}");
1459   fCommonConstants->Fill(2.5,fPhiMax);
1460   fCommonConstants->GetXaxis()->SetBinLabel(4,"nBins (p_{t})");
1461   fCommonConstants->Fill(3.5,fnBinsPt);
1462   fCommonConstants->GetXaxis()->SetBinLabel(5,"(p_{t})_{min}");
1463   fCommonConstants->Fill(4.5,fPtMin);
1464   fCommonConstants->GetXaxis()->SetBinLabel(6,"(p_{t})_{max}");
1465   fCommonConstants->Fill(5.5,fPtMax);
1466   fCommonConstants->GetXaxis()->SetBinLabel(7,"nBins (#eta)");
1467   fCommonConstants->Fill(6.5,fnBinsEta);
1468   fCommonConstants->GetXaxis()->SetBinLabel(8,"#eta_{min}");
1469   fCommonConstants->Fill(7.5,fEtaMin);
1470   fCommonConstants->GetXaxis()->SetBinLabel(9,"#eta_{max}");
1471   fCommonConstants->Fill(8.5,fEtaMax);
1472   fHistList->Add(fCommonConstants);
1473  } // end of if(method == "Init")
1474  else if(method == "Finish")
1475  {
1476   // c) If this method was called in Finish() access common constants from TProfile booked and filled in b):
1477   if(!fCommonConstants)
1478   {
1479    printf("\n WARNING (MH): fCommonConstants is NULL in AFAWMH::AC(\"%s\") !!!!\n\n",method.Data());
1480    exit(0);
1481   }
1482   fnBinsPhi = (Int_t)fCommonConstants->GetBinContent(1);
1483   fPhiMin = fCommonConstants->GetBinContent(2);
1484   fPhiMax = fCommonConstants->GetBinContent(3);
1485   if(fnBinsPhi){fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;}
1486   fnBinsPt = (Int_t)fCommonConstants->GetBinContent(4);
1487   fPtMin = fCommonConstants->GetBinContent(5);
1488   fPtMax = fCommonConstants->GetBinContent(6);
1489   if(fnBinsPt){fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;}
1490   fnBinsEta = (Int_t)fCommonConstants->GetBinContent(7);
1491   fEtaMin = fCommonConstants->GetBinContent(8);
1492   fEtaMax = fCommonConstants->GetBinContent(9);
1493   if(fnBinsEta){fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;}
1494  } // end of else if(method == "Finish")
1495
1496 } // end of void AliFlowAnalysisWithMixedHarmonics::AccessConstants(TString method)
1497
1498 //================================================================================================================
1499
1500 void AliFlowAnalysisWithMixedHarmonics::CrossCheckSettings()
1501 {
1502  // Cross-check if the user settings make sense.
1503
1504  // ...
1505
1506 } // end of void AliFlowAnalysisWithMixedHarmonics::CrossCheckSettings()
1507
1508 //================================================================================================================
1509
1510 void AliFlowAnalysisWithMixedHarmonics::BookAndFillWeightsHistograms()
1511 {
1512  // Book and fill (by accessing file "weights.root") histograms which hold phi, pt and eta weights.
1513
1514  if(!fWeightsList)
1515  {
1516   cout<<endl;
1517   cout<<" WARNING (MH): fWeightsList is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1518   cout<<endl;
1519   exit(0);
1520  }
1521  // Profile to hold flags for weights:
1522  TString fUseParticleWeightsName = "fUseParticleWeightsMH";
1523  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
1524  fUseParticleWeights->SetStats(kFALSE);
1525  fUseParticleWeights->SetLabelSize(0.06);
1526  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
1527  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
1528  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
1529  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
1530  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
1531  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
1532  fWeightsList->Add(fUseParticleWeights);
1533  // Phi-weights:
1534  if(fUsePhiWeights)
1535  {
1536   if(fWeightsList->FindObject("phi_weights"))
1537   {
1538    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
1539    if (!fPhiWeights)
1540    {
1541      printf("WARNING: no phi weights. bye!\n");
1542      exit(0);
1543    }
1544    if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
1545    {
1546     cout<<endl;
1547     cout<<" WARNING (MH): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
1548     cout<<endl;
1549     exit(0);
1550    }
1551   } else
1552     {
1553      cout<<endl;
1554      cout<<" WARNING (MH): fWeightsList->FindObject(\"phi_weights\") is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1555      cout<<endl;
1556      exit(0);
1557     }
1558  } // end of if(fUsePhiWeights)
1559  // Pt-weights:
1560  if(fUsePtWeights)
1561  {
1562   if(fWeightsList->FindObject("pt_weights"))
1563   {
1564    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
1565    if (!fPtWeights)
1566    {
1567      printf("WARNING: no pt weights. bye!\n");
1568      exit(0);
1569    }
1570    if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
1571    {
1572     cout<<endl;
1573     cout<<" WARNING (MH): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
1574     cout<<endl;
1575     exit(0);
1576    }
1577   } else
1578     {
1579      cout<<endl;
1580      cout<<" WARNING (MH): fWeightsList->FindObject(\"pt_weights\") is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1581      cout<<endl;
1582      exit(0);
1583     }
1584  } // end of if(fUsePtWeights)
1585  // Eta-weights:
1586  if(fUseEtaWeights)
1587  {
1588   if(fWeightsList->FindObject("eta_weights"))
1589   {
1590    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
1591    if (!fEtaWeights)
1592    {
1593      printf("WARNING: no pt weights. bye!\n");
1594      exit(0);
1595    }
1596    if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
1597    {
1598     cout<<endl;
1599     cout<<" WARNING (MH): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
1600     cout<<endl;
1601     exit(0);
1602    }
1603   } else
1604     {
1605      cout<<endl;
1606      cout<<" WARNING (MH): fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1607      cout<<endl;
1608      exit(0);
1609     }
1610  } // end of if(fUseEtaWeights)
1611
1612 } // end of AliFlowAnalysisWithMixedHarmonics::BookAndFillWeightsHistograms()
1613
1614 //================================================================================================================
1615
1616 void AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInMake()
1617 {
1618  // Check pointers used in method Make().
1619
1620  if(!fReQnk || !fImQnk || !fSpk )
1621  {
1622   cout<<endl;
1623   cout<<" WARNING (MH): fReQnk || fImQnk || fSpk is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1624   cout<<endl;
1625   exit(0);
1626  }
1627  if(!f3pCorrelatorPro)
1628  {
1629   cout<<endl;
1630   cout<<" WARNING (MH): f3pCorrelatorPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1631   cout<<endl;
1632   exit(0);
1633  }
1634  if(!f5pCorrelatorPro)
1635  {
1636   cout<<endl;
1637   cout<<" WARNING (MH): f5pCorrelatorPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1638   cout<<endl;
1639   exit(0);
1640  }
1641  if(!fNonIsotropicTermsPro)
1642  {
1643   cout<<endl;
1644   cout<<" WARNING (MH): fNonIsotropicTermsPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1645   cout<<endl;
1646   exit(0);
1647  }
1648  if(!f3pCorrelatorVsMPro && fCalculateVsM)
1649  {
1650   cout<<endl;
1651   cout<<" WARNING (MH): f3pCorrelatorVsMPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1652   cout<<endl;
1653   exit(0);
1654  }
1655  if(!f3pPOICorrelatorVsM && fCalculateVsM)
1656  {
1657   cout<<endl;
1658   cout<<" WARNING (MH): f3pPOICorrelatorVsM is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1659   cout<<endl;
1660   exit(0);
1661  }
1662  if(!fNonIsotropicTermsVsMPro && fCalculateVsM)
1663  {
1664   cout<<endl;
1665   cout<<" WARNING (MH): fNonIsotropicTermsVsMPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1666   cout<<endl;
1667   exit(0);
1668  }
1669
1670  // Differential correlators:
1671  if(!fEvaluateDifferential3pCorrelator){return;}
1672  for(Int_t sd=0;sd<2;sd++)
1673  {
1674   if(!(f3pCorrelatorVsPtSumDiffPro[sd]))
1675   {
1676    cout<<endl;
1677    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsPtSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1678    cout<<endl;
1679    exit(0);
1680   }
1681   if(!(f3pCorrelatorVsEtaSumDiffPro[sd]))
1682   {
1683    cout<<endl;
1684    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsEtaSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1685    cout<<endl;
1686    exit(0);
1687   }
1688   for(Int_t t=0;t<10;t++)
1689   {
1690    if(!(fNonIsotropicTermsVsPtSumDiffPro[sd][t]))
1691    {
1692     cout<<endl;
1693     cout<<" WARNING (MH): "<<Form("fNonIsotropicTermsVsPtSumDiffPro[%d][%d]",sd,t)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1694     cout<<endl;
1695     exit(0);
1696    }
1697    if(!(fNonIsotropicTermsVsEtaSumDiffPro[sd][t]))
1698    {
1699     cout<<endl;
1700     cout<<" WARNING (MH): "<<Form("fNonIsotropicTermsVsEtaSumDiffPro[%d][%d]",sd,t)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1701     cout<<endl;
1702     exit(0);
1703    }
1704   } // end of for(Int_t t=0;t<10;t++)
1705  } // end of for(Int_t sd=0;sd<2;sd++)
1706  for(Int_t sd=0;sd<2;sd++)
1707  {
1708   if(!fRePEBE[sd]||!fImPEBE[sd])
1709   {
1710    cout<<endl;
1711    cout<<" WARNING (MH): "<<Form("!fRePEBE[%d]||!fImPEBE[%d]",sd,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1712    cout<<endl;
1713    exit(0);
1714   }
1715   if(!fReEtaEBE[sd]||!fImEtaEBE[sd])
1716   {
1717    cout<<endl;
1718    cout<<" WARNING (MH): "<<Form("!fReEtaEBE[%d]||!fImEtaEBE[%d]",sd,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1719    cout<<endl;
1720    exit(0);
1721   }
1722   for(Int_t fs=0;fs<2;fs++)
1723   {
1724    if(!fOverlapEBE[fs][sd]||!fOverlapEBE[fs][sd])
1725    {
1726     cout<<endl;
1727     cout<<" WARNING (MH): "<<Form("!fOverlapEBE[%d][%d]||!fOverlapEBE[%d][%d]",fs,sd,fs,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1728     cout<<endl;
1729     exit(0);
1730    }
1731    if(!fOverlapEBE2[fs][sd]||!fOverlapEBE2[fs][sd])
1732    {
1733     cout<<endl;
1734     cout<<" WARNING (MH): "<<Form("!fOverlapEBE2[%d][%d]||!fOverlapEBE2[%d][%d]",fs,sd,fs,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1735     cout<<endl;
1736     exit(0);
1737    }
1738   } // end of for(Int_t fs=0;fs<2;fs++)
1739  } // end of for(Int_t sd=0;sd<2;sd++)
1740  for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
1741  {
1742   for(Int_t ao=0;ao<2;ao++) // all/overlap
1743   {
1744    for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
1745    {
1746     if(!fReNITEBE[p12][ao][pe]||!fImNITEBE[p12][ao][pe])
1747     {
1748      cout<<endl;
1749      cout<<" WARNING (MH): "<<Form("!fReNITEBE[%d][%d][%d]||!fImNITEBE[%d][%d][%d]",p12,ao,pe,p12,ao,pe)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1750      cout<<endl;
1751      exit(0);
1752     }
1753    } // end of for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
1754   } // end of for(Int_t ao=0;ao<2;ao++) // all/overlap
1755  } // end of for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
1756
1757  //2p correlator vs |Pt1-Pt2|
1758  if(!f2pCorrelatorCosPsiDiffPtDiff) {
1759   cout<<endl;
1760   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffPtDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1761   cout<<endl;
1762   exit(0);
1763  }
1764  if(!f2pCorrelatorCosPsiSumPtDiff) {
1765   cout<<endl;
1766   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumPtDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1767   cout<<endl;
1768   exit(0);
1769  }
1770  if(!f2pCorrelatorSinPsiDiffPtDiff) {
1771   cout<<endl;
1772   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffPtDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1773   cout<<endl;
1774   exit(0);
1775  }
1776  if(!f2pCorrelatorSinPsiSumPtDiff) {
1777   cout<<endl;
1778   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumPtDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1779   cout<<endl;
1780   exit(0);
1781  }
1782
1783  //2p correlator vs (Pt1+Pt2)/2
1784  if(!f2pCorrelatorCosPsiDiffPtSum) {
1785   cout<<endl;
1786   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffPtSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1787   cout<<endl;
1788   exit(0);
1789  }
1790  if(!f2pCorrelatorCosPsiSumPtSum) {
1791   cout<<endl;
1792   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumPtSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1793   cout<<endl;
1794   exit(0);
1795  }
1796  if(!f2pCorrelatorSinPsiDiffPtSum) {
1797   cout<<endl;
1798   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffPtSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1799   cout<<endl;
1800   exit(0);
1801  }
1802  if(!f2pCorrelatorSinPsiSumPtSum) {
1803   cout<<endl;
1804   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumPtSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1805   cout<<endl;
1806   exit(0);
1807  }
1808
1809  //2p correlator vs |eta1-eta2|
1810  if(!f2pCorrelatorCosPsiDiffEtaDiff) {
1811   cout<<endl;
1812   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffEtaDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1813   cout<<endl;
1814   exit(0);
1815  }
1816  if(!f2pCorrelatorCosPsiSumEtaDiff) {
1817   cout<<endl;
1818   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumEtaDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1819   cout<<endl;
1820   exit(0);
1821  }
1822  if(!f2pCorrelatorSinPsiDiffEtaDiff) {
1823   cout<<endl;
1824   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffEtaDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1825   cout<<endl;
1826   exit(0);
1827  }
1828  if(!f2pCorrelatorSinPsiSumEtaDiff) {
1829   cout<<endl;
1830   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumEtaDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1831   cout<<endl;
1832   exit(0);
1833  }
1834
1835  //2p correlator vs (eta1+eta2)/2
1836  if(!f2pCorrelatorCosPsiDiffEtaSum) {
1837   cout<<endl;
1838   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffEtaSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1839   cout<<endl;
1840   exit(0);
1841  }
1842  if(!f2pCorrelatorCosPsiSumEtaSum) {
1843   cout<<endl;
1844   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumEtaSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1845   cout<<endl;
1846   exit(0);
1847  }
1848  if(!f2pCorrelatorSinPsiDiffEtaSum) {
1849   cout<<endl;
1850   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffEtaSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1851   cout<<endl;
1852   exit(0);
1853  }
1854  if(!f2pCorrelatorSinPsiSumEtaSum) {
1855   cout<<endl;
1856   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumEtaSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1857   cout<<endl;
1858   exit(0);
1859  }
1860
1861 } // end of AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInMake()
1862
1863 //================================================================================================================
1864
1865 void AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInFinish()
1866 {
1867  // Check pointers used in method Finish().
1868
1869  if(!fAnalysisSettings)
1870  {
1871   cout<<endl;
1872   cout<<" WARNING (MH): fAnalysisSettings is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1873   cout<<endl;
1874   exit(0);
1875  }
1876  if(!f3pCorrelatorPro)
1877  {
1878   cout<<endl;
1879   cout<<" WARNING (MH): f3pCorrelatorPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1880   cout<<endl;
1881   exit(0);
1882  }
1883  if(!fNonIsotropicTermsPro)
1884  {
1885   cout<<endl;
1886   cout<<" WARNING (MH): fNonIsotropicTermsPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1887   cout<<endl;
1888   exit(0);
1889  }
1890  if(!f3pPOICorrelatorVsM && fCalculateVsM)
1891  {
1892   cout<<endl;
1893   cout<<" WARNING (MH): f3pPOICorrelatorVsM is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1894   cout<<endl;
1895   exit(0);
1896  }
1897  if(!f3pCorrelatorVsMPro && fCalculateVsM)
1898  {
1899   cout<<endl;
1900   cout<<" WARNING (MH): f3pCorrelatorVsMPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1901   cout<<endl;
1902   exit(0);
1903  }
1904  if(!f3pCorrelatorVsMHist && fCalculateVsM)
1905  {
1906   cout<<endl;
1907   cout<<" WARNING (MH): f3pCorrelatorVsMHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1908   cout<<endl;
1909   exit(0);
1910  }
1911  if(!fNonIsotropicTermsVsMPro && fCalculateVsM)
1912  {
1913   cout<<endl;
1914   cout<<" WARNING (MH): fNonIsotropicTermsVsMPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1915   cout<<endl;
1916   exit(0);
1917  }
1918
1919  if(!f3pCorrelatorHist)
1920  {
1921   cout<<endl;
1922   cout<<" WARNING (MH): f3pCorrelatorHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1923   cout<<endl;
1924   exit(0);
1925  }
1926  if(!fDetectorBiasHist)
1927  {
1928   cout<<endl;
1929   cout<<" WARNING (MH): fDetectorBiasHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1930   cout<<endl;
1931   exit(0);
1932  }
1933  /* to be improved - enabled eventually
1934  if(!fDetectorBiasVsMHist && fCalculateVsM)
1935  {
1936   cout<<endl;
1937   cout<<" WARNING (MH): !fDetectorBiasVsMHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1938   cout<<endl;
1939   exit(0);
1940  }
1941  */
1942
1943  // Differential correlators:
1944  if(!fEvaluateDifferential3pCorrelator){return;}
1945  for(Int_t sd=0;sd<2;sd++)
1946  {
1947   if(!(f3pCorrelatorVsPtSumDiffPro[sd]))
1948   {
1949    cout<<endl;
1950    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsPtSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1951    cout<<endl;
1952    exit(0);
1953   }
1954   if(!(f3pCorrelatorVsEtaSumDiffPro[sd]))
1955   {
1956    cout<<endl;
1957    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsEtaSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1958    cout<<endl;
1959    exit(0);
1960   }
1961   if(!(f3pCorrelatorVsPtSumDiffHist[sd]))
1962   {
1963    cout<<endl;
1964    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsPtSumDiffHist[%d]",sd)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1965    cout<<endl;
1966    exit(0);
1967   }
1968   if(!(f3pCorrelatorVsEtaSumDiffHist[sd]))
1969   {
1970    cout<<endl;
1971    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsEtaSumDiffHist[%d]",sd)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1972    cout<<endl;
1973    exit(0);
1974   }
1975   for(Int_t t=0;t<10;t++)
1976   {
1977    if(!(fNonIsotropicTermsVsPtSumDiffPro[sd][t]))
1978    {
1979     cout<<endl;
1980     cout<<" WARNING (MH): "<<Form("fNonIsotropicTermsVsPtSumDiffPro[%d][%d]",sd,t)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1981     cout<<endl;
1982     exit(0);
1983    }
1984    if(!(fNonIsotropicTermsVsEtaSumDiffPro[sd][t]))
1985    {
1986     cout<<endl;
1987     cout<<" WARNING (MH): "<<Form("fNonIsotropicTermsVsEtaSumDiffPro[%d][%d]",sd,t)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1988     cout<<endl;
1989     exit(0);
1990    }
1991   } // end of for(Int_t t=0;t<10;t++)
1992  } // end of for(Int_t sd=0;sd<2;sd++)
1993  //2p correlator vs |Pt1-Pt2|
1994  if(!f2pCorrelatorCosPsiDiffPtDiff) {
1995   cout<<endl;
1996   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffPtDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1997   cout<<endl;
1998   exit(0);
1999  }
2000  if(!f2pCorrelatorCosPsiSumPtDiff) {
2001   cout<<endl;
2002   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumPtDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2003   cout<<endl;
2004   exit(0);
2005  }
2006  if(!f2pCorrelatorSinPsiDiffPtDiff) {
2007   cout<<endl;
2008   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffPtDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2009   cout<<endl;
2010   exit(0);
2011  }
2012  if(!f2pCorrelatorSinPsiSumPtDiff) {
2013   cout<<endl;
2014   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumPtDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2015   cout<<endl;
2016   exit(0);
2017  }
2018
2019  //2p correlator vs (Pt1+Pt2)/2
2020  if(!f2pCorrelatorCosPsiDiffPtSum) {
2021   cout<<endl;
2022   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffPtSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2023   cout<<endl;
2024   exit(0);
2025  }
2026  if(!f2pCorrelatorCosPsiSumPtSum) {
2027   cout<<endl;
2028   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumPtSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2029   cout<<endl;
2030   exit(0);
2031  }
2032  if(!f2pCorrelatorSinPsiDiffPtSum) {
2033   cout<<endl;
2034   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffPtSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2035   cout<<endl;
2036   exit(0);
2037  }
2038  if(!f2pCorrelatorSinPsiSumPtSum) {
2039   cout<<endl;
2040   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumPtSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2041   cout<<endl;
2042   exit(0);
2043  }
2044
2045  //2p correlator vs |eta1-eta2|
2046  if(!f2pCorrelatorCosPsiDiffEtaDiff) {
2047   cout<<endl;
2048   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffEtaDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2049   cout<<endl;
2050   exit(0);
2051  }
2052  if(!f2pCorrelatorCosPsiSumEtaDiff) {
2053   cout<<endl;
2054   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumEtaDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2055   cout<<endl;
2056   exit(0);
2057  }
2058  if(!f2pCorrelatorSinPsiDiffEtaDiff) {
2059   cout<<endl;
2060   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffEtaDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2061   cout<<endl;
2062   exit(0);
2063  }
2064  if(!f2pCorrelatorSinPsiSumEtaDiff) {
2065   cout<<endl;
2066   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumEtaDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2067   cout<<endl;
2068   exit(0);
2069  }
2070
2071  //2p correlator vs (eta1+eta2)/2
2072  if(!f2pCorrelatorCosPsiDiffEtaSum) {
2073   cout<<endl;
2074   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiffEtaSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2075   cout<<endl;
2076   exit(0);
2077  }
2078  if(!f2pCorrelatorCosPsiSumEtaSum) {
2079   cout<<endl;
2080   cout<<" WARNING (MH): f2pCorrelatorCosPsiSumEtaSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2081   cout<<endl;
2082   exit(0);
2083  }
2084  if(!f2pCorrelatorSinPsiDiffEtaSum) {
2085   cout<<endl;
2086   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiffEtaSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2087   cout<<endl;
2088   exit(0);
2089  }
2090  if(!f2pCorrelatorSinPsiSumEtaSum) {
2091   cout<<endl;
2092   cout<<" WARNING (MH): f2pCorrelatorSinPsiSumEtaSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
2093   cout<<endl;
2094   exit(0);
2095  }
2096
2097 } // end of AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInFinish()
2098
2099 //================================================================================================================
2100
2101 void AliFlowAnalysisWithMixedHarmonics::PrintOnTheScreen()
2102 {
2103  // Print the final results on the screen.
2104
2105  // <<cos[n*(phi1+phi2-2phi3)]>>:
2106  Double_t d3pCorrelator = 0.;
2107  Double_t d3pCorrelatorError = 0.;
2108  if(!fCorrectForDetectorEffects)
2109  {
2110   d3pCorrelator = f3pCorrelatorPro->GetBinContent(1);
2111   d3pCorrelatorError = f3pCorrelatorPro->GetBinError(1);
2112  } else
2113    {
2114     d3pCorrelator = f3pCorrelatorHist->GetBinContent(1);
2115     d3pCorrelatorError = f3pCorrelatorHist->GetBinError(1);
2116    }
2117
2118  // <<cos[n*(psi1+psi2-2phi3)]>>:
2119  Double_t d3pCorrelatorPoi = 0.;
2120  Double_t d3pCorrelatorPoiError = 0.;
2121
2122  if(fEvaluateDifferential3pCorrelator)
2123  {
2124   GetCorrelatorAndError(f3pCorrelatorVsPtSumDiffPro[0],
2125                        d3pCorrelatorPoi,
2126                        d3pCorrelatorPoiError);
2127  }
2128  cout<<endl;
2129  cout<<"*******************************************************"<<endl;
2130  cout<<"*******************************************************"<<endl;
2131  cout<<"                    Mixed Harmonics                      "<<endl;
2132  cout<<endl;
2133  if(fHarmonic!=1)
2134  {
2135   cout<<"  cos["<<fHarmonic<<"(phi1+phi2-2phi3)] = "<<d3pCorrelator<<" +/- "<<d3pCorrelatorError<<endl;
2136   cout<<"  cos["<<fHarmonic<<"(psi1+psi2-2phi3)] = "<<d3pCorrelatorPoi<<" +/- "<<d3pCorrelatorPoiError<<endl;
2137  } else
2138    {
2139     cout<<"  cos(phi1+phi2-2phi3) = "<<d3pCorrelator<<" +/- "<<d3pCorrelatorError<<endl;
2140     cout<<"  cos(psi1+psi2-2phi3) = "<<d3pCorrelatorPoi<<" +/- "<<d3pCorrelatorPoiError<<endl;
2141    }
2142  if(!fCorrectForDetectorEffects)
2143  {
2144   cout<<"  Detector Bias = "<<fDetectorBiasHist->GetBinContent(1)<<" (not corrected for)"<<endl;
2145  } else
2146    {
2147     cout<<"  Detector Bias = "<<fDetectorBiasHist->GetBinContent(1)<<" (corrected for)"<<endl;
2148    }
2149  cout<<endl;
2150  cout<<"             nEvts = "<<(Int_t)fCommonHists->GetHistMultRP()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()<<endl;
2151  cout<<"*******************************************************"<<endl;
2152  cout<<"*******************************************************"<<endl;
2153
2154 } // end of void AliFlowAnalysisWithMixedHarmonics::PrintOnTheScreen()
2155
2156 //================================================================================================================
2157
2158 void AliFlowAnalysisWithMixedHarmonics::AccessSettings()
2159 {
2160  // Access the settings for analysis with mixed harmonics.
2161
2162  fCorrectForDetectorEffects = (Bool_t)fAnalysisSettings->GetBinContent(1);
2163  fNoOfMultipicityBins = (Int_t)fAnalysisSettings->GetBinContent(2);
2164  fMultipicityBinWidth = (Double_t)fAnalysisSettings->GetBinContent(3);
2165  fMinMultiplicity = (Double_t)fAnalysisSettings->GetBinContent(4);
2166  fPrintOnTheScreen = (Bool_t)fAnalysisSettings->GetBinContent(5);
2167  fHarmonic = (Int_t)fAnalysisSettings->GetBinContent(6);
2168  fOppositeChargesPOI = (Bool_t)fAnalysisSettings->GetBinContent(7);
2169  fEvaluateDifferential3pCorrelator = (Bool_t)fAnalysisSettings->GetBinContent(8);
2170  fCalculateVsM = (Bool_t)fAnalysisSettings->GetBinContent(9);
2171  fShowBinLabelsVsM = (Bool_t)fAnalysisSettings->GetBinContent(10);
2172
2173 } // end of AliFlowAnalysisWithMixedHarmonics::AccessSettings()
2174
2175 //================================================================================================================
2176
2177 void AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffects()
2178 {
2179  // a.) Correct integrated 3-p correlator cos[n(phi1+phi2-2phi3)] for detector effects;
2180  // b.) Correct differential 3-p correlator cos[n(psi1+psi2-2phi3)] for detector effects.
2181
2182  // a.) Correct integrated 3-p correlator cos[n(phi1+phi2-2phi3)] for detector effects:
2183  Double_t measured3pCorrelator = f3pCorrelatorPro->GetBinContent(1); // biased by detector effects
2184  Double_t corrected3pCorrelator = 0.; // corrected for detector effects
2185  Double_t nonIsotropicTerms[10] = {0.}; // there are 10 distinct non-isotropic terms
2186  for(Int_t nit=0;nit<10;nit++)
2187  {
2188   nonIsotropicTerms[nit] = fNonIsotropicTermsPro->GetBinContent(nit+1);
2189  }
2190  // Calculate corrected 3-p correlator:
2191  corrected3pCorrelator = measured3pCorrelator
2192                        - nonIsotropicTerms[2]*nonIsotropicTerms[4]
2193                        - nonIsotropicTerms[3]*nonIsotropicTerms[5]
2194                        - 2.*nonIsotropicTerms[0]*nonIsotropicTerms[6]
2195                        - 2.*nonIsotropicTerms[1]*nonIsotropicTerms[7]
2196                        + 2.*nonIsotropicTerms[2]*(pow(nonIsotropicTerms[0],2.)-pow(nonIsotropicTerms[1],2.))
2197                        + 4.*nonIsotropicTerms[3]*nonIsotropicTerms[0]*nonIsotropicTerms[1];
2198  // Store corrected correlator:
2199  if(fCorrectForDetectorEffects)
2200  {
2201   f3pCorrelatorHist->SetBinContent(1,corrected3pCorrelator);
2202   f3pCorrelatorHist->SetBinError(1,f3pCorrelatorPro->GetBinError(1)); // to be improved (propagate error for non-isotropic terms)
2203  }
2204  // Quantify bias from detector inefficiences to 3-p correlator. Remark: Bias is quantified as a
2205  // ratio between corrected and measured 3-p correlator:
2206  //              bias = corrected/measured
2207  // This bias is stored in histogram fDetectorBias.
2208  Double_t bias = 0.;
2209  if(TMath::Abs(measured3pCorrelator)>1.e-44)
2210  {
2211   bias = corrected3pCorrelator/measured3pCorrelator;
2212   fDetectorBiasHist->SetBinContent(1,bias);
2213  }
2214
2215  if(!fEvaluateDifferential3pCorrelator){return;}
2216  // b.) Correct differential 3-p correlator cos[n(psi1+psi2-2phi3)] for detector effects:
2217  for(Int_t sd=0;sd<2;sd++)
2218  {
2219   // [(p1+p2)/2,|p1-p2|]
2220   // looping over all bins and calculating reduced correlations:
2221   for(Int_t b=1;b<=fnBinsPt;b++)
2222   {
2223    Double_t measured = f3pCorrelatorVsPtSumDiffPro[sd]->GetBinContent(b);
2224    Double_t measuredErr = f3pCorrelatorVsPtSumDiffPro[sd]->GetBinError(b);
2225    Double_t corrected = 0.; // 3-p correlator corrected for detector effects
2226    Double_t correctedErr = measuredErr; // to be improved - propagate error also for non-isotropic terms
2227    // non-isotropic terms:
2228    Double_t cosPsiPOI1 = fNonIsotropicTermsVsPtSumDiffPro[sd][0]->GetBinContent(b); // <<cos(#psi_{POI_1})>>
2229    Double_t sinPsiPOI1 = fNonIsotropicTermsVsPtSumDiffPro[sd][1]->GetBinContent(b); // <<sin(#psi_{POI_1})>>
2230    Double_t cosPsiPOI2 = fNonIsotropicTermsVsPtSumDiffPro[sd][2]->GetBinContent(b); // <<cos(#psi_{POI_2})>>
2231    Double_t sinPsiPOI2 = fNonIsotropicTermsVsPtSumDiffPro[sd][3]->GetBinContent(b); // <<sin(#psi_{POI_2})>>
2232    Double_t cosPsiPOI1m2PhiRP = fNonIsotropicTermsVsPtSumDiffPro[sd][4]->GetBinContent(b); // <<cos(#psi_{POI_1}-2*phi_{RP})>>
2233    Double_t sinPsiPOI1m2PhiRP = fNonIsotropicTermsVsPtSumDiffPro[sd][5]->GetBinContent(b); // <<sin(#psi_{POI_1}-2*phi_{RP})>>
2234    Double_t cosPsiPOI2m2PhiRP = fNonIsotropicTermsVsPtSumDiffPro[sd][6]->GetBinContent(b); // <<cos(#psi_{POI_2}-2*phi_{RP})>>
2235    Double_t sinPsiPOI2m2PhiRP = fNonIsotropicTermsVsPtSumDiffPro[sd][7]->GetBinContent(b); // <<sin(#psi_{POI_2}-2*phi_{RP})>>
2236    Double_t cosPsiPOI1pPsiPOI2 = fNonIsotropicTermsVsPtSumDiffPro[sd][8]->GetBinContent(b); // <<cos(#psi_{POI_1}+#psi_{POI_2})>>
2237    Double_t sinPsiPOI1pPsiPOI2 = fNonIsotropicTermsVsPtSumDiffPro[sd][9]->GetBinContent(b); // <<sin(#psi_{POI_1}+#psi_{POI_2})>>
2238    Double_t cos2PhiRP = fNonIsotropicTermsPro->GetBinContent(3); // <<cos(2n*(phi_{RP}))>>
2239    Double_t sin2PhiRP = fNonIsotropicTermsPro->GetBinContent(4); // <<sin(2n*(phi_{RP}))>>
2240    // apply correction:
2241    corrected = measured
2242              - (cosPsiPOI1*cosPsiPOI2m2PhiRP-sinPsiPOI1*sinPsiPOI2m2PhiRP
2243              + cosPsiPOI2*cosPsiPOI1m2PhiRP-sinPsiPOI2*sinPsiPOI1m2PhiRP
2244              + cos2PhiRP*cosPsiPOI1pPsiPOI2+sin2PhiRP*sinPsiPOI1pPsiPOI2)
2245              + 2.*cos2PhiRP*(cosPsiPOI1*cosPsiPOI2-sinPsiPOI1*sinPsiPOI2)
2246              + 2.*sin2PhiRP*(cosPsiPOI1*sinPsiPOI2+sinPsiPOI1*cosPsiPOI2);
2247    if(fCorrectForDetectorEffects)
2248    {
2249     f3pCorrelatorVsPtSumDiffHist[sd]->SetBinContent(b,corrected);
2250     f3pCorrelatorVsPtSumDiffHist[sd]->SetBinError(b,correctedErr);
2251    }
2252   } // end of for(Int_t b=1;b<=fnBinsPt;b++)
2253
2254   // [(eta1+eta2)/2,|eta1-eta2|]
2255   // looping over all bins and calculating reduced correlations:
2256   for(Int_t b=1;b<=fnBinsEta;b++)
2257   {
2258    Double_t measured = f3pCorrelatorVsEtaSumDiffPro[sd]->GetBinContent(b);
2259    Double_t measuredErr = f3pCorrelatorVsEtaSumDiffPro[sd]->GetBinError(b);
2260    Double_t corrected = 0.; // 3-p correlator corrected for detector effects
2261    Double_t correctedErr = measuredErr; // to be improved - propagate error also for non-isotropic terms
2262    // non-isotropic terms:
2263    Double_t cosPsiPOI1 = fNonIsotropicTermsVsEtaSumDiffPro[sd][0]->GetBinContent(b); // <<cos(#psi_{POI_1})>>
2264    Double_t sinPsiPOI1 = fNonIsotropicTermsVsEtaSumDiffPro[sd][1]->GetBinContent(b); // <<sin(#psi_{POI_1})>>
2265    Double_t cosPsiPOI2 = fNonIsotropicTermsVsEtaSumDiffPro[sd][2]->GetBinContent(b); // <<cos(#psi_{POI_2})>>
2266    Double_t sinPsiPOI2 = fNonIsotropicTermsVsEtaSumDiffPro[sd][3]->GetBinContent(b); // <<sin(#psi_{POI_2})>>
2267    Double_t cosPsiPOI1m2PhiRP = fNonIsotropicTermsVsEtaSumDiffPro[sd][4]->GetBinContent(b); // <<cos(#psi_{POI_1}-2*phi_{RP})>>
2268    Double_t sinPsiPOI1m2PhiRP = fNonIsotropicTermsVsEtaSumDiffPro[sd][5]->GetBinContent(b); // <<sin(#psi_{POI_1}-2*phi_{RP})>>
2269    Double_t cosPsiPOI2m2PhiRP = fNonIsotropicTermsVsEtaSumDiffPro[sd][6]->GetBinContent(b); // <<cos(#psi_{POI_2}-2*phi_{RP})>>
2270    Double_t sinPsiPOI2m2PhiRP = fNonIsotropicTermsVsEtaSumDiffPro[sd][7]->GetBinContent(b); // <<sin(#psi_{POI_2}-2*phi_{RP})>>
2271    Double_t cosPsiPOI1pPsiPOI2 = fNonIsotropicTermsVsEtaSumDiffPro[sd][8]->GetBinContent(b); // <<cos(#psi_{POI_1}+#psi_{POI_2})>>
2272    Double_t sinPsiPOI1pPsiPOI2 = fNonIsotropicTermsVsEtaSumDiffPro[sd][9]->GetBinContent(b); // <<sin(#psi_{POI_1}+#psi_{POI_2})>>
2273    Double_t cos2PhiRP = fNonIsotropicTermsPro->GetBinContent(3); // <<cos(2n*(phi_{RP}))>>
2274    Double_t sin2PhiRP = fNonIsotropicTermsPro->GetBinContent(4); // <<sin(2n*(phi_{RP}))>>
2275    // apply correction:
2276    corrected = measured
2277              - (cosPsiPOI1*cosPsiPOI2m2PhiRP-sinPsiPOI1*sinPsiPOI2m2PhiRP
2278              + cosPsiPOI2*cosPsiPOI1m2PhiRP-sinPsiPOI2*sinPsiPOI1m2PhiRP
2279              + cos2PhiRP*cosPsiPOI1pPsiPOI2+sin2PhiRP*sinPsiPOI1pPsiPOI2)
2280              + 2.*cos2PhiRP*(cosPsiPOI1*cosPsiPOI2-sinPsiPOI1*sinPsiPOI2)
2281              + 2.*sin2PhiRP*(cosPsiPOI1*sinPsiPOI2+sinPsiPOI1*cosPsiPOI2);
2282    if(fCorrectForDetectorEffects)
2283    {
2284     f3pCorrelatorVsEtaSumDiffHist[sd]->SetBinContent(b,corrected);
2285     f3pCorrelatorVsEtaSumDiffHist[sd]->SetBinError(b,correctedErr);
2286    }
2287   } // end of for(Int_t b=1;b<=fnBinsEta;b++)
2288  } // end of for(Int_t sd=0;sd<2;sd++)
2289
2290 } // end of AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffects()
2291
2292 //================================================================================================================
2293
2294 void AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffectsVsM()
2295 {
2296  // Correct measured 3-p correlator cos[n(phi1+phi2-2phi3)] vs M for detector effects.
2297
2298  for(Int_t b=1;b<=fNoOfMultipicityBins+2;b++)
2299  {
2300   Double_t measured3pCorrelator = f3pCorrelatorVsMPro->GetBinContent(b); // biased by detector effects
2301   Double_t corrected3pCorrelator = 0.; // corrected for detector effects
2302   Double_t nonIsotropicTerms[10] = {0.}; // there are 10 distinct non-isotropic terms
2303   for(Int_t nit=0;nit<10;nit++)
2304   {
2305    nonIsotropicTerms[nit] = fNonIsotropicTermsVsMPro->GetBinContent(fNonIsotropicTermsVsMPro->GetBin(nit+1,b));
2306   }
2307   // Calculate corrected 3-p correlator:
2308   corrected3pCorrelator = measured3pCorrelator
2309                         - nonIsotropicTerms[2]*nonIsotropicTerms[4]
2310                         - nonIsotropicTerms[3]*nonIsotropicTerms[5]
2311                         - 2.*nonIsotropicTerms[0]*nonIsotropicTerms[6]
2312                         - 2.*nonIsotropicTerms[1]*nonIsotropicTerms[7]
2313                         + 2.*nonIsotropicTerms[2]*(pow(nonIsotropicTerms[0],2.)-pow(nonIsotropicTerms[1],2.))
2314                         + 4.*nonIsotropicTerms[3]*nonIsotropicTerms[0]*nonIsotropicTerms[1];
2315   // Store corrected correlator:
2316   if(fCorrectForDetectorEffects)
2317   {
2318    f3pCorrelatorVsMHist->SetBinContent(b,corrected3pCorrelator);
2319    f3pCorrelatorVsMHist->SetBinError(b,f3pCorrelatorVsMPro->GetBinError(b)); // to be improved (propagate error for non-isotropic terms)
2320   }
2321   // Quantify bias from detector inefficiences to 3-p correlator. Remark: Bias is quantified as a
2322   // ratio between corrected and measured 3-p correlator:
2323   //              bias = corrected/measured
2324   // This bias is stored in histogram fDetectorBias.
2325   Double_t bias = 0.;
2326   if(measured3pCorrelator)
2327   {
2328    bias = corrected3pCorrelator/measured3pCorrelator;
2329    fDetectorBiasVsMHist->SetBinContent(b,bias);
2330   }
2331  } // end of for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
2332
2333 } // end of AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffectsVsM()
2334
2335 //================================================================================================================
2336
2337 void AliFlowAnalysisWithMixedHarmonics::ResetEventByEventQuantities()
2338 {
2339  // Reset all event-by-event quantities.
2340
2341  fReQnk->Zero();
2342  fImQnk->Zero();
2343  fSpk->Zero();
2344
2345  if(!fEvaluateDifferential3pCorrelator){return;}
2346  for(Int_t sd=0;sd<2;sd++)
2347  {
2348   fRePEBE[sd]->Reset();
2349   fImPEBE[sd]->Reset();
2350   fReEtaEBE[sd]->Reset();
2351   fImEtaEBE[sd]->Reset();
2352  }
2353  for(Int_t fs=0;fs<2;fs++)
2354  {
2355   for(Int_t sd=0;sd<2;sd++)
2356   {
2357    fOverlapEBE[fs][sd]->Reset();
2358    fOverlapEBE2[fs][sd]->Reset();
2359   }
2360  }
2361  for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
2362  {
2363   for(Int_t ao=0;ao<2;ao++) // all/overlap
2364   {
2365    for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
2366    {
2367     fReNITEBE[p12][ao][pe]->Reset();
2368     fImNITEBE[p12][ao][pe]->Reset();
2369    } // end of for(Int_t pe=0;pe<4;pe++) // [(p1+p2)/2,|p1-p2|,(eta1+eta2)/2,|eta1-eta2|]
2370   }
2371  } // end of for(Int_t p12=0;p12<2;p12++) // 1st/2nd POI
2372
2373 } // end of void AliFlowAnalysisWithMixedHarmonics::ResetEventByEventQuantities()
2374
2375 //================================================================================================================
2376
2377 void AliFlowAnalysisWithMixedHarmonics::Calculate3pCorrelator()
2378 {
2379  // Calculate 3-p azimuthal correlator cos[n(phi1+phi2-2phi3)] in terms of Q_{n,k} and S_{p,k}.
2380
2381  // a) Calculate 3-p correlator without using particle weights;
2382  // b) Calculate 3-p correlator with using particle weights.
2383
2384  // a) Calculate 3-p correlator without using particle weights:
2385  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2386  {
2387   // Multiplicity (number of RPs):
2388   Double_t dMult = (*fSpk)(0,0);
2389   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n and 2n:
2390   Double_t dReQ1n = (*fReQnk)(0,0);
2391   Double_t dReQ2n = (*fReQnk)(1,0);
2392   Double_t dImQ1n = (*fImQnk)(0,0);
2393   Double_t dImQ2n = (*fImQnk)(1,0);
2394   // 3-particle azimuthal correlator <cos(n*(phi1+phi2-2phi3))>:
2395   Double_t three1n1n2n = (pow(dReQ1n,2.)*dReQ2n + 2.*dReQ1n*dImQ1n*dImQ2n - pow(dImQ1n,2.)*dReQ2n
2396                        - 2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
2397                        - (pow(dReQ2n,2.)+pow(dImQ2n,2.))+2.*dMult)
2398                        / (dMult*(dMult-1.)*(dMult-2.));
2399
2400   // Fill all-events profile:
2401   f3pCorrelatorPro->Fill(0.5,three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
2402
2403   // 3-particle azimuthal correlator <cos(n*(phi1+phi2-2phi3))> vs multiplicity:
2404   if(fCalculateVsM)
2405   {
2406    if(dMult<fMinMultiplicity)
2407    {
2408     f3pCorrelatorVsMPro->Fill(0.5,three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
2409    } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
2410      {
2411       f3pCorrelatorVsMPro->Fill(0.5+fNoOfMultipicityBins+1,three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
2412      } else
2413        {
2414         f3pCorrelatorVsMPro->Fill(1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
2415        }
2416   } // end of if(fCalculateVsM)
2417
2418  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2419
2420  // b) Calculate 3-p correlator with using particle weights:
2421  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2422  {
2423   // ...
2424  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2425
2426 } // end of void AliFlowAnalysisWithMixedHarmonics::Calculate3pCorrelator()
2427
2428 //================================================================================================================
2429
2430 void AliFlowAnalysisWithMixedHarmonics::Calculate5pCorrelator()
2431 {
2432  // Calculate 5-p azimuthal correlator cos[n(2phi1+2phi2+2phi3-3phi4-3phi5)] in terms of Q_{n,k} and S_{p,k}.
2433
2434  // a) Calculate 5-p correlator without using particle weights;
2435  // b) Calculate 5-p correlator with using particle weights.
2436
2437  // a) Calculate 5-p correlator without using particle weights:
2438  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2439  {
2440   // Multiplicity (number of RPs):
2441   Double_t dMult = (*fSpk)(0,0);
2442   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n,2n,...,6n:
2443   Double_t dReQ1n = (*fReQnk)(0,0);
2444   Double_t dReQ2n = (*fReQnk)(1,0);
2445   Double_t dReQ3n = (*fReQnk)(2,0);
2446   Double_t dReQ4n = (*fReQnk)(3,0);
2447   //Double_t dReQ5n = (*fReQnk)(4,0); // not needed
2448   Double_t dReQ6n = (*fReQnk)(5,0);
2449   Double_t dImQ1n = (*fImQnk)(0,0);
2450   Double_t dImQ2n = (*fImQnk)(1,0);
2451   Double_t dImQ3n = (*fImQnk)(2,0);
2452   Double_t dImQ4n = (*fImQnk)(3,0);
2453   //Double_t dImQ5n = (*fImQnk)(4,0); // not needed
2454   Double_t dImQ6n = (*fImQnk)(5,0);
2455
2456   // 5-particle azimuthal correlator:
2457   Double_t five2n2n2n3n3n = 0.; // <cos[n(2phi1+2phi2+2phi3-3phi4-3phi5)]>
2458   Double_t reQ2nQ2nQ2nQ3nstarQ3nstar = pow(dReQ2n,3.)*pow(dReQ3n,2.)
2459                                      - 3.*dReQ2n*pow(dReQ3n,2.)*pow(dImQ2n,2.)
2460                                      + 6.*pow(dReQ2n,2.)*dReQ3n*dImQ2n*dImQ3n
2461                                      - 2.*dReQ3n*pow(dImQ2n,3.)*dImQ3n-pow(dReQ2n,3.)*pow(dImQ3n,2.)
2462                                      + 3.*dReQ2n*pow(dImQ2n,2.)*pow(dImQ3n,2.);
2463   Double_t reQ2nQ2nQ2nQ6nstar = dReQ6n*pow(dReQ2n,3)-3.*dReQ2n*dReQ6n*pow(dImQ2n,2)
2464                                + 3.*dImQ2n*dImQ6n*pow(dReQ2n,2)-dImQ6n*pow(dImQ2n,3);
2465   Double_t reQ4nQ2nQ3nstarQ3nstar = (dReQ4n*dReQ2n-dImQ4n*dImQ2n)*(dReQ3n*dReQ3n-dImQ3n*dImQ3n)
2466                                   + 2.*(dReQ4n*dImQ2n+dImQ4n*dReQ2n)*dReQ3n*dImQ3n;
2467   Double_t reQ2nQ2nQ1nstarQ3nstar = (pow(dReQ2n,2.)-pow(dImQ2n,2.))*(dReQ3n*dReQ1n-dImQ3n*dImQ1n)
2468                                   + 2.*dReQ2n*dImQ2n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
2469   Double_t reQ6nQ3nstarQ3nstar = pow(dReQ3n,2.)*dReQ6n + 2.*dReQ3n*dImQ3n*dImQ6n
2470                                - pow(dImQ3n,2.)*dReQ6n;
2471   Double_t reQ4nQ2nQ6nstar = dReQ6n*dReQ4n*dReQ2n-dReQ6n*dImQ4n*dImQ2n+dImQ6n*dReQ4n*dImQ2n
2472                            + dImQ6n*dImQ4n*dReQ2n;
2473   Double_t reQ4nQ1nstarQ3nstar = dReQ4n*(dReQ3n*dReQ1n-dImQ3n*dImQ1n)+dImQ4n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
2474   Double_t reQ2nQ2nQ4nstar = pow(dReQ2n,2.)*dReQ4n+2.*dReQ2n*dImQ2n*dImQ4n-pow(dImQ2n,2.)*dReQ4n;
2475   Double_t reQ2nQ1nQ3nstar = dReQ3n*dReQ2n*dReQ1n-dReQ3n*dImQ2n*dImQ1n+dImQ3n*dReQ2n*dImQ1n
2476                            + dImQ3n*dImQ2n*dReQ1n;
2477   Double_t reQ2nQ1nstarQ1nstar = pow(dReQ1n,2.)*dReQ2n + 2.*dReQ1n*dImQ1n*dImQ2n - pow(dImQ1n,2.)*dReQ2n;
2478   // Analytic expression for 5-particle azimuthal correlator:
2479   five2n2n2n3n3n = (reQ2nQ2nQ2nQ3nstarQ3nstar-reQ2nQ2nQ2nQ6nstar-3.*reQ4nQ2nQ3nstarQ3nstar
2480                  - 6.*reQ2nQ2nQ1nstarQ3nstar+2.*reQ6nQ3nstarQ3nstar+3.*reQ4nQ2nQ6nstar
2481                  + 6.*reQ4nQ1nstarQ3nstar+6.*reQ2nQ2nQ4nstar
2482                  + 12.*reQ2nQ1nQ3nstar+6.*reQ2nQ1nstarQ1nstar
2483                  - 2.*((pow(dReQ6n,2.)+pow(dImQ6n,2.))
2484                  + 3.*(pow(dReQ4n,2.)+pow(dImQ4n,2.))
2485                  + 6.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2486                  + 9.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2487                  + 6.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))-12.*dMult))
2488                  /(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2489   // Fill all-events profile:
2490   f5pCorrelatorPro->Fill(0.5,five2n2n2n3n3n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2491  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2492
2493  // b) Calculate 5-p correlator with using particle weights:
2494  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2495  {
2496   // ...
2497  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2498
2499 } // end of void AliFlowAnalysisWithMixedHarmonics::Calculate5pCorrelator()
2500
2501 //================================================================================================================
2502
2503 void AliFlowAnalysisWithMixedHarmonics::CalculateNonIsotropicTerms()
2504 {
2505  // Calculate non-isotropic terms which appear in the decomposition of 3-p correlator <cos[n(phi1+phi2-2phi3)]>.
2506
2507  // a) Calculate without using particle weights;
2508  // b) Calculate using particle weights;
2509
2510  // For detector with uniform acceptance all these terms vanish. These non-isotropic terms are stored in fNonIsotropicTermsPro.
2511  // Binning of fNonIsotropicTermsPro is organized as follows:
2512  //  1st bin: <<cos(n*phi1)>>
2513  //  2nd bin: <<sin(n*phi1)>>
2514  //  3rd bin: <<cos(2n*phi1)>>
2515  //  4th bin: <<sin(2n*phi1)>>
2516  //  5th bin: <<cos(n*(phi1+phi2)>>
2517  //  6th bin: <<sin(n*(phi1+phi2)>>
2518  //  7th bin: <<cos(n*(2phi1-phi2)>>
2519  //  8th bin: <<sin(n*(2phi1-phi2)>>
2520  //  9th bin: <<cos(n*(phi1-phi2-phi3)>> // not needed
2521  // 10th bin: <<sin(n*(phi1-phi2-phi3)>> // not needed
2522
2523  // a) Calculate without using particle weights:
2524  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2525  {
2526   // Multiplicity (number of RPs):
2527   Double_t dMult = (*fSpk)(0,0);
2528   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n and 2n:
2529   Double_t dReQ1n = (*fReQnk)(0,0);
2530   Double_t dReQ2n = (*fReQnk)(1,0);
2531   Double_t dImQ1n = (*fImQnk)(0,0);
2532   Double_t dImQ2n = (*fImQnk)(1,0);
2533   // 1-particle terms:
2534   Double_t cosP1n = 0.; // <cos(n*(phi1))>
2535   Double_t sinP1n = 0.; // <sin(n*(phi1))>
2536   Double_t cosP2n = 0.; // <cos(2n*(phi1))>
2537   Double_t sinP2n = 0.; // <sin(2n*(phi1))>
2538   if(dMult>0)
2539   {
2540    cosP1n = dReQ1n/dMult;
2541    sinP1n = dImQ1n/dMult;
2542    cosP2n = dReQ2n/dMult;
2543    sinP2n = dImQ2n/dMult;
2544    // All-event avarages:
2545    fNonIsotropicTermsPro->Fill(0.5,cosP1n,dMult); // <<cos(n*(phi1))>>
2546    fNonIsotropicTermsPro->Fill(1.5,sinP1n,dMult); // <<sin(n*(phi1))>>
2547    fNonIsotropicTermsPro->Fill(2.5,cosP2n,dMult); // <<cos(2n*(phi1))>>
2548    fNonIsotropicTermsPro->Fill(3.5,sinP2n,dMult); // <<sin(2n*(phi1))>>
2549    // All-event avarages vs M:
2550    if(fCalculateVsM)
2551    {
2552     if(dMult<fMinMultiplicity)
2553     {
2554      fNonIsotropicTermsVsMPro->Fill(0.5,0.5,cosP1n,dMult); // <<cos(n*(phi1))>>
2555      fNonIsotropicTermsVsMPro->Fill(1.5,0.5,sinP1n,dMult); // <<sin(n*(phi1))>>
2556      fNonIsotropicTermsVsMPro->Fill(2.5,0.5,cosP2n,dMult); // <<cos(2n*(phi1))>>
2557      fNonIsotropicTermsVsMPro->Fill(3.5,0.5,sinP2n,dMult); // <<sin(2n*(phi1))>>
2558     } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
2559       {
2560        fNonIsotropicTermsVsMPro->Fill(0.5,0.5+fNoOfMultipicityBins+1,cosP1n,dMult); // <<cos(n*(phi1))>>
2561        fNonIsotropicTermsVsMPro->Fill(1.5,0.5+fNoOfMultipicityBins+1,sinP1n,dMult); // <<sin(n*(phi1))>>
2562        fNonIsotropicTermsVsMPro->Fill(2.5,0.5+fNoOfMultipicityBins+1,cosP2n,dMult); // <<cos(2n*(phi1))>>
2563        fNonIsotropicTermsVsMPro->Fill(3.5,0.5+fNoOfMultipicityBins+1,sinP2n,dMult); // <<sin(2n*(phi1))>>
2564       } else
2565         {
2566          fNonIsotropicTermsVsMPro->Fill(0.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP1n,dMult); // <<cos(n*(phi1))>>
2567          fNonIsotropicTermsVsMPro->Fill(1.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP1n,dMult); // <<sin(n*(phi1))>>
2568          fNonIsotropicTermsVsMPro->Fill(2.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP2n,dMult); // <<cos(2n*(phi1))>>
2569          fNonIsotropicTermsVsMPro->Fill(3.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP2n,dMult); // <<sin(2n*(phi1))>>
2570         }
2571    } // end of if(fCalculateVsM)
2572   } // end of if(dMult>0)
2573   // 2-particle terms:
2574   Double_t cosP1nP1n = 0.; // <cos(n*(phi1+phi2))>
2575   Double_t sinP1nP1n = 0.; // <sin(n*(phi1+phi2))>
2576   Double_t cosP2nM1n = 0.; // <cos(n*(2phi1-phi2))>
2577   Double_t sinP2nM1n = 0.; // <sin(n*(2phi1-phi2))>
2578   if(dMult>1)
2579   {
2580    cosP1nP1n = (pow(dReQ1n,2)-pow(dImQ1n,2)-dReQ2n)/(dMult*(dMult-1));
2581    sinP1nP1n = (2.*dReQ1n*dImQ1n-dImQ2n)/(dMult*(dMult-1));
2582    cosP2nM1n = (dReQ2n*dReQ1n+dImQ2n*dImQ1n-dReQ1n)/(dMult*(dMult-1));
2583    sinP2nM1n = (dImQ2n*dReQ1n-dReQ2n*dImQ1n-dImQ1n)/(dMult*(dMult-1));
2584    // All-event avarages:
2585    fNonIsotropicTermsPro->Fill(4.5,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
2586    fNonIsotropicTermsPro->Fill(5.5,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
2587    fNonIsotropicTermsPro->Fill(6.5,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
2588    fNonIsotropicTermsPro->Fill(7.5,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
2589    // All-event avarages vs M:
2590    if(fCalculateVsM)
2591    {
2592     if(dMult<fMinMultiplicity)
2593     {
2594      fNonIsotropicTermsVsMPro->Fill(4.5,0.5,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
2595      fNonIsotropicTermsVsMPro->Fill(5.5,0.5,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
2596      fNonIsotropicTermsVsMPro->Fill(6.5,0.5,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
2597      fNonIsotropicTermsVsMPro->Fill(7.5,0.5,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
2598     } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
2599       {
2600        fNonIsotropicTermsVsMPro->Fill(4.5,0.5+fNoOfMultipicityBins+1,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
2601        fNonIsotropicTermsVsMPro->Fill(5.5,0.5+fNoOfMultipicityBins+1,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
2602        fNonIsotropicTermsVsMPro->Fill(6.5,0.5+fNoOfMultipicityBins+1,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
2603        fNonIsotropicTermsVsMPro->Fill(7.5,0.5+fNoOfMultipicityBins+1,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
2604       } else
2605         {
2606          fNonIsotropicTermsVsMPro->Fill(4.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
2607          fNonIsotropicTermsVsMPro->Fill(5.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
2608          fNonIsotropicTermsVsMPro->Fill(6.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
2609          fNonIsotropicTermsVsMPro->Fill(7.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
2610         }
2611    } // end of if(fCalculateVsM)
2612   } // end of if(dMult>1)
2613   // 3-particle: correct and ready but not needed, hence commented out.
2614   /*
2615   Double_t cosP1nM1nM1n = 0.; // <cos(n*(phi1-phi2-phi3))>
2616   Double_t sinP1nM1nM1n = 0.; // <sin(n*(phi1-phi2-phi3))>
2617   if(dMult>2)
2618   {
2619    cosP1nM1nM1n = (dReQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))-dReQ1n*dReQ2n-dImQ1n*dImQ2n-2.*(dMult-1)*dReQ1n)
2620                 / (dMult*(dMult-1)*(dMult-2));
2621    sinP1nM1nM1n = (-dImQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))+dReQ1n*dImQ2n-dImQ1n*dReQ2n+2.*(dMult-1)*dImQ1n)
2622                 / (dMult*(dMult-1)*(dMult-2));
2623    // All-events avarages:
2624    fNonIsotropicTermsPro->Fill(8.5,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
2625    fNonIsotropicTermsPro->Fill(9.5,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
2626    // All-events avarages vs M:
2627    if(fCalculateVsM)
2628    {
2629     if(dMult<fMinMultiplicity)
2630     {
2631      fNonIsotropicTermsVsMPro->Fill(8.5,0.5,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
2632      fNonIsotropicTermsVsMPro->Fill(9.5,0.5,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
2633     } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
2634       {
2635        fNonIsotropicTermsVsMPro->Fill(8.5,0.5+fNoOfMultipicityBins+1,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
2636        fNonIsotropicTermsVsMPro->Fill(9.5,0.5+fNoOfMultipicityBins+1,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
2637       } else
2638         {
2639          fNonIsotropicTermsVsMPro->Fill(8.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),
2640                                         cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
2641          fNonIsotropicTermsVsMPro->Fill(9.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),
2642                                         sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
2643         }
2644    } // end of if(fCalculateVsM)
2645   } // end of if(dMult>2)
2646   */
2647  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2648
2649  // b) Calculate using particle weights:
2650  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2651  {
2652   // ...
2653  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2654
2655 } // end of void AliFlowAnalysisWithMixedHarmonics::CalculateNonIsotropicTerms()
2656
2657 //================================================================================================================
2658
2659 void AliFlowAnalysisWithMixedHarmonics::CalculateDifferential3pCorrelator(Double_t &gIntegratedValue)
2660 {
2661  // Calculate differential 3-p azimuthal correlator cos[n(psi1+psi2-2phi3)] in terms of Q_{2n}, p_{n}, q1_{n} and q2_{n}.
2662
2663  // a) Calculate differential 3-p correlator without using particle weights;
2664  // b) Calculate differential 3-p correlator with using particle weights;
2665  // c) Calculate non-isotropic terms for 3-p correlator.
2666
2667  // a) Calculate differential 3-p correlator without using particle weights:
2668  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2669  {
2670    Int_t iBinCounter = 0;
2671    Double_t gSumBinContentTimesWeight = 0., gSumWeight = 0.;
2672    Double_t gSumBinContentTimesWeightSquared = 0.;
2673
2674   // Multiplicity (number of RPs):
2675   Double_t dMult = (*fSpk)(0,0);
2676   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonic 2n:
2677   Double_t dReQ2n = (*fReQnk)(1,0);
2678   Double_t dImQ2n = (*fImQnk)(1,0);
2679   for(Int_t sd=0;sd<2;sd++)
2680   {
2681    // [(p1+p2)/2,|p1-p2|]
2682    // looping over all bins and calculating reduced correlations:
2683    for(Int_t b=1;b<=fnBinsPt;b++)
2684    {
2685     // real and imaginary parts of p_{n}:
2686     Double_t p1nRe = fRePEBE[sd]->GetBinContent(b)*fRePEBE[sd]->GetBinEntries(b);
2687     Double_t p1nIm = fImPEBE[sd]->GetBinContent(b)*fImPEBE[sd]->GetBinEntries(b);
2688     // overlap 1: to be improved (terminology)
2689     Double_t overlap1 = fOverlapEBE[0][sd]->GetBinContent(b)*fOverlapEBE[0][sd]->GetBinEntries(b);
2690     // overlap 2: to be improved (terminology)
2691     Double_t overlap2 = fOverlapEBE[1][sd]->GetBinContent(b)*fOverlapEBE[1][sd]->GetBinEntries(b);
2692     // number of pairs of POIs in particular (p1+p2)/2 or |p1-p2| bin:
2693     Double_t mp = fRePEBE[sd]->GetBinEntries(b);
2694     // number of pairs of POI1/RP and POI2 in particular (p1+p2)/2 or |p1-p2| bin:
2695     Double_t mOverlap1 = fOverlapEBE[0][sd]->GetBinEntries(b);
2696     // number of pairs of POI2/RP and POI1 in particular (p1+p2)/2 or |p1-p2| bin:
2697     Double_t mOverlap2 = fOverlapEBE[1][sd]->GetBinEntries(b);
2698     // e-b-e weight for cos[n(psi1+psi2-2phi3)]:
2699     Double_t weight = mp*dMult-mOverlap1-mOverlap2;
2700
2701     Double_t cosP2nphi1M1npsi2M1npsi2 = 0; // cos[n(psi1+psi2-2phi3)]
2702     if(weight>0.)
2703     {
2704      cosP2nphi1M1npsi2M1npsi2 = (p1nRe*dReQ2n+p1nIm*dImQ2n-overlap1-overlap2)/(weight);
2705     }
2706     f3pCorrelatorVsPtSumDiffPro[sd]->Fill(fPtMin+(b-1)*fPtBinWidth,cosP2nphi1M1npsi2M1npsi2,weight);
2707     if(sd == 0) {
2708       iBinCounter += 1;
2709
2710       gSumBinContentTimesWeight += f3pCorrelatorVsPtSumDiffPro[sd]->GetBinContent(b)*f3pCorrelatorVsPtSumDiffPro[sd]->GetBinEntries(b);
2711       gSumWeight += f3pCorrelatorVsPtSumDiffPro[sd]->GetBinEntries(b);
2712       gSumBinContentTimesWeightSquared += TMath::Power(gSumBinContentTimesWeight,2);
2713     }
2714
2715     // non-isotropic terms, 1st POI:
2716     Double_t p1nRePOI1 = fReNITEBE[0][0][sd]->GetBinContent(b)*fReNITEBE[0][0][sd]->GetBinEntries(b);
2717     Double_t p1nImPOI1 = fImNITEBE[0][0][sd]->GetBinContent(b)*fImNITEBE[0][0][sd]->GetBinEntries(b);
2718     Double_t mpPOI1 = fReNITEBE[0][0][sd]->GetBinEntries(b);
2719     Double_t q1nRePOI1 = fReNITEBE[0][1][sd]->GetBinContent(b)*fReNITEBE[0][1][sd]->GetBinEntries(b);
2720     Double_t q1nImPOI1 = fImNITEBE[0][1][sd]->GetBinContent(b)*fImNITEBE[0][1][sd]->GetBinEntries(b);
2721     Double_t mqPOI1 = fReNITEBE[0][1][sd]->GetBinEntries(b);
2722     // non-isotropic terms, 2nd POI:
2723     Double_t p1nRePOI2 = fReNITEBE[1][0][sd]->GetBinContent(b)*fReNITEBE[1][0][sd]->GetBinEntries(b);
2724     Double_t p1nImPOI2 = fImNITEBE[1][0][sd]->GetBinContent(b)*fImNITEBE[1][0][sd]->GetBinEntries(b);
2725     Double_t mpPOI2 = fReNITEBE[1][0][sd]->GetBinEntries(b);
2726     Double_t q1nRePOI2 = fReNITEBE[1][1][sd]->GetBinContent(b)*fReNITEBE[1][1][sd]->GetBinEntries(b);
2727     Double_t q1nImPOI2 = fImNITEBE[1][1][sd]->GetBinContent(b)*fImNITEBE[1][1][sd]->GetBinEntries(b);
2728     Double_t mqPOI2 = fReNITEBE[1][1][sd]->GetBinEntries(b);
2729     // Fill all-event profiles:
2730     if(weight>0. && mpPOI1>0.)
2731     {
2732      fNonIsotropicTermsVsPtSumDiffPro[sd][0]->Fill(fPtMin+(b-1)*fPtBinWidth,p1nRePOI1/mpPOI1,mpPOI1); // <<cos(#psi_{POI_1})>>
2733      fNonIsotropicTermsVsPtSumDiffPro[sd][1]->Fill(fPtMin+(b-1)*fPtBinWidth,p1nImPOI1/mpPOI1,mpPOI1); // <<sin(#psi_{POI_1})>>
2734     }
2735     if(weight>0. && mpPOI2>0.)
2736     {
2737      fNonIsotropicTermsVsPtSumDiffPro[sd][2]->Fill(fPtMin+(b-1)*fPtBinWidth,p1nRePOI2/mpPOI2,mpPOI2); // <<cos(#psi_{POI_2})>>
2738      fNonIsotropicTermsVsPtSumDiffPro[sd][3]->Fill(fPtMin+(b-1)*fPtBinWidth,p1nImPOI2/mpPOI2,mpPOI2); // <<sin(#psi_{POI_2})>>
2739     }
2740     if(weight>0. && mpPOI1*dMult-mqPOI1>0.)
2741     {
2742      fNonIsotropicTermsVsPtSumDiffPro[sd][4]->Fill(fPtMin+(b-1)*fPtBinWidth,
2743       (p1nRePOI1*dReQ2n+p1nImPOI1*dImQ2n-q1nRePOI1)/(mpPOI1*dMult-mqPOI1),mpPOI1*dMult-mqPOI1); // <<cos(#psi_{POI_1}-2*phi)>>
2744      fNonIsotropicTermsVsPtSumDiffPro[sd][5]->Fill(fPtMin+(b-1)*fPtBinWidth,
2745       (p1nImPOI1*dReQ2n-p1nRePOI1*dImQ2n+q1nImPOI1)/(mpPOI1*dMult-mqPOI1),mpPOI1*dMult-mqPOI1); // <<sin(#psi_{POI_1}-2*phi)>>
2746     }
2747     if(weight>0. && mpPOI2*dMult-mqPOI2>0.)
2748     {
2749      fNonIsotropicTermsVsPtSumDiffPro[sd][6]->Fill(fPtMin+(b-1)*fPtBinWidth,
2750       (p1nRePOI2*dReQ2n+p1nImPOI2*dImQ2n-q1nRePOI2)/(mpPOI2*dMult-mqPOI2),mpPOI2*dMult-mqPOI2); // <<cos(#psi_{POI_2}-2*phi)>>
2751      fNonIsotropicTermsVsPtSumDiffPro[sd][7]->Fill(fPtMin+(b-1)*fPtBinWidth,
2752       (p1nImPOI2*dReQ2n-p1nRePOI2*dImQ2n+q1nImPOI2)/(mpPOI2*dMult-mqPOI2),mpPOI2*dMult-mqPOI2); // <<sin(#psi_{POI_2}-2*phi)>>
2753     }
2754     if(weight>0. && mp>0.)
2755     {
2756      fNonIsotropicTermsVsPtSumDiffPro[sd][8]->Fill(fPtMin+(b-1)*fPtBinWidth,p1nRe/mp,mp); // <<cos(#psi_{POI_1}+#psi_{POI_2})>>
2757      fNonIsotropicTermsVsPtSumDiffPro[sd][9]->Fill(fPtMin+(b-1)*fPtBinWidth,p1nIm/mp,mp); // <<sin(#psi_{POI_1}+#psi_{POI_2})>>
2758     }
2759    } // end of for(Int_t b=1;b<=fnBinsPt;b++)
2760
2761    // [(eta1+eta2)/2,|eta1-eta2|]
2762    // looping over all bins and calculating reduced correlations:
2763    for(Int_t k=1;k<=fnBinsEta;k++)
2764    {
2765     // real and imaginary parts of p_{n}:
2766     Double_t p1nRe = fReEtaEBE[sd]->GetBinContent(k)*fReEtaEBE[sd]->GetBinEntries(k);
2767     Double_t p1nIm = fImEtaEBE[sd]->GetBinContent(k)*fImEtaEBE[sd]->GetBinEntries(k);
2768     // overlap 1: to be improved (terminology)
2769     Double_t overlap1 = fOverlapEBE2[0][sd]->GetBinContent(k)*fOverlapEBE2[0][sd]->GetBinEntries(k);
2770     // overlap 2: to be improved (terminology)
2771     Double_t overlap2 = fOverlapEBE2[1][sd]->GetBinContent(k)*fOverlapEBE2[1][sd]->GetBinEntries(k);
2772     // number of pairs of POIs in particular (eta1+eta2)/2 or |eta1-eta2| bin:
2773     Double_t mp = fReEtaEBE[sd]->GetBinEntries(k);
2774     // number of pairs of POI1/RP and POI2 in particular (eta1+eta2)/2 or |eta1-eta2| bin:
2775     Double_t mOverlap1 = fOverlapEBE2[0][sd]->GetBinEntries(k);
2776     // number of pairs of POI2/RP and POI1 in particular (eta1+eta2)/2 or |eta1-eta2| bin:
2777     Double_t mOverlap2 = fOverlapEBE2[1][sd]->GetBinEntries(k);
2778     // e-b-e weight for cos[n(psi1+psi2-2phi3)]:
2779     Double_t weight = mp*dMult-mOverlap1-mOverlap2;
2780
2781     Double_t cosP2nphi1M1npsi2M1npsi2 = 0; // cos[n(psi1+psi2-2phi3)]
2782     if(weight>0.)
2783     {
2784      cosP2nphi1M1npsi2M1npsi2 = (p1nRe*dReQ2n+p1nIm*dImQ2n-overlap1-overlap2)/(weight);
2785     }
2786     f3pCorrelatorVsEtaSumDiffPro[sd]->Fill(fEtaMin+(k-1)*fEtaBinWidth,cosP2nphi1M1npsi2M1npsi2,weight);
2787
2788     // non-isotropic terms, 1st POI:
2789     Double_t p1nRePOI1 = fReNITEBE[0][0][sd+2]->GetBinContent(k)*fReNITEBE[0][0][sd+2]->GetBinEntries(k);
2790     Double_t p1nImPOI1 = fImNITEBE[0][0][sd+2]->GetBinContent(k)*fImNITEBE[0][0][sd+2]->GetBinEntries(k);
2791     Double_t mpPOI1 = fReNITEBE[0][0][sd+2]->GetBinEntries(k);
2792     Double_t q1nRePOI1 = fReNITEBE[0][1][sd+2]->GetBinContent(k)*fReNITEBE[0][1][sd+2]->GetBinEntries(k);
2793     Double_t q1nImPOI1 = fImNITEBE[0][1][sd+2]->GetBinContent(k)*fImNITEBE[0][1][sd+2]->GetBinEntries(k);
2794     Double_t mqPOI1 = fReNITEBE[0][1][sd+2]->GetBinEntries(k);
2795     // non-isotropic terms, 2nd POI:
2796     Double_t p1nRePOI2 = fReNITEBE[1][0][sd+2]->GetBinContent(k)*fReNITEBE[1][0][sd+2]->GetBinEntries(k);
2797     Double_t p1nImPOI2 = fImNITEBE[1][0][sd+2]->GetBinContent(k)*fImNITEBE[1][0][sd+2]->GetBinEntries(k);
2798     Double_t mpPOI2 = fReNITEBE[1][0][sd+2]->GetBinEntries(k);
2799     Double_t q1nRePOI2 = fReNITEBE[1][1][sd+2]->GetBinContent(k)*fReNITEBE[1][1][sd+2]->GetBinEntries(k);
2800     Double_t q1nImPOI2 = fImNITEBE[1][1][sd+2]->GetBinContent(k)*fImNITEBE[1][1][sd+2]->GetBinEntries(k);
2801     Double_t mqPOI2 = fReNITEBE[1][1][sd+2]->GetBinEntries(k);
2802     // Fill all-event profiles:
2803     if(weight>0. && mpPOI1>0.)
2804     {
2805      fNonIsotropicTermsVsEtaSumDiffPro[sd][0]->Fill(fEtaMin+(k-1)*fEtaBinWidth,p1nRePOI1/mpPOI1,mpPOI1); // <<cos(#psi_{POI_1})>>
2806      fNonIsotropicTermsVsEtaSumDiffPro[sd][1]->Fill(fEtaMin+(k-1)*fEtaBinWidth,p1nImPOI1/mpPOI1,mpPOI1); // <<sin(#psi_{POI_1})>>
2807     }
2808     if(weight>0. && mpPOI2>0.)
2809     {
2810      fNonIsotropicTermsVsEtaSumDiffPro[sd][2]->Fill(fEtaMin+(k-1)*fEtaBinWidth,p1nRePOI2/mpPOI2,mpPOI2); // <<cos(#psi_{POI_2})>>
2811      fNonIsotropicTermsVsEtaSumDiffPro[sd][3]->Fill(fEtaMin+(k-1)*fEtaBinWidth,p1nImPOI2/mpPOI2,mpPOI2); // <<sin(#psi_{POI_2})>>
2812     }
2813     if(weight>0. && mpPOI1*dMult-mqPOI1>0.)
2814     {
2815      fNonIsotropicTermsVsEtaSumDiffPro[sd][4]->Fill(fEtaMin+(k-1)*fEtaBinWidth,
2816       (p1nRePOI1*dReQ2n+p1nImPOI1*dImQ2n-q1nRePOI1)/(mpPOI1*dMult-mqPOI1),mpPOI1*dMult-mqPOI1); // <<cos(#psi_{POI_1}-2*phi)>>
2817      fNonIsotropicTermsVsEtaSumDiffPro[sd][5]->Fill(fEtaMin+(k-1)*fEtaBinWidth,
2818       (p1nImPOI1*dReQ2n-p1nRePOI1*dImQ2n+q1nImPOI1)/(mpPOI1*dMult-mqPOI1),mpPOI1*dMult-mqPOI1); // <<sin(#psi_{POI_1}-2*phi)>>
2819     }
2820     if(weight>0. && mpPOI2*dMult-mqPOI2>0.)
2821     {
2822      fNonIsotropicTermsVsEtaSumDiffPro[sd][6]->Fill(fEtaMin+(k-1)*fEtaBinWidth,
2823       (p1nRePOI2*dReQ2n+p1nImPOI2*dImQ2n-q1nRePOI2)/(mpPOI2*dMult-mqPOI2),mpPOI2*dMult-mqPOI2); // <<cos(#psi_{POI_2}-2*phi)>>
2824      fNonIsotropicTermsVsEtaSumDiffPro[sd][7]->Fill(fEtaMin+(k-1)*fEtaBinWidth,
2825       (p1nImPOI2*dReQ2n-p1nRePOI2*dImQ2n+q1nImPOI2)/(mpPOI2*dMult-mqPOI2),mpPOI2*dMult-mqPOI2); // <<sin(#psi_{POI_2}-2*phi)>>
2826     }
2827     if(weight>0. && mp>0.)
2828     {
2829      fNonIsotropicTermsVsEtaSumDiffPro[sd][8]->Fill(fEtaMin+(k-1)*fEtaBinWidth,p1nRe/mp,mp); // <<cos(#psi_{POI_1}+#psi_{POI_2})>>
2830      fNonIsotropicTermsVsEtaSumDiffPro[sd][9]->Fill(fEtaMin+(k-1)*fEtaBinWidth,p1nIm/mp,mp); // <<sin(#psi_{POI_1}+#psi_{POI_2})>>
2831     }
2832    } // end of for(Int_t k=1;k<=fnBinsEta;k++)
2833   } // end of for(Int_t sd=0;sd<2;sd++)
2834
2835   gIntegratedValue = -1000.;
2836   if((gSumWeight)&&(iBinCounter))
2837     gIntegratedValue = gSumBinContentTimesWeight/gSumWeight;
2838  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2839
2840  // b) Calculate differential 3-p correlator by using particle weights:
2841  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2842  {
2843   // ...
2844  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2845
2846 } // end of void AliFlowAnalysisWithMixedHarmonics::CalculateDifferential3pCorrelator()
2847
2848 //================================================================================================================
2849
2850 void AliFlowAnalysisWithMixedHarmonics::GetCorrelatorAndError(TProfile *g3pCorrelatorVsPt, Double_t &g3pCorrelatorValue, Double_t &g3pCorrelatorError) {
2851   //Retrieves the 3p correlator <<cos[n(psi1+psi2-2phi3)]>>
2852   //and its error
2853   Double_t gSumXi = 0.;
2854   Double_t gSumYi = 0.;
2855   Double_t gSumXiYi = 0.;
2856   Double_t gSumXiYi2 = 0.;
2857   Double_t gSumXi2Yi2 = 0.;
2858   Double_t gSumDeltaXi2 = 0.;
2859   Double_t gSumYi2DeltaXi2 = 0.;
2860
2861   for(Int_t iBin = 1; iBin <= g3pCorrelatorVsPt->GetNbinsX(); iBin++) {
2862     gSumXi += g3pCorrelatorVsPt->GetBinEntries(iBin);
2863     gSumYi += g3pCorrelatorVsPt->GetBinContent(iBin);
2864     gSumXiYi += g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin);
2865     gSumXiYi2 += g3pCorrelatorVsPt->GetBinEntries(iBin)*TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2);
2866     gSumXi2Yi2 += TMath::Power(g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin),2);
2867     gSumDeltaXi2 += TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);
2868     gSumYi2DeltaXi2 += TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2) + TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);
2869   }
2870
2871   g3pCorrelatorValue = -1000.;
2872   g3pCorrelatorError = 1000.;
2873
2874   if(gSumXi != 0.)
2875     g3pCorrelatorValue = gSumXiYi/gSumXi;
2876   if((gSumXi != 0.)&&(gSumXiYi != 0.))
2877     g3pCorrelatorError = TMath::Abs((gSumXiYi/gSumXi))*TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumYi2DeltaXi2)/gSumXiYi),2) + TMath::Power((gSumDeltaXi2/gSumXi),2));
2878 }
2879 //================================================================================================================
2880