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