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 fWeightsList(NULL),
  82 fUsePhiWeights(kFALSE),
  83 fUsePtWeights(kFALSE),
  84 fUseEtaWeights(kFALSE),
  85 fUseParticleWeights(NULL),
  86 fPhiWeights(NULL),
  87 fPtWeights(NULL),
  88 fEtaWeights(NULL),
  89 fReQnk(NULL),
  90 fImQnk(NULL),
  91 fSpk(NULL),
  92 fProfileList(NULL),
  93 f3pCorrelatorPro(NULL),
  94 f5pCorrelatorPro(NULL),
  95 fNonIsotropicTermsPro(NULL),
  96 f3pCorrelatorVsMPro(NULL),
  97 f3pPOICorrelatorVsM(NULL),
  98 fNonIsotropicTermsVsMPro(NULL),
  99 f2pCorrelatorCosPsiDiff(NULL),
 100 f2pCorrelatorCosPsiSum(NULL),
 101 f2pCorrelatorSinPsiDiff(NULL),
 102 f2pCorrelatorSinPsiSum(NULL),
 103 fResultsList(NULL),
 104 f3pCorrelatorHist(NULL),
 105 fDetectorBiasHist(NULL),
 106 f3pCorrelatorVsMHist(NULL),
 107 fDetectorBiasVsMHist(NULL)
 108 {
 109  // Constructor.
 110
 111  // Base list to hold all output objects:
 112  fHistList = new TList();
 113  fHistListName = new TString("cobjMH");
 114  fHistList->SetName(fHistListName->Data());
 115  fHistList->SetOwner(kTRUE);
 116
 117  // List to hold histograms with phi, pt and eta weights:
 118  fWeightsList = new TList();
 119
 120  // List to hold all all-event profiles:
 121  fProfileList = new TList();
 122
 123  // List to hold objects with final results:
 124  fResultsList = new TList();
 125
 126  // Initialize all arrays:
 127  this->InitializeArrays();
 128
 129 } // AliFlowAnalysisWithMixedHarmonics::AliFlowAnalysisWithMixedHarmonics()
 130
 131 //================================================================================================================
 132
 133 AliFlowAnalysisWithMixedHarmonics::~AliFlowAnalysisWithMixedHarmonics()
 134 {
 135  // Destructor.
 136
 137  delete fHistList;
 138
 139 } // end of AliFlowAnalysisWithMixedHarmonics::~AliFlowAnalysisWithMixedHarmonics()
 140
 141 //================================================================================================================
 142
 143 void AliFlowAnalysisWithMixedHarmonics::Init()
 144 {
 145  // Initialize and book all objects.
 146
 147  // a) Cross check if the user settings make sense before starting;
 148  // b) Access all common constants;
 149  // c) Book and nest all lists in the base list fHistList;
 150  // d) Book common control histograms;
 151  // e) Book all event-by-event quantities;
 152  // f) Book all all-event quantities;
 153  // g) Book and fill histograms to hold phi, pt and eta weights;
 154  // h) Store harmonic n used in cos[n*(phi1+phi2-2phi3)] and cos[n*(psi1+psi2-2phi3)].
 155
 156  //save old value and prevent histograms from being added to directory
 157  //to avoid name clashes in case multiple analaysis objects are used
 158  //in an analysis
 159  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 160  TH1::AddDirectory(kFALSE);
 161
 162  TH1::SetDefaultSumw2();
 163
 164  this->CrossCheckSettings();
 165  this->AccessConstants();
 166  this->BookAndNestAllLists();
 167  this->BookProfileHoldingSettings();
 168  this->BookCommonHistograms();
 169  this->BookAllEventByEventQuantities();
 170  this->BookAllAllEventQuantities();
 171  this->BookAndFillWeightsHistograms();
 172  this->StoreHarmonic();
 173
 174  TH1::AddDirectory(oldHistAddStatus);
 175
 176 } // end of void AliFlowAnalysisWithMixedHarmonics::Init()
 177
 178 //================================================================================================================
 179
 180 void AliFlowAnalysisWithMixedHarmonics::Make(AliFlowEventSimple* anEvent)
 181 {
 182  // Running over data only in this method.
 183
 184  // a) Check all pointers used in this method;
 185  // b) Define local variables;
 186  // c) Fill common control histograms;
 187  // d) Loop over data and calculate e-b-e quantities Q_{n,k} and S_{p,k};
 188  // 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};
 189  // f) Calculate differential 3-p azimuthal correlator cos[n(psi1+psi2-2*phi3)] in terms of Q_{2n} and p_{n}:
 190  // g) Reset all event-by-event quantities.
 191
 192  // a) Check all pointers used in this method:
 193  this->CheckPointersUsedInMake();
 194
 195  // b) Define local variables:
 196  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 197  Double_t dPt  = 0.; // transverse momentum
 198  Double_t dEta = 0.; // pseudorapidity
 199  Double_t wPhi = 1.; // phi weight
 200  Double_t wPt  = 1.; // pt weight
 201  Double_t wEta = 1.; // eta weight
 202  AliFlowTrackSimple *aftsTrack = NULL; // simple track
 203
 204  // c) Fill common control histograms:
 205  fCommonHists->FillControlHistograms(anEvent);
 206
 207  // d) Loop over data and calculate e-b-e quantities:
 208  Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI
 209                                            // nRP   = # of particles used to determine the reaction plane ("Reference Particles");
 210                                            // nPOI  = # of particles of interest for a detailed flow analysis ("Particles of Interest");
 211
 212  Int_t nRefMult = anEvent->GetReferenceMultiplicity();
 213
 214  // Start loop over data:
 215  for(Int_t i=0;i<nPrim;i++)
 216  {
 217   aftsTrack=anEvent->GetTrack(i);
 218   if(aftsTrack)
 219   {
 220    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue; // consider only tracks which are either RPs or POIs
 221    Int_t n = fHarmonic;
 222    if(aftsTrack->InRPSelection()) // checking RP condition:
 223    {
 224     dPhi = aftsTrack->Phi();
 225     dPt  = aftsTrack->Pt();
 226     dEta = aftsTrack->Eta();
 227     if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi-weight for this particle:
 228     {
 229      wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 230     }
 231     if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt-weight for this particle:
 232     {
 233      wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 234     }
 235     if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta-weight for this particle:
 236     {
 237      wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 238     }
 239     // 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:
 240     for(Int_t m=0;m<6;m++)
 241     {
 242      for(Int_t k=0;k<4;k++) // to be improved (what is the maximum k that I need?)
 243      {
 244       (*fReQnk)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1)*n*dPhi);
 245       (*fImQnk)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1)*n*dPhi);
 246      }
 247     }
 248     // Calculate partially S_{p,k} for this event (final calculation of S_{p,k} follows after the loop over data bellow):
 249     for(Int_t p=0;p<4;p++) // to be improved (what is maximum p that I need?)
 250     {
 251      for(Int_t k=0;k<4;k++) // to be improved (what is maximum k that I need?)
 252      {
 253       (*fSpk)(p,k)+=pow(wPhi*wPt*wEta,k);
 254      }
 255     }
 256    } // end of if(aftsTrack->InRPSelection())
 257    // POIs:
 258    if(fEvaluateDifferential3pCorrelator)
 259    {
 260     if(aftsTrack->InPOISelection()) // 1st POI
 261     {
 262      Double_t dPsi1 = aftsTrack->Phi();
 263      Double_t dPt1 = aftsTrack->Pt();
 264      Double_t dEta1 = aftsTrack->Eta();
 265      Int_t iCharge1 = aftsTrack->Charge();
 266      Bool_t b1stPOIisAlsoRP = kFALSE;
 267      if(aftsTrack->InRPSelection()){b1stPOIisAlsoRP = kTRUE;}
 268      for(Int_t j=0;j<nPrim;j++)
 269      {
 270       if(j==i){continue;}
 271       aftsTrack=anEvent->GetTrack(j);
 272       if(aftsTrack->InPOISelection()) // 2nd POI
 273       {
 274        Double_t dPsi2 = aftsTrack->Phi();
 275        Double_t dPt2 = aftsTrack->Pt();
 276        Double_t dEta2 = aftsTrack->Eta();
 277        Int_t iCharge2 = aftsTrack->Charge();
 278        if(fOppositeChargesPOI && iCharge1 == iCharge2){continue;}
 279        Bool_t b2ndPOIisAlsoRP = kFALSE;
 280        if(aftsTrack->InRPSelection()){b2ndPOIisAlsoRP = kTRUE;}
 281
 282        // Fill:Pt
 283        fRePEBE[0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1+dPsi2)),1.);
 284        fImPEBE[0]->Fill((dPt1+dPt2)/2.,TMath::Sin(n*(dPsi1+dPsi2)),1.);
 285        fRePEBE[1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1+dPsi2)),1.);
 286        fImPEBE[1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1+dPsi2)),1.);
 287
 288        // Fill:Eta
 289        fReEtaEBE[0]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1+dPsi2)),1.);
 290        fImEtaEBE[0]->Fill((dEta1+dEta2)/2.,TMath::Sin(n*(dPsi1+dPsi2)),1.);
 291        fReEtaEBE[1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1+dPsi2)),1.);
 292        fImEtaEBE[1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Sin(n*(dPsi1+dPsi2)),1.);
 293
 294        //2particle correlator <cos(n*(psi1 - ps12))>
 295        f2pCorrelatorCosPsiDiff->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1-dPsi2)));
 296        f2pCorrelatorCosPsiSum->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1+dPsi2)));
 297        f2pCorrelatorSinPsiDiff->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1-dPsi2)));
 298        f2pCorrelatorSinPsiSum->Fill(TMath::Abs(dPt1-dPt2),TMath::Sin(n*(dPsi1+dPsi2)));
 299
 300        if(b1stPOIisAlsoRP)
 301        {
 302         fOverlapEBE[0][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 303         fOverlapEBE[0][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 304         fOverlapEBE2[0][0]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 305         fOverlapEBE2[0][1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 306        }
 307        if(b2ndPOIisAlsoRP)
 308        {
 309         fOverlapEBE[1][0]->Fill((dPt1+dPt2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 310         fOverlapEBE[1][1]->Fill(TMath::Abs(dPt1-dPt2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 311         fOverlapEBE2[1][0]->Fill((dEta1+dEta2)/2.,TMath::Cos(n*(dPsi1-dPsi2)),1.);
 312         fOverlapEBE2[1][1]->Fill(TMath::Abs(dEta1-dEta2),TMath::Cos(n*(dPsi1-dPsi2)),1.);
 313        }
 314       } // end of if(aftsTrack->InPOISelection()) // 2nd POI
 315      } // end of for(Int_t j=i+1;j<nPrim;j++)
 316     } // end of if(aftsTrack->InPOISelection()) // 1st POI
 317    } // end of if(fEvaluateDifferential3pCorrelator)
 318   } else // to if(aftsTrack)
 319     {
 320      cout<<endl;
 321      cout<<" WARNING (MH): No particle! (i.e. aftsTrack is a NULL pointer in Make().)"<<endl;
 322      cout<<endl;
 323     }
 324  } // end of for(Int_t i=0;i<nPrim;i++)
 325
 326  // Calculate the final expressions for S_{p,k}:
 327  for(Int_t p=0;p<4;p++) // to be improved (what is maximum p that I need?)
 328  {
 329   for(Int_t k=0;k<4;k++) // to be improved (what is maximum k that I need?)
 330   {
 331    (*fSpk)(p,k)=pow((*fSpk)(p,k),p+1);
 332   }
 333  }
 334
 335  // e) Calculate 3-p correlator cos[n(phi1+phi2-2*phi3)] in terms of Q_{n,k} and S_{p,k}:
 336  if(anEvent->GetEventNSelTracksRP() >= 3)
 337  {
 338   this->Calculate3pCorrelator();
 339   this->CalculateNonIsotropicTerms();
 340   if(anEvent->GetEventNSelTracksRP() >= 5)
 341   {
 342    this->Calculate5pCorrelator();
 343   } // end of if(anEvent->GetEventNSelTracksRP() >= 5)
 344  } // end of if(anEvent->GetEventNSelTracksRP() >= 3)
 345
 346  // f) Calculate differential 3-p azimuthal correlator cos[n(psi1+psi2-2*phi3)] in terms of Q_{2n} and p_{n}:
 347  if(fEvaluateDifferential3pCorrelator && anEvent->GetEventNSelTracksRP() >= 1)
 348  {
 349    Double_t gIntegrated3pCorrelator = 0.;
 350    this->CalculateDifferential3pCorrelator(gIntegrated3pCorrelator); // to be improved - add relevant if statements for the min # POIs as well
 351
 352   //3particle correlator vs ref. mult
 353    if(fCalculateVsM)
 354      f3pPOICorrelatorVsM->Fill(nRefMult,gIntegrated3pCorrelator);
 355  }
 356
 357  // g) Reset all event-by-event quantities:
 358  this->ResetEventByEventQuantities();
 359
 360 } // end of AliFlowAnalysisWithMixedHarmonics::Make(AliFlowEventSimple* anEvent)
 361
 362 //================================================================================================================
 363
 364 void AliFlowAnalysisWithMixedHarmonics::Finish()
 365 {
 366  // Calculate the final results.
 367
 368  // a) Check all pointers used in this method;
 369  // b) Access settings for analysis with mixed harmonics;
 370  // c) Correct for detector effects;
 371  // d) Print on the screen the final results.
 372
 373  this->CheckPointersUsedInFinish();
 374  this->AccessSettings();
 375  this->CorrectForDetectorEffects();
 376  if(fCalculateVsM){this->CorrectForDetectorEffectsVsM();}
 377  if(fPrintOnTheScreen){this->PrintOnTheScreen();}
 378
 379 } // end of AliFlowAnalysisWithMixedHarmonics::Finish()
 380
 381 //================================================================================================================
 382
 383 void AliFlowAnalysisWithMixedHarmonics::GetOutputHistograms(TList *outputListHistos)
 384 {
 385  // Get pointers to all objects saved in the output file.
 386
 387  // a) Get pointers for common control histograms.
 388  if(outputListHistos)
 389  {
 390   this->SetHistList(outputListHistos);
 391   if(!fHistList)
 392   {
 393    cout<<endl;
 394    cout<<" WARNING (MH): fHistList is NULL in GetOutputHistograms() !!!!"<<endl;
 395    cout<<endl;
 396    exit(0);
 397   }
 398   this->GetPointersForBaseHistograms();
 399   this->GetPointersForCommonHistograms();
 400   this->GetPointersForAllEventProfiles();
 401   this->GetPointersForResultsHistograms();
 402  } else
 403    {
 404     cout<<endl;
 405     cout<<" WARNING (MH): outputListHistos is NULL in GetOutputHistograms() !!!!"<<endl;
 406     cout<<endl;
 407     exit(0);
 408    }
 409
 410 } // end of void AliFlowAnalysisWithMixedHarmonics::GetOutputHistograms(TList *outputListHistos)
 411
 412 //================================================================================================================
 413
 414 void AliFlowAnalysisWithMixedHarmonics::GetPointersForBaseHistograms()
 415 {
 416  // Get pointers to base histograms.
 417
 418  TString analysisSettingsName = "fAnalysisSettings";
 419  TProfile *analysisSettings = dynamic_cast<TProfile*>(fHistList->FindObject(analysisSettingsName.Data()));
 420  if(analysisSettings)
 421  {
 422   this->SetAnalysisSettings(analysisSettings);
 423  } else
 424    {
 425     cout<<endl;
 426     cout<<" WARNING (MH): analysisSettings is NULL in GetPointersForBaseHistograms() !!!!"<<endl;
 427     cout<<endl;
 428     exit(0);
 429    }
 430
 431 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForBaseHistograms()
 432
 433 //================================================================================================================
 434
 435 void AliFlowAnalysisWithMixedHarmonics::GetPointersForCommonHistograms()
 436 {
 437  // Get pointers to common control histograms.
 438
 439  TString commonHistsName = "AliFlowCommonHistMH";
 440  AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(fHistList->FindObject(commonHistsName.Data()));
 441  if(commonHist)
 442  {
 443   this->SetCommonHists(commonHist);
 444  } else
 445    {
 446     cout<<endl;
 447     cout<<" WARNING (MH): commonHist is NULL in GetPointersForCommonHistograms() !!!!"<<endl;
 448     cout<<endl;
 449     exit(0);
 450    }
 451
 452 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForCommonHistograms()
 453
 454 //================================================================================================================
 455
 456 void AliFlowAnalysisWithMixedHarmonics::GetPointersForAllEventProfiles()
 457 {
 458  // Get pointers to profiles holding final results.
 459
 460  TList *profileList = NULL;
 461  profileList = dynamic_cast<TList*>(fHistList->FindObject("Profiles"));
 462  if(!profileList)
 463  {
 464   cout<<endl;
 465   cout<<" WARNING (MH): profileList is NULL in GetPointersForAllEventProfiles() !!!!"<<endl;
 466   cout<<endl;
 467   exit(0);
 468  }
 469
 470  TString s3pCorrelatorProName = "f3pCorrelatorPro";
 471  TProfile *p3pCorrelatorPro = dynamic_cast<TProfile*>(profileList->FindObject(s3pCorrelatorProName.Data()));
 472  if(p3pCorrelatorPro)
 473  {
 474   this->Set3pCorrelatorPro(p3pCorrelatorPro);
 475  }
 476  TString s3pCorrelatorVsMProName = "f3pCorrelatorVsMPro";
 477  TProfile *p3pCorrelatorVsMPro = dynamic_cast<TProfile*>(profileList->FindObject(s3pCorrelatorVsMProName.Data()));
 478  if(p3pCorrelatorVsMPro)
 479  {
 480   this->Set3pCorrelatorVsMPro(p3pCorrelatorVsMPro);
 481  }
 482  TString s3pPOICorrelatorVsMName = "f3pPOICorrelatorVsM";
 483  TProfile *p3pPOICorrelatorVsM = dynamic_cast<TProfile*>(profileList->FindObject(s3pPOICorrelatorVsMName.Data()));
 484  if(p3pPOICorrelatorVsM)
 485  {
 486   this->Set3pPOICorrelatorVsM(p3pPOICorrelatorVsM);
 487  }
 488  TString nonIsotropicTermsProName = "fNonIsotropicTermsPro";
 489  TProfile *nonIsotropicTermsPro = dynamic_cast<TProfile*>(profileList->FindObject(nonIsotropicTermsProName.Data()));
 490  if(nonIsotropicTermsPro)
 491  {
 492   this->SetNonIsotropicTermsPro(nonIsotropicTermsPro);
 493  }
 494  TString nonIsotropicTermsVsMProName = "fNonIsotropicTermsVsMPro";
 495  TProfile2D *nonIsotropicTermsVsMPro = dynamic_cast<TProfile2D*>(profileList->FindObject(nonIsotropicTermsVsMProName.Data()));
 496  if(nonIsotropicTermsVsMPro)
 497  {
 498   this->SetNonIsotropicTermsVsMPro(nonIsotropicTermsVsMPro);
 499  }
 500  TString psdFlag[2] = {"PtSum","PtDiff"};
 501  for(Int_t sd=0;sd<2;sd++)
 502  {
 503   TProfile *p3pCorrelatorVsPtSumDiffPro = dynamic_cast<TProfile*>(profileList->FindObject(Form("f3pCorrelatorVs%sPro",psdFlag[sd].Data())));
 504   if(p3pCorrelatorVsPtSumDiffPro)
 505   {
 506    this->Set3pCorrelatorVsPtSumDiffPro(p3pCorrelatorVsPtSumDiffPro,sd);
 507   }
 508  }
 509
 510  TProfile *g2pCorrelatorCosPsiDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiDiff"));
 511  if(g2pCorrelatorCosPsiDiff)
 512    this->Set2pCorrelatorCosPsiDiff(g2pCorrelatorCosPsiDiff);
 513  TProfile *g2pCorrelatorCosPsiSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorCosPsiSum"));
 514  if(g2pCorrelatorCosPsiSum)
 515    this->Set2pCorrelatorCosPsiSum(g2pCorrelatorCosPsiSum);
 516  TProfile *g2pCorrelatorSinPsiDiff = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiDiff"));
 517  if(g2pCorrelatorSinPsiDiff)
 518    this->Set2pCorrelatorSinPsiDiff(g2pCorrelatorSinPsiDiff);
 519  TProfile *g2pCorrelatorSinPsiSum = dynamic_cast<TProfile *>(profileList->FindObject("f2pCorrelatorSinPsiSum"));
 520  if(g2pCorrelatorSinPsiSum)
 521    this->Set2pCorrelatorSinPsiSum(g2pCorrelatorSinPsiSum);
 522
 523  TString s5pCorrelatorProName = "f5pCorrelatorPro";
 524  TProfile *p5pCorrelatorPro = dynamic_cast<TProfile*>(profileList->FindObject(s5pCorrelatorProName.Data()));
 525  if(p5pCorrelatorPro)
 526  {
 527   this->Set5pCorrelatorPro(p5pCorrelatorPro);
 528  }
 529
 530 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForAllEventProfiles()
 531
 532 //================================================================================================================
 533
 534 void AliFlowAnalysisWithMixedHarmonics::GetPointersForResultsHistograms()
 535 {
 536  // Get pointers to histograms holding final results.
 537
 538  TList *resultsList = NULL;
 539  resultsList = dynamic_cast<TList*>(fHistList->FindObject("Results"));
 540  if(!resultsList)
 541  {
 542   cout<<endl;
 543   cout<<" WARNING (MH): resultsList is NULL in GetPointersForResultsHistograms() !!!!"<<endl;
 544   cout<<endl;
 545   exit(0);
 546  }
 547  TString s3pCorrelatorHistName = "f3pCorrelatorHist";
 548  TH1D *h3pCorrelatorHist = dynamic_cast<TH1D*>(resultsList->FindObject(s3pCorrelatorHistName.Data()));
 549  if(h3pCorrelatorHist)
 550  {
 551   this->Set3pCorrelatorHist(h3pCorrelatorHist);
 552  }
 553  TString s3pCorrelatorVsMHistName = "f3pCorrelatorVsMHist";
 554  TH1D *h3pCorrelatorVsMHist = dynamic_cast<TH1D*>(resultsList->FindObject(s3pCorrelatorVsMHistName.Data()));
 555  if(h3pCorrelatorVsMHist)
 556  {
 557   this->Set3pCorrelatorVsMHist(h3pCorrelatorVsMHist);
 558  }
 559  TString detectorBiasHistName = "fDetectorBiasHist";
 560  TH1D *detectorBiasHist = dynamic_cast<TH1D*>(resultsList->FindObject(detectorBiasHistName.Data()));
 561  if(detectorBiasHist)
 562  {
 563   this->SetDetectorBiasHist(detectorBiasHist);
 564  }
 565  TString detectorBiasVsMHistName = "fDetectorBiasVsMHist";
 566  TH1D *detectorBiasVsMHist = dynamic_cast<TH1D*>(resultsList->FindObject(detectorBiasVsMHistName.Data()));
 567  if(detectorBiasVsMHist)
 568  {
 569   this->SetDetectorBiasVsMHist(detectorBiasVsMHist);
 570  }
 571
 572 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForResultsHistograms()
 573
 574 //================================================================================================================
 575
 576 void AliFlowAnalysisWithMixedHarmonics::WriteHistograms(TString outputFileName)
 577 {
 578  // Store the final results in output .root file.
 579  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 580  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
 581  delete output;
 582 }
 583
 584 //================================================================================================================
 585
 586 void AliFlowAnalysisWithMixedHarmonics::WriteHistograms(TDirectoryFile *outputFileName)
 587 {
 588  // Store the final results in output .root file.
 589  fHistList->SetName("cobjMH");
 590  fHistList->SetOwner(kTRUE);
 591  outputFileName->Add(fHistList);
 592  outputFileName->Write(outputFileName->GetName(),TObject::kSingleKey);
 593 }
 594
 595 //================================================================================================================
 596
 597 void AliFlowAnalysisWithMixedHarmonics::StoreHarmonic()
 598 {
 599  // Store harmonic n used in cos[n*(phi1+phi2-2phi3)] and cos[n*(psi1+psi2-2phi3)].
 600
 601  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
 602
 603 } // end of void AliFlowAnalysisWithMixedHarmonics::StoreHarmonic()
 604
 605 //================================================================================================================
 606
 607 void AliFlowAnalysisWithMixedHarmonics::InitializeArrays()
 608 {
 609  // Initialize arrays.
 610
 611  for(Int_t sd=0;sd<2;sd++)
 612  {
 613   fRePEBE[sd] = NULL;
 614   fImPEBE[sd] = NULL;
 615   fReEtaEBE[sd] = NULL;
 616   fImEtaEBE[sd] = NULL;
 617   f3pCorrelatorVsPtSumDiffPro[sd] = NULL;
 618   f3pCorrelatorVsEtaSumDiffPro[sd] = NULL;
 619  }
 620  for(Int_t fs=0;fs<2;fs++) // 1st/2nd POI which is also RP
 621  {
 622   for(Int_t sd=0;sd<2;sd++)
 623   {
 624    fOverlapEBE[fs][sd] = NULL;
 625    fOverlapEBE2[fs][sd] = NULL;
 626   }
 627  }
 628
 629 } // end of AliFlowAnalysisWithMixedHarmonics::InitializeArrays()
 630
 631 //================================================================================================================
 632
 633 void AliFlowAnalysisWithMixedHarmonics::BookAndNestAllLists()
 634 {
 635  // Book and nest all list in base list fHistList.
 636
 637  // Weights:
 638  fWeightsList->SetName("Weights");
 639  fWeightsList->SetOwner(kTRUE);
 640  fHistList->Add(fWeightsList);
 641  // Profiles:
 642  fProfileList->SetName("Profiles");
 643  fProfileList->SetOwner(kTRUE);
 644  fHistList->Add(fProfileList);
 645  // Results:
 646  fResultsList->SetName("Results");
 647  fResultsList->SetOwner(kTRUE);
 648  fHistList->Add(fResultsList);
 649
 650 } // end of void AliFlowAnalysisWithMixedHarmonics::BookAndNestAllLists()
 651
 652 //================================================================================================================
 653
 654 void AliFlowAnalysisWithMixedHarmonics::BookProfileHoldingSettings()
 655 {
 656  // Book profile to hold all analysis settings.
 657
 658  TString analysisSettingsName = "fAnalysisSettings";
 659  fAnalysisSettings = new TProfile(analysisSettingsName.Data(),"Settings for analysis with mixed harmonics",10,0,10);
 660  fAnalysisSettings->SetStats(kFALSE);
 661  fAnalysisSettings->GetXaxis()->SetLabelSize(0.03);
 662  fAnalysisSettings->GetXaxis()->SetBinLabel(1,"Corr. for det. effects?");
 663  fAnalysisSettings->Fill(0.5,(Int_t)fCorrectForDetectorEffects);
 664  fAnalysisSettings->GetXaxis()->SetBinLabel(2,"# of mult. bins");
 665  fAnalysisSettings->Fill(1.5,fNoOfMultipicityBins);
 666  fAnalysisSettings->GetXaxis()->SetBinLabel(3,"Width of mult. bins");
 667  fAnalysisSettings->Fill(2.5,fMultipicityBinWidth);
 668  fAnalysisSettings->GetXaxis()->SetBinLabel(4,"Minimal mult.");
 669  fAnalysisSettings->Fill(3.5,fMinMultiplicity);
 670  fAnalysisSettings->GetXaxis()->SetBinLabel(5,"Print on the screen?");
 671  fAnalysisSettings->Fill(4.5,(Int_t)fPrintOnTheScreen);
 672  fAnalysisSettings->GetXaxis()->SetBinLabel(6,"fHarmonic");
 673  fAnalysisSettings->Fill(5.5,(Int_t)fHarmonic);
 674  fAnalysisSettings->GetXaxis()->SetBinLabel(7,"fOppositeChargesPOI");
 675  fAnalysisSettings->Fill(6.5,(Int_t)fOppositeChargesPOI);
 676  fAnalysisSettings->GetXaxis()->SetBinLabel(8,"fEvaluateDifferential3pCorrelator");
 677  fAnalysisSettings->Fill(7.5,(Int_t)fOppositeChargesPOI);
 678  fAnalysisSettings->GetXaxis()->SetBinLabel(9,"fCalculateVsM");
 679  fAnalysisSettings->Fill(8.5,(Int_t)fCalculateVsM);
 680  fAnalysisSettings->GetXaxis()->SetBinLabel(10,"fShowBinLabelsVsM");
 681  fAnalysisSettings->Fill(9.5,(Int_t)fShowBinLabelsVsM);
 682  fHistList->Add(fAnalysisSettings);
 683
 684 } // end of void AliFlowAnalysisWithMixedHarmonics::BookProfileHoldingSettings()
 685
 686 //================================================================================================================
 687
 688 void AliFlowAnalysisWithMixedHarmonics::BookCommonHistograms()
 689 {
 690  // Book common control histograms and common histograms for final results.
 691
 692  TString commonHistsName = "AliFlowCommonHistMH";
 693  fCommonHists = new AliFlowCommonHist(commonHistsName.Data());
 694  fHistList->Add(fCommonHists);
 695
 696 } // end of void AliFlowAnalysisWithMixedHarmonics::BookCommonHistograms()
 697
 698 //================================================================================================================
 699
 700 void AliFlowAnalysisWithMixedHarmonics::BookAllEventByEventQuantities()
 701 {
 702  // Book all event-by-event quantitites.
 703
 704   // Q_{n,k} and S{p,k}:
 705  fReQnk = new TMatrixD(6,9); // to be improved (check bound on k!)
 706  fImQnk = new TMatrixD(6,9); // to be improved (check bound on k!)
 707  fSpk = new TMatrixD(4,4); // to be improved (check bound on p and k!)
 708
 709  // p_n vs [(p1+p2)/2,|p1-p2|]
 710  if(!fEvaluateDifferential3pCorrelator){return;}
 711  TString psdFlag[2] = {"PtSum","PtDiff"};
 712  TString p2sdFlag[2] = {"PtSum","PtDiff"};
 713  TString fsFlag[2] = {"1st","2nd"};
 714  for(Int_t sd=0;sd<2;sd++)
 715  {
 716   fRePEBE[sd] = new TProfile(Form("fRePEBE%s",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
 717   fImPEBE[sd] = new TProfile(Form("fImPEBE%s",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
 718   fReEtaEBE[sd] = new TProfile(Form("fReEtaEBE%s",p2sdFlag[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
 719   fImEtaEBE[sd] = new TProfile(Form("fImEtaEBE%s",p2sdFlag[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
 720  }
 721  for(Int_t fs=0;fs<2;fs++)
 722  {
 723   for(Int_t sd=0;sd<2;sd++)
 724   {
 725    fOverlapEBE[fs][sd] = new TProfile(Form("%s POI, %s",fsFlag[sd].Data(),psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
 726    fOverlapEBE2[fs][sd] = new TProfile(Form("%s POI 2, %s",fsFlag[sd].Data(),p2sdFlag[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
 727   }
 728  }
 729
 730 } // end fo void AliFlowAnalysisWithMixedHarmonics::BookAllEventByEventQuantities()
 731
 732 //================================================================================================================
 733
 734 void AliFlowAnalysisWithMixedHarmonics::BookAllAllEventQuantities()
 735 {
 736  // Book all all-event quantitites.
 737
 738  // a) Book histos and profiles without any binning in multiplicity, pt or eta;
 739  // b) Book quantites with multiplicity binning;
 740  // c) Book quantites with binning in (p1+p2)/2 and |p1-p2|.
 741
 742  this->BookDefault();
 743  if(fCalculateVsM){this->BookVsM();}
 744  if(fEvaluateDifferential3pCorrelator){this->BookDifferential();}
 745
 746 } // end of void AliFlowAnalysisWithMixedHarmonics::BookAllAllEventQuantities()
 747
 748 //================================================================================================================
 749
 750 void AliFlowAnalysisWithMixedHarmonics::BookDefault()
 751 {
 752  // Book histos and profiles without any binning in multiplicity, pt or eta.
 753
 754  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects);
 755  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>>;
 756  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects;
 757  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>>;
 758  // e) 5-p correlator <<cos[n*(2phi1+2phi2+2phi3-3phi4-3phi5)]>> for all events (not corrected for detector effects - not supported yet).
 759
 760  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects);
 761  TString s3pCorrelatorProName = "f3pCorrelatorPro";
 762  f3pCorrelatorPro = new TProfile(s3pCorrelatorProName.Data(),"",1,0,1);
 763  f3pCorrelatorPro->SetStats(kFALSE);
 764  f3pCorrelatorPro->GetXaxis()->SetLabelOffset(0.01);
 765  f3pCorrelatorPro->GetXaxis()->SetLabelSize(0.05);
 766  if(fHarmonic == 1)
 767  {
 768   f3pCorrelatorPro->GetXaxis()->SetBinLabel(1,"#LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
 769  } else
 770    {
 771     f3pCorrelatorPro->GetXaxis()->SetBinLabel(1,Form("#LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
 772    }
 773  fProfileList->Add(f3pCorrelatorPro);
 774
 775  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>>:
 776  TString nonIsotropicTermsProName = "fNonIsotropicTermsPro";
 777  fNonIsotropicTermsPro = new TProfile(nonIsotropicTermsProName.Data(),"",8,0,8);
 778  fNonIsotropicTermsPro->SetStats(kFALSE);
 779  if(fHarmonic == 1)
 780  {
 781   fNonIsotropicTermsPro->SetTitle("Non-isotropic terms in decomposition of #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
 782   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(1,"cos(#phi_{1})");
 783   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(2,"sin(#phi_{1})");
 784   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(3,"cos(2#phi_{1})");
 785   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(4,"sin(2#phi_{1})");
 786   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(5,"cos(#phi_{1}+#phi_{2})");
 787   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(6,"sin(#phi_{1}+#phi_{2})");
 788   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(7,"cos(2#phi_{1}-#phi_{2})");
 789   fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(8,"sin(2#phi_{1}-#phi_{2})");
 790   // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(9,"cos(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
 791   // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(10,"sin(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
 792  } else
 793    {
 794     fNonIsotropicTermsPro->SetTitle(Form("Non-isotropic terms in decomposition of #LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
 795     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(1,Form("cos(%d#phi_{1})",fHarmonic));
 796     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(2,Form("sin(%d#phi_{1})",fHarmonic));
 797     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(3,Form("cos(%d#phi_{1})",2*fHarmonic));
 798     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(4,Form("sin(%d#phi_{1})",2*fHarmonic));
 799     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(5,Form("cos[%d(#phi_{1}+#phi_{2})]",fHarmonic));
 800     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(6,Form("sin[%d(#phi_{1}+#phi_{2})]",fHarmonic));
 801     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(7,Form("cos[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
 802     fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(8,Form("sin[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
 803     // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(9,Form("cos(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
 804     // fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(10,Form("sin(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
 805    }
 806  fProfileList->Add(fNonIsotropicTermsPro);
 807
 808  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects:
 809  TString s3pCorrelatorHistName = "f3pCorrelatorHist";
 810  f3pCorrelatorHist = new TH1D(s3pCorrelatorHistName.Data(),"",1,0,1);
 811  f3pCorrelatorHist->SetStats(kFALSE);
 812  f3pCorrelatorHist->GetXaxis()->SetLabelOffset(0.01);
 813  f3pCorrelatorHist->GetXaxis()->SetLabelSize(0.05);
 814  if(fHarmonic == 1)
 815  {
 816   f3pCorrelatorHist->GetXaxis()->SetBinLabel(1,"#LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
 817  } else
 818    {
 819     f3pCorrelatorHist->GetXaxis()->SetBinLabel(1,Form("#LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
 820    }
 821  fResultsList->Add(f3pCorrelatorHist);
 822
 823  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>>:
 824  TString detectorBiasHistName = "fDetectorBiasHist";
 825  fDetectorBiasHist = new TH1D(detectorBiasHistName.Data(),"Bias coming from detector inefficiences",1,0,1);
 826  fDetectorBiasHist->SetStats(kFALSE);
 827  if(fHarmonic == 1)
 828  {
 829   fDetectorBiasHist->GetXaxis()->SetBinLabel(1,"#frac{corrected}{measured} #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT");
 830  } else
 831    {
 832     fDetectorBiasHist->GetXaxis()->SetBinLabel(1,Form("#frac{corrected}{measured} #LT#LTcos[%i(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
 833    }
 834  fResultsList->Add(fDetectorBiasHist);
 835
 836  // e) 5-p correlator <<cos[n*(2phi1+2phi2+2phi3-3phi4-3phi5)]>> for all events (not corrected for detector effects - not supported yet):
 837  TString s5pCorrelatorProName = "f5pCorrelatorPro";
 838  f5pCorrelatorPro = new TProfile(s5pCorrelatorProName.Data(),"",1,0,1);
 839  f5pCorrelatorPro->SetStats(kFALSE);
 840  f5pCorrelatorPro->GetXaxis()->SetLabelOffset(0.01);
 841  f5pCorrelatorPro->GetXaxis()->SetLabelSize(0.05);
 842  if(fHarmonic == 1)
 843  {
 844   f5pCorrelatorPro->GetXaxis()->SetBinLabel(1,"#LT#LTcos(2#phi_{1}+2#phi_{2}+2#phi_{3}-3#phi_{4}-3#phi_{5})#GT#GT");
 845  } else
 846    {
 847     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));
 848    }
 849  fProfileList->Add(f5pCorrelatorPro);
 850
 851 } // end of void AliFlowAnalysisWithMixedHarmonics::BookDefault()
 852
 853 //================================================================================================================
 854
 855 void AliFlowAnalysisWithMixedHarmonics::BookVsM()
 856 {
 857  // Book histos and profiles holding results vs multiplicity.
 858
 859  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects) vs M;
 860  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>> vs M;
 861  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects vs M;
 862  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> vs M.
 863
 864  // a) 3-p correlator <<cos[n*(phi1+phi2-2phi3)]>> for all events (not corrected for detector effects) vs M:
 865  TString s3pCorrelatorVsMProName = "f3pCorrelatorVsMPro";
 866  f3pCorrelatorVsMPro = new TProfile(s3pCorrelatorVsMProName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
 867  f3pCorrelatorVsMPro->SetStats(kFALSE);
 868  if(fHarmonic == 1)
 869  {
 870   f3pCorrelatorVsMPro->SetTitle("#LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
 871  } else
 872    {
 873     f3pCorrelatorVsMPro->SetTitle(Form("#LT#LTcos[%d(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} M",fHarmonic));
 874    }
 875  if(fShowBinLabelsVsM)
 876  {
 877   f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
 878   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
 879   {
 880    f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
 881   }
 882   f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 883  } else
 884    {
 885     f3pCorrelatorVsMPro->GetXaxis()->SetTitle("M");
 886    }
 887  fProfileList->Add(f3pCorrelatorVsMPro);
 888
 889  TString s3pPOICorrelatorVsMName = "f3pPOICorrelatorVsM";
 890  f3pPOICorrelatorVsM = new TProfile(s3pPOICorrelatorVsMName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
 891  f3pPOICorrelatorVsM->SetStats(kFALSE);
 892  if(fHarmonic == 1)
 893  {
 894   f3pPOICorrelatorVsM->SetTitle("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
 895  } else
 896    {
 897     f3pPOICorrelatorVsM->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} M",fHarmonic));
 898    }
 899  if(fShowBinLabelsVsM)
 900  {
 901   f3pPOICorrelatorVsM->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
 902   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
 903   {
 904    f3pPOICorrelatorVsM->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
 905   }
 906   f3pPOICorrelatorVsM->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 907  } else
 908    {
 909     f3pPOICorrelatorVsM->GetXaxis()->SetTitle("M");
 910    }
 911  fProfileList->Add(f3pPOICorrelatorVsM);
 912
 913  // b) Non-isotropic terms in the decomposition of <<cos[n(phi1+phi2-2phi3)]>> vs M:
 914  TString s3pCorrelatorVsMHistName = "f3pCorrelatorVsMHist";
 915  f3pCorrelatorVsMHist = new TH1D(s3pCorrelatorVsMHistName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
 916  f3pCorrelatorVsMHist->SetStats(kFALSE);
 917  if(fHarmonic == 1)
 918  {
 919   f3pCorrelatorVsMHist->SetTitle("cos(#phi_{1}+#phi_{2}-2#phi_{3}) #font[72]{vs} M");
 920  } else
 921    {
 922     f3pCorrelatorVsMHist->SetTitle(Form("cos[%d(#phi_{1}+#phi_{2}-2#phi_{3})] #font[72]{vs} M",fHarmonic));
 923    }
 924  if(fShowBinLabelsVsM)
 925  {
 926   f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
 927   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
 928   {
 929    f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
 930   }
 931   f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 932  } else
 933    {
 934     f3pCorrelatorVsMHist->GetXaxis()->SetTitle("M");
 935    }
 936  fResultsList->Add(f3pCorrelatorVsMHist);
 937
 938  // c) 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> corrected for detector effects vs M:
 939  TString nonIsotropicTermsVsMProName = "fNonIsotropicTermsVsMPro";
 940  fNonIsotropicTermsVsMPro = new TProfile2D(nonIsotropicTermsVsMProName.Data(),"",8,0,8,fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
 941  fNonIsotropicTermsVsMPro->SetStats(kFALSE);
 942  if(fHarmonic == 1)
 943  {
 944   fNonIsotropicTermsVsMPro->SetTitle("Non-isotropic terms in decomposition of #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
 945   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(1,"cos(#phi_{1})");
 946   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(2,"sin(#phi_{1})");
 947   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(3,"cos(2#phi_{1})");
 948   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(4,"sin(2#phi_{1})");
 949   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(5,"cos(#phi_{1}+#phi_{2})");
 950   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(6,"sin(#phi_{1}+#phi_{2})");
 951   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(7,"cos(2#phi_{1}-#phi_{2})");
 952   fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(8,"sin(2#phi_{1}-#phi_{2})");
 953   // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(9,"cos(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
 954   // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(10,"sin(#phi_{1}-#phi_{2}-#phi_{3})"); // not needed
 955  } else
 956    {
 957     fNonIsotropicTermsVsMPro->SetTitle(Form("Non-isotropic terms in decomposition of #LT#LTcos[%d(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT",fHarmonic));
 958     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(1,Form("cos(%d#phi_{1})",fHarmonic));
 959     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(2,Form("sin(%d#phi_{1})",fHarmonic));
 960     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(3,Form("cos(%d#phi_{1})",2*fHarmonic));
 961     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(4,Form("sin(%d#phi_{1})",2*fHarmonic));
 962     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(5,Form("cos[%d(#phi_{1}+#phi_{2})]",fHarmonic));
 963     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(6,Form("sin[%d(#phi_{1}+#phi_{2})]",fHarmonic));
 964     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(7,Form("cos[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
 965     fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(8,Form("sin[%d(2#phi_{1}-#phi_{2})]",fHarmonic));
 966     // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(9,Form("cos(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
 967     // fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(10,Form("sin(%d(#phi_{1}-#phi_{2}-#phi_{3}))",fHarmonic)); // not needed
 968    }
 969  if(fShowBinLabelsVsM)
 970  {
 971   fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
 972   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
 973   {
 974    fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
 975   }
 976   fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 977  } else
 978    {
 979     fNonIsotropicTermsVsMPro->GetYaxis()->SetTitle("M");
 980    }
 981  fProfileList->Add(fNonIsotropicTermsVsMPro);
 982
 983  // d) Histogram which quantifies bias coming from detector inefficiencies to 3-p correlator <<cos[n(phi1+phi2-2phi3)]>> vs M:
 984  TString detectorBiasVsMHistName = "fDetectorBiasVsMHist";
 985  fDetectorBiasVsMHist = new TH1D(detectorBiasVsMHistName.Data(),"",fNoOfMultipicityBins+2,0,fNoOfMultipicityBins+2);
 986  fDetectorBiasVsMHist->SetStats(kFALSE);
 987  if(fHarmonic == 1)
 988  {
 989   fDetectorBiasVsMHist->SetTitle("#frac{corrected}{measured} #LT#LTcos(#phi_{1}+#phi_{2}-2#phi_{3})#GT#GT #font[72]{vs} M");
 990  } else
 991    {
 992     fDetectorBiasVsMHist->SetTitle(Form("#frac{corrected}{measured} #LT#LTcos[%d(#phi_{1}+#phi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} M",fHarmonic));
 993    }
 994  if(fShowBinLabelsVsM)
 995  {
 996   fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(1,Form("M < %d",(Int_t)fMinMultiplicity));
 997   for(Int_t b=2;b<=fNoOfMultipicityBins+1;b++)
 998   {
 999    fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-2)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth)));
1000   }
1001   fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+2,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
1002  } else
1003    {
1004     fDetectorBiasVsMHist->GetXaxis()->SetTitle("M");
1005    }
1006  fResultsList->Add(fDetectorBiasVsMHist);
1007
1008 } // end of void AliFlowAnalysisWithMixedHarmonics::BookVsM()
1009
1010 //================================================================================================================
1011
1012 void AliFlowAnalysisWithMixedHarmonics::BookDifferential()
1013 {
1014  // Book histos and profiles holding results vs (p1+p2)/2 and |p1-p2|.
1015
1016  TString psdFlag[2] = {"PtSum","PtDiff"};
1017  TString psdTitleFlag[2] = {"(p_{T,1}+ p_{T,2})/2","#left|p_{T,1}- p_{T,2}#right|"};
1018  TString psdFlag2[2] = {"EtaSum","EtaDiff"};
1019  TString psdTitleFlag2[2] = {"(#eta_{1}+ #eta_{2})/2","#left|#eta_{1}- #eta_{2}#right|"};
1020  //TString s3pCorrelatorVsPtSumDiffProName = "f3pCorrelatorVsPtSumDiffPro";
1021  for(Int_t sd=0;sd<2;sd++)
1022  {
1023   f3pCorrelatorVsPtSumDiffPro[sd] = new TProfile(Form("f3pCorrelatorVs%sPro",psdFlag[sd].Data()),"",fnBinsPt,0.,fPtMax);
1024   f3pCorrelatorVsPtSumDiffPro[sd]->SetStats(kFALSE);
1025   f3pCorrelatorVsEtaSumDiffPro[sd] = new TProfile(Form("f3pCorrelatorVs%sPro",psdFlag2[sd].Data()),"",fnBinsEta,fEtaMin,fEtaMax);
1026   f3pCorrelatorVsEtaSumDiffPro[sd]->SetStats(kFALSE);
1027   //f3pCorrelatorVsPtSumDiffPro[sd]->SetLabelSize(0.05);
1028   //f3pCorrelatorVsPtSumDiffPro[sd]->SetMarkerStyle(25);
1029   if(fHarmonic == 1)
1030   {
1031    f3pCorrelatorVsPtSumDiffPro[sd]->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} %s",psdTitleFlag[sd].Data()));
1032    f3pCorrelatorVsEtaSumDiffPro[sd]->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3})#GT#GT #font[72]{vs} %s",psdTitleFlag2[sd].Data()));
1033   } else
1034     {
1035      f3pCorrelatorVsPtSumDiffPro[sd]->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag[sd].Data()));
1036      f3pCorrelatorVsEtaSumDiffPro[sd]->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2}-2#phi_{3})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[sd].Data()));
1037     }
1038   f3pCorrelatorVsPtSumDiffPro[sd]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
1039   fProfileList->Add(f3pCorrelatorVsPtSumDiffPro[sd]);
1040   f3pCorrelatorVsEtaSumDiffPro[sd]->GetXaxis()->SetTitle(psdTitleFlag2[sd].Data());
1041   fProfileList->Add(f3pCorrelatorVsEtaSumDiffPro[sd]);
1042  }
1043
1044  f2pCorrelatorCosPsiDiff = new TProfile("f2pCorrelatorCosPsiDiff","",fnBinsPt,0.,fPtMax);
1045  f2pCorrelatorCosPsiDiff->SetStats(kFALSE);
1046  f2pCorrelatorCosPsiSum = new TProfile("f2pCorrelatorCosPsiSum","",fnBinsPt,0.,fPtMax);
1047  f2pCorrelatorCosPsiSum->SetStats(kFALSE);
1048  f2pCorrelatorSinPsiDiff = new TProfile("f2pCorrelatorSinPsiDiff","",fnBinsPt,0.,fPtMax);
1049  f2pCorrelatorSinPsiDiff->SetStats(kFALSE);
1050  f2pCorrelatorSinPsiSum = new TProfile("f2pCorrelatorSinPsiSum","",fnBinsPt,0.,fPtMax);
1051  f2pCorrelatorSinPsiSum->SetStats(kFALSE);
1052  if(fHarmonic == 1) {
1053    f2pCorrelatorCosPsiDiff->SetTitle(Form("#LT#LTcos(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1054    f2pCorrelatorCosPsiSum->SetTitle(Form("#LT#LTcos(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1055    f2pCorrelatorSinPsiDiff->SetTitle(Form("#LT#LTsin(#psi_{1}-#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1056    f2pCorrelatorSinPsiSum->SetTitle(Form("#LT#LTsin(#psi_{1}+#psi_{2})#GT#GT #font[72]{vs} %s",psdTitleFlag[1].Data()));
1057  }
1058  else {
1059    f2pCorrelatorCosPsiDiff->SetTitle(Form("#LT#LTcos[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1060    f2pCorrelatorCosPsiSum->SetTitle(Form("#LT#LTcos[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1061    f2pCorrelatorSinPsiDiff->SetTitle(Form("#LT#LTsin[%d(#psi_{1}-#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1062    f2pCorrelatorSinPsiSum->SetTitle(Form("#LT#LTsin[%d(#psi_{1}+#psi_{2})]#GT#GT #font[72]{vs} %s",fHarmonic,psdTitleFlag2[1].Data()));
1063  }
1064  fProfileList->Add(f2pCorrelatorCosPsiDiff);
1065  fProfileList->Add(f2pCorrelatorCosPsiSum);
1066  fProfileList->Add(f2pCorrelatorSinPsiDiff);
1067  fProfileList->Add(f2pCorrelatorSinPsiSum);
1068
1069 } // end of void AliFlowAnalysisWithMixedHarmonics::BookDifferential()
1070
1071 //================================================================================================================
1072
1073 void AliFlowAnalysisWithMixedHarmonics::AccessConstants()
1074 {
1075  // Access needed common constants from AliFlowCommonConstants.
1076
1077  fnBinsPhi = AliFlowCommonConstants::GetMaster()->GetNbinsPhi();
1078  fPhiMin = AliFlowCommonConstants::GetMaster()->GetPhiMin();
1079  fPhiMax = AliFlowCommonConstants::GetMaster()->GetPhiMax();
1080  if(fnBinsPhi) fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;
1081  fnBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
1082  fPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
1083  fPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
1084  if(fnBinsPt) fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;
1085  fnBinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
1086  fEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
1087  fEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
1088  if(fnBinsEta) fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;
1089
1090 } // end of void AliFlowAnalysisWithMixedHarmonics::AccessConstants()
1091
1092 //================================================================================================================
1093
1094 void AliFlowAnalysisWithMixedHarmonics::CrossCheckSettings()
1095 {
1096  // Cross-check if the user settings make sense.
1097
1098  // ...
1099
1100 } // end of void AliFlowAnalysisWithMixedHarmonics::CrossCheckSettings()
1101
1102 //================================================================================================================
1103
1104 void AliFlowAnalysisWithMixedHarmonics::BookAndFillWeightsHistograms()
1105 {
1106  // Book and fill (by accessing file "weights.root") histograms which hold phi, pt and eta weights.
1107
1108  if(!fWeightsList)
1109  {
1110   cout<<endl;
1111   cout<<" WARNING (MH): fWeightsList is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1112   cout<<endl;
1113   exit(0);
1114  }
1115  // Profile to hold flags for weights:
1116  TString fUseParticleWeightsName = "fUseParticleWeightsMH";
1117  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
1118  fUseParticleWeights->SetStats(kFALSE);
1119  fUseParticleWeights->SetLabelSize(0.06);
1120  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
1121  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
1122  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
1123  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
1124  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
1125  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
1126  fWeightsList->Add(fUseParticleWeights);
1127  // Phi-weights:
1128  if(fUsePhiWeights)
1129  {
1130   if(fWeightsList->FindObject("phi_weights"))
1131   {
1132    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
1133    if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
1134    {
1135     cout<<endl;
1136     cout<<" WARNING (MH): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
1137     cout<<endl;
1138     exit(0);
1139    }
1140   } else
1141     {
1142      cout<<endl;
1143      cout<<" WARNING (MH): fWeightsList->FindObject(\"phi_weights\") is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1144      cout<<endl;
1145      exit(0);
1146     }
1147  } // end of if(fUsePhiWeights)
1148  // Pt-weights:
1149  if(fUsePtWeights)
1150  {
1151   if(fWeightsList->FindObject("pt_weights"))
1152   {
1153    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
1154    if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
1155    {
1156     cout<<endl;
1157     cout<<" WARNING (MH): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
1158     cout<<endl;
1159     exit(0);
1160    }
1161   } else
1162     {
1163      cout<<endl;
1164      cout<<" WARNING (MH): fWeightsList->FindObject(\"pt_weights\") is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1165      cout<<endl;
1166      exit(0);
1167     }
1168  } // end of if(fUsePtWeights)
1169  // Eta-weights:
1170  if(fUseEtaWeights)
1171  {
1172   if(fWeightsList->FindObject("eta_weights"))
1173   {
1174    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
1175    if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
1176    {
1177     cout<<endl;
1178     cout<<" WARNING (MH): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
1179     cout<<endl;
1180     exit(0);
1181    }
1182   } else
1183     {
1184      cout<<endl;
1185      cout<<" WARNING (MH): fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in BookAndFillWeightsHistograms() !!!!"<<endl;
1186      cout<<endl;
1187      exit(0);
1188     }
1189  } // end of if(fUseEtaWeights)
1190
1191 } // end of AliFlowAnalysisWithMixedHarmonics::BookAndFillWeightsHistograms()
1192
1193 //================================================================================================================
1194
1195 void AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInMake()
1196 {
1197  // Check pointers used in method Make().
1198
1199  if(!fReQnk || !fImQnk || !fSpk )
1200  {
1201   cout<<endl;
1202   cout<<" WARNING (MH): fReQnk || fImQnk || fSpk is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1203   cout<<endl;
1204   exit(0);
1205  }
1206  if(!f3pCorrelatorPro)
1207  {
1208   cout<<endl;
1209   cout<<" WARNING (MH): f3pCorrelatorPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1210   cout<<endl;
1211   exit(0);
1212  }
1213  if(!f5pCorrelatorPro)
1214  {
1215   cout<<endl;
1216   cout<<" WARNING (MH): f5pCorrelatorPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1217   cout<<endl;
1218   exit(0);
1219  }
1220  if(!fNonIsotropicTermsPro)
1221  {
1222   cout<<endl;
1223   cout<<" WARNING (MH): fNonIsotropicTermsPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1224   cout<<endl;
1225   exit(0);
1226  }
1227  if(!f3pCorrelatorVsMPro && fCalculateVsM)
1228  {
1229   cout<<endl;
1230   cout<<" WARNING (MH): f3pCorrelatorVsMPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1231   cout<<endl;
1232   exit(0);
1233  }
1234  if(!f3pPOICorrelatorVsM && fCalculateVsM)
1235  {
1236   cout<<endl;
1237   cout<<" WARNING (MH): f3pPOICorrelatorVsM is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1238   cout<<endl;
1239   exit(0);
1240  }
1241  if(!fNonIsotropicTermsVsMPro && fCalculateVsM)
1242  {
1243   cout<<endl;
1244   cout<<" WARNING (MH): fNonIsotropicTermsVsMPro is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1245   cout<<endl;
1246   exit(0);
1247  }
1248
1249  // Differential correlators:
1250  if(!fEvaluateDifferential3pCorrelator){return;}
1251  for(Int_t sd=0;sd<2;sd++)
1252  {
1253   if(!(f3pCorrelatorVsPtSumDiffPro[sd]))
1254   {
1255    cout<<endl;
1256    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsPtSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1257    cout<<endl;
1258    exit(0);
1259   }
1260   if(!(f3pCorrelatorVsEtaSumDiffPro[sd]))
1261   {
1262    cout<<endl;
1263    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsEtaSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1264    cout<<endl;
1265    exit(0);
1266   }
1267  }
1268  for(Int_t sd=0;sd<2;sd++)
1269  {
1270   if(!fRePEBE[sd]||!fImPEBE[sd])
1271   {
1272    cout<<endl;
1273    cout<<" WARNING (MH): "<<Form("!fRePEBE[%d]||!fImPEBE[%d]",sd,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1274    cout<<endl;
1275    exit(0);
1276   }
1277   if(!fReEtaEBE[sd]||!fImEtaEBE[sd])
1278   {
1279    cout<<endl;
1280    cout<<" WARNING (MH): "<<Form("!fReEtaEBE[%d]||!fImEtaEBE[%d]",sd,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1281    cout<<endl;
1282    exit(0);
1283   }
1284 for(Int_t fs=0;fs<2;fs++)
1285   {
1286    if(!fOverlapEBE[fs][sd]||!fOverlapEBE[fs][sd])
1287    {
1288     cout<<endl;
1289     cout<<" WARNING (MH): "<<Form("!fOverlapEBE[%d][%d]||!fOverlapEBE[%d][%d]",fs,sd,fs,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1290     cout<<endl;
1291     exit(0);
1292    }
1293    if(!fOverlapEBE2[fs][sd]||!fOverlapEBE2[fs][sd])
1294    {
1295     cout<<endl;
1296     cout<<" WARNING (MH): "<<Form("!fOverlapEBE2[%d][%d]||!fOverlapEBE2[%d][%d]",fs,sd,fs,sd)<<" is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1297     cout<<endl;
1298     exit(0);
1299    }
1300   } // end of for(Int_t fs=0;fs<2;fs++)
1301  } // end of for(Int_t sd=0;sd<2;sd++)
1302
1303  if(!f2pCorrelatorCosPsiDiff) {
1304   cout<<endl;
1305   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1306   cout<<endl;
1307   exit(0);
1308  }
1309  if(!f2pCorrelatorCosPsiSum) {
1310   cout<<endl;
1311   cout<<" WARNING (MH): f2pCorrelatorCosPsiSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1312   cout<<endl;
1313   exit(0);
1314  }
1315  if(!f2pCorrelatorSinPsiDiff) {
1316   cout<<endl;
1317   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiff is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1318   cout<<endl;
1319   exit(0);
1320  }
1321  if(!f2pCorrelatorSinPsiSum) {
1322   cout<<endl;
1323   cout<<" WARNING (MH): f2pCorrelatorSinPsiSum is NULL in CheckPointersUsedInMake() !!!!"<<endl;
1324   cout<<endl;
1325   exit(0);
1326  }
1327
1328 } // end of AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInMake()
1329
1330 //================================================================================================================
1331
1332 void AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInFinish()
1333 {
1334  // Check pointers used in method Finish().
1335
1336  if(!fAnalysisSettings)
1337  {
1338   cout<<endl;
1339   cout<<" WARNING (MH): fAnalysisSettings is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1340   cout<<endl;
1341   exit(0);
1342  }
1343  if(!f3pCorrelatorPro)
1344  {
1345   cout<<endl;
1346   cout<<" WARNING (MH): f3pCorrelatorPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1347   cout<<endl;
1348   exit(0);
1349  }
1350  if(!fNonIsotropicTermsPro)
1351  {
1352   cout<<endl;
1353   cout<<" WARNING (MH): fNonIsotropicTermsPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1354   cout<<endl;
1355   exit(0);
1356  }
1357  if(!f3pPOICorrelatorVsM && fCalculateVsM)
1358  {
1359   cout<<endl;
1360   cout<<" WARNING (MH): f3pPOICorrelatorVsM is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1361   cout<<endl;
1362   exit(0);
1363  }
1364  if(!f3pCorrelatorVsMPro && fCalculateVsM)
1365  {
1366   cout<<endl;
1367   cout<<" WARNING (MH): f3pCorrelatorVsMPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1368   cout<<endl;
1369   exit(0);
1370  }
1371  if(!f3pCorrelatorVsMHist && fCalculateVsM)
1372  {
1373   cout<<endl;
1374   cout<<" WARNING (MH): f3pCorrelatorVsMHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1375   cout<<endl;
1376   exit(0);
1377  }
1378  if(!fNonIsotropicTermsVsMPro && fCalculateVsM)
1379  {
1380   cout<<endl;
1381   cout<<" WARNING (MH): fNonIsotropicTermsVsMPro is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1382   cout<<endl;
1383   exit(0);
1384  }
1385
1386  if(!f3pCorrelatorHist)
1387  {
1388   cout<<endl;
1389   cout<<" WARNING (MH): f3pCorrelatorHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1390   cout<<endl;
1391   exit(0);
1392  }
1393  if(!fDetectorBiasHist)
1394  {
1395   cout<<endl;
1396   cout<<" WARNING (MH): fDetectorBiasHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1397   cout<<endl;
1398   exit(0);
1399  }
1400  /* to be improved - enabled eventually
1401  if(!fDetectorBiasVsMHist && fCalculateVsM)
1402  {
1403   cout<<endl;
1404   cout<<" WARNING (MH): !fDetectorBiasVsMHist is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1405   cout<<endl;
1406   exit(0);
1407  }
1408  */
1409
1410  // Differential correlators:
1411  if(!fEvaluateDifferential3pCorrelator){return;}
1412  for(Int_t sd=0;sd<2;sd++)
1413  {
1414   if(!(f3pCorrelatorVsPtSumDiffPro[sd]))
1415   {
1416    cout<<endl;
1417    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsPtSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1418    cout<<endl;
1419    exit(0);
1420   }
1421   if(!(f3pCorrelatorVsEtaSumDiffPro[sd]))
1422   {
1423    cout<<endl;
1424    cout<<" WARNING (MH): "<<Form("f3pCorrelatorVsEtaSumDiffPro[%d]",sd)<<" is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1425    cout<<endl;
1426    exit(0);
1427   }
1428  }
1429
1430  if(!f2pCorrelatorCosPsiDiff) {
1431   cout<<endl;
1432   cout<<" WARNING (MH): f2pCorrelatorCosPsiDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1433   cout<<endl;
1434   exit(0);
1435  }
1436  if(!f2pCorrelatorCosPsiSum) {
1437   cout<<endl;
1438   cout<<" WARNING (MH): f2pCorrelatorCosPsiSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1439   cout<<endl;
1440   exit(0);
1441  }
1442  if(!f2pCorrelatorSinPsiDiff) {
1443   cout<<endl;
1444   cout<<" WARNING (MH): f2pCorrelatorSinPsiDiff is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1445   cout<<endl;
1446   exit(0);
1447  }
1448  if(!f2pCorrelatorSinPsiSum) {
1449   cout<<endl;
1450   cout<<" WARNING (MH): f2pCorrelatorSinPsiSum is NULL in CheckPointersUsedInFinish() !!!!"<<endl;
1451   cout<<endl;
1452   exit(0);
1453  }
1454
1455 } // end of AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInFinish()
1456
1457 //================================================================================================================
1458
1459 void AliFlowAnalysisWithMixedHarmonics::PrintOnTheScreen()
1460 {
1461  // Print the final results on the screen.
1462
1463  // <<cos[n*(phi1+phi2-2phi3)]>>:
1464  Double_t d3pCorrelator = 0.;
1465  Double_t d3pCorrelatorError = 0.;
1466  if(!fCorrectForDetectorEffects)
1467  {
1468   d3pCorrelator = f3pCorrelatorPro->GetBinContent(1);
1469   d3pCorrelatorError = f3pCorrelatorPro->GetBinError(1);
1470  } else
1471    {
1472     d3pCorrelator = f3pCorrelatorHist->GetBinContent(1);
1473     d3pCorrelatorError = f3pCorrelatorHist->GetBinError(1);
1474    }
1475
1476  // <<cos[n*(psi1+psi2-2phi3)]>>:
1477  Double_t d3pCorrelatorPoi = 0.;
1478  Double_t d3pCorrelatorPoiError = 0.;
1479  if(fEvaluateDifferential3pCorrelator)
1480  {
1481   GetCorrelatorAndError(f3pCorrelatorVsPtSumDiffPro[0],
1482                        d3pCorrelatorPoi,
1483                        d3pCorrelatorPoiError);
1484  }
1485  cout<<endl;
1486  cout<<"*******************************************************"<<endl;
1487  cout<<"*******************************************************"<<endl;
1488  cout<<"                    Mixed Harmonics                      "<<endl;
1489  cout<<endl;
1490  if(fHarmonic!=1)
1491  {
1492   cout<<"  cos["<<fHarmonic<<"(phi1+phi2-2phi3)] = "<<d3pCorrelator<<" +/- "<<d3pCorrelatorError<<endl;
1493   cout<<"  cos["<<fHarmonic<<"(psi1+psi2-2phi3)] = "<<d3pCorrelatorPoi<<" +/- "<<d3pCorrelatorPoiError<<endl;
1494  } else
1495    {
1496     cout<<"  cos(phi1+phi2-2phi3) = "<<d3pCorrelator<<" +/- "<<d3pCorrelatorError<<endl;
1497     cout<<"  cos(psi1+psi2-2phi3) = "<<d3pCorrelatorPoi<<" +/- "<<d3pCorrelatorPoiError<<endl;
1498    }
1499  if(!fCorrectForDetectorEffects)
1500  {
1501   cout<<"  Detector Bias = "<<fDetectorBiasHist->GetBinContent(1)<<" (not corrected for)"<<endl;
1502  } else
1503    {
1504     cout<<"  Detector Bias = "<<fDetectorBiasHist->GetBinContent(1)<<" (corrected for)"<<endl;
1505    }
1506  cout<<endl;
1507  cout<<"             nEvts = "<<(Int_t)fCommonHists->GetHistMultRP()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()<<endl;
1508  cout<<"*******************************************************"<<endl;
1509  cout<<"*******************************************************"<<endl;
1510
1511 } // end of void AliFlowAnalysisWithMixedHarmonics::PrintOnTheScreen()
1512
1513 //================================================================================================================
1514
1515 void AliFlowAnalysisWithMixedHarmonics::AccessSettings()
1516 {
1517  // Access the settings for analysis with mixed harmonics.
1518
1519  fCorrectForDetectorEffects = (Bool_t)fAnalysisSettings->GetBinContent(1);
1520  fNoOfMultipicityBins = (Int_t)fAnalysisSettings->GetBinContent(2);
1521  fMultipicityBinWidth = (Double_t)fAnalysisSettings->GetBinContent(3);
1522  fMinMultiplicity = (Double_t)fAnalysisSettings->GetBinContent(4);
1523  fPrintOnTheScreen = (Bool_t)fAnalysisSettings->GetBinContent(5);
1524  fHarmonic = (Int_t)fAnalysisSettings->GetBinContent(6);
1525  fOppositeChargesPOI = (Bool_t)fAnalysisSettings->GetBinContent(7);
1526  fEvaluateDifferential3pCorrelator = (Bool_t)fAnalysisSettings->GetBinContent(8);
1527  fCalculateVsM = (Bool_t)fAnalysisSettings->GetBinContent(9);
1528  fShowBinLabelsVsM = (Bool_t)fAnalysisSettings->GetBinContent(10);
1529
1530 } // end of AliFlowAnalysisWithMixedHarmonics::AccessSettings()
1531
1532 //================================================================================================================
1533
1534 void AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffects()
1535 {
1536  // Correct measured 3-p correlator cos[n(phi1+phi2-2phi3)] for detector effects.
1537
1538  Double_t measured3pCorrelator = f3pCorrelatorPro->GetBinContent(1); // biased by detector effects
1539  Double_t corrected3pCorrelator = 0.; // corrected for detector effects
1540  Double_t nonIsotropicTerms[10] = {0.}; // there are 10 distinct non-isotropic terms
1541  for(Int_t nit=0;nit<10;nit++)
1542  {
1543   nonIsotropicTerms[nit] = fNonIsotropicTermsPro->GetBinContent(nit+1);
1544  }
1545  // Calculate corrected 3-p correlator:
1546  corrected3pCorrelator = measured3pCorrelator
1547                        - nonIsotropicTerms[2]*nonIsotropicTerms[4]
1548                        - nonIsotropicTerms[3]*nonIsotropicTerms[5]
1549                        - 2.*nonIsotropicTerms[0]*nonIsotropicTerms[6]
1550                        - 2.*nonIsotropicTerms[1]*nonIsotropicTerms[7]
1551                        + 2.*nonIsotropicTerms[2]*(pow(nonIsotropicTerms[0],2.)-pow(nonIsotropicTerms[1],2.))
1552                        + 4.*nonIsotropicTerms[3]*nonIsotropicTerms[0]*nonIsotropicTerms[1];
1553  // Store corrected correlator:
1554  if(fCorrectForDetectorEffects)
1555  {
1556   f3pCorrelatorHist->SetBinContent(1,corrected3pCorrelator);
1557   f3pCorrelatorHist->SetBinError(1,f3pCorrelatorPro->GetBinError(1)); // to be improved (propagate error for non-isotropic terms)
1558  }
1559  // Quantify bias from detector inefficiences to 3-p correlator. Remark: Bias is quantified as a
1560  // ratio between corrected and measured 3-p correlator:
1561  //              bias = corrected/measured
1562  // This bias is stored in histogram fDetectorBias.
1563  Double_t bias = 0.;
1564  if(TMath::Abs(measured3pCorrelator)>1.e-44)
1565  {
1566   bias = corrected3pCorrelator/measured3pCorrelator;
1567   fDetectorBiasHist->SetBinContent(1,bias);
1568  }
1569
1570 } // end of AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffects()
1571
1572 //================================================================================================================
1573
1574 void AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffectsVsM()
1575 {
1576  // Correct measured 3-p correlator cos[n(phi1+phi2-2phi3)] vs M for detector effects.
1577
1578  for(Int_t b=1;b<=fNoOfMultipicityBins+2;b++)
1579  {
1580   Double_t measured3pCorrelator = f3pCorrelatorVsMPro->GetBinContent(b); // biased by detector effects
1581   Double_t corrected3pCorrelator = 0.; // corrected for detector effects
1582   Double_t nonIsotropicTerms[10] = {0.}; // there are 10 distinct non-isotropic terms
1583   for(Int_t nit=0;nit<10;nit++)
1584   {
1585    nonIsotropicTerms[nit] = fNonIsotropicTermsVsMPro->GetBinContent(fNonIsotropicTermsVsMPro->GetBin(nit+1,b));
1586   }
1587   // Calculate corrected 3-p correlator:
1588   corrected3pCorrelator = measured3pCorrelator
1589                         - nonIsotropicTerms[2]*nonIsotropicTerms[4]
1590                         - nonIsotropicTerms[3]*nonIsotropicTerms[5]
1591                         - 2.*nonIsotropicTerms[0]*nonIsotropicTerms[6]
1592                         - 2.*nonIsotropicTerms[1]*nonIsotropicTerms[7]
1593                         + 2.*nonIsotropicTerms[2]*(pow(nonIsotropicTerms[0],2.)-pow(nonIsotropicTerms[1],2.))
1594                         + 4.*nonIsotropicTerms[3]*nonIsotropicTerms[0]*nonIsotropicTerms[1];
1595   // Store corrected correlator:
1596   if(fCorrectForDetectorEffects)
1597   {
1598    f3pCorrelatorVsMHist->SetBinContent(b,corrected3pCorrelator);
1599    f3pCorrelatorVsMHist->SetBinError(b,f3pCorrelatorVsMPro->GetBinError(b)); // to be improved (propagate error for non-isotropic terms)
1600   }
1601   // Quantify bias from detector inefficiences to 3-p correlator. Remark: Bias is quantified as a
1602   // ratio between corrected and measured 3-p correlator:
1603   //              bias = corrected/measured
1604   // This bias is stored in histogram fDetectorBias.
1605   Double_t bias = 0.;
1606   if(measured3pCorrelator)
1607   {
1608    bias = corrected3pCorrelator/measured3pCorrelator;
1609    fDetectorBiasVsMHist->SetBinContent(b,bias);
1610   }
1611  } // end of for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
1612
1613 } // end of AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffectsVsM()
1614
1615 //================================================================================================================
1616
1617 void AliFlowAnalysisWithMixedHarmonics::ResetEventByEventQuantities()
1618 {
1619  // Reset all event-by-event quantities.
1620
1621  fReQnk->Zero();
1622  fImQnk->Zero();
1623  fSpk->Zero();
1624
1625  if(!fEvaluateDifferential3pCorrelator){return;}
1626  for(Int_t sd=0;sd<2;sd++)
1627  {
1628   fRePEBE[sd]->Reset();
1629   fImPEBE[sd]->Reset();
1630   fReEtaEBE[sd]->Reset();
1631   fImEtaEBE[sd]->Reset();
1632  }
1633  for(Int_t fs=0;fs<2;fs++)
1634  {
1635   for(Int_t sd=0;sd<2;sd++)
1636   {
1637    fOverlapEBE[fs][sd]->Reset();
1638    fOverlapEBE2[fs][sd]->Reset();
1639   }
1640  }
1641
1642 } // end of void AliFlowAnalysisWithMixedHarmonics::ResetEventByEventQuantities()
1643
1644 //================================================================================================================
1645
1646 void AliFlowAnalysisWithMixedHarmonics::Calculate3pCorrelator()
1647 {
1648  // Calculate 3-p azimuthal correlator cos[n(phi1+phi2-2phi3)] in terms of Q_{n,k} and S_{p,k}.
1649
1650  // a) Calculate 3-p correlator without using particle weights;
1651  // b) Calculate 3-p correlator with using particle weights.
1652
1653  // a) Calculate 3-p correlator without using particle weights:
1654  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1655  {
1656   // Multiplicity (number of RPs):
1657   Double_t dMult = (*fSpk)(0,0);
1658   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n and 2n:
1659   Double_t dReQ1n = (*fReQnk)(0,0);
1660   Double_t dReQ2n = (*fReQnk)(1,0);
1661   Double_t dImQ1n = (*fImQnk)(0,0);
1662   Double_t dImQ2n = (*fImQnk)(1,0);
1663   // 3-particle azimuthal correlator <cos(n*(phi1+phi2-2phi3))>:
1664   Double_t three1n1n2n = (pow(dReQ1n,2.)*dReQ2n + 2.*dReQ1n*dImQ1n*dImQ2n - pow(dImQ1n,2.)*dReQ2n
1665                        - 2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
1666                        - (pow(dReQ2n,2.)+pow(dImQ2n,2.))+2.*dMult)
1667                        / (dMult*(dMult-1.)*(dMult-2.));
1668
1669   // Fill all-events profile:
1670   f3pCorrelatorPro->Fill(0.5,three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
1671
1672   // 3-particle azimuthal correlator <cos(n*(phi1+phi2-2phi3))> vs multiplicity:
1673   if(fCalculateVsM)
1674   {
1675    if(dMult<fMinMultiplicity)
1676    {
1677     f3pCorrelatorVsMPro->Fill(0.5,three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
1678    } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
1679      {
1680       f3pCorrelatorVsMPro->Fill(0.5+fNoOfMultipicityBins+1,three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
1681      } else
1682        {
1683         f3pCorrelatorVsMPro->Fill(1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),three1n1n2n,dMult*(dMult-1.)*(dMult-2.));
1684        }
1685   } // end of if(fCalculateVsM)
1686
1687  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1688
1689  // b) Calculate 3-p correlator with using particle weights:
1690  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1691  {
1692   // ...
1693  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1694
1695 } // end of void AliFlowAnalysisWithMixedHarmonics::Calculate3pCorrelator()
1696
1697 //================================================================================================================
1698
1699 void AliFlowAnalysisWithMixedHarmonics::Calculate5pCorrelator()
1700 {
1701  // Calculate 5-p azimuthal correlator cos[n(2phi1+2phi2+2phi3-3phi4-3phi5)] in terms of Q_{n,k} and S_{p,k}.
1702
1703  // a) Calculate 5-p correlator without using particle weights;
1704  // b) Calculate 5-p correlator with using particle weights.
1705
1706  // a) Calculate 5-p correlator without using particle weights:
1707  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1708  {
1709   // Multiplicity (number of RPs):
1710   Double_t dMult = (*fSpk)(0,0);
1711   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n,2n,...,6n:
1712   Double_t dReQ1n = (*fReQnk)(0,0);
1713   Double_t dReQ2n = (*fReQnk)(1,0);
1714   Double_t dReQ3n = (*fReQnk)(2,0);
1715   Double_t dReQ4n = (*fReQnk)(3,0);
1716   //Double_t dReQ5n = (*fReQnk)(4,0); // not needed
1717   Double_t dReQ6n = (*fReQnk)(5,0);
1718   Double_t dImQ1n = (*fImQnk)(0,0);
1719   Double_t dImQ2n = (*fImQnk)(1,0);
1720   Double_t dImQ3n = (*fImQnk)(2,0);
1721   Double_t dImQ4n = (*fImQnk)(3,0);
1722   //Double_t dImQ5n = (*fImQnk)(4,0); // not needed
1723   Double_t dImQ6n = (*fImQnk)(5,0);
1724
1725   // 5-particle azimuthal correlator:
1726   Double_t five2n2n2n3n3n = 0.; // <cos[n(2phi1+2phi2+2phi3-3phi4-3phi5)]>
1727   Double_t reQ2nQ2nQ2nQ3nstarQ3nstar = pow(dReQ2n,3.)*pow(dReQ3n,2.)
1728                                      - 3.*dReQ2n*pow(dReQ3n,2.)*pow(dImQ2n,2.)
1729                                      + 6.*pow(dReQ2n,2.)*dReQ3n*dImQ2n*dImQ3n
1730                                      - 2.*dReQ3n*pow(dImQ2n,3.)*dImQ3n-pow(dReQ2n,3.)*pow(dImQ3n,2.)
1731                                      + 3.*dReQ2n*pow(dImQ2n,2.)*pow(dImQ3n,2.);
1732   Double_t reQ2nQ2nQ2nQ6nstar = dReQ6n*pow(dReQ2n,3)-3.*dReQ2n*dReQ6n*pow(dImQ2n,2)
1733                                + 3.*dImQ2n*dImQ6n*pow(dReQ2n,2)-dImQ6n*pow(dImQ2n,3);
1734   Double_t reQ4nQ2nQ3nstarQ3nstar = (dReQ4n*dReQ2n-dImQ4n*dImQ2n)*(dReQ3n*dReQ3n-dImQ3n*dImQ3n)
1735                                   + 2.*(dReQ4n*dImQ2n+dImQ4n*dReQ2n)*dReQ3n*dImQ3n;
1736   Double_t reQ2nQ2nQ1nstarQ3nstar = (pow(dReQ2n,2.)-pow(dImQ2n,2.))*(dReQ3n*dReQ1n-dImQ3n*dImQ1n)
1737                                   + 2.*dReQ2n*dImQ2n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
1738   Double_t reQ6nQ3nstarQ3nstar = pow(dReQ3n,2.)*dReQ6n + 2.*dReQ3n*dImQ3n*dImQ6n
1739                                - pow(dImQ3n,2.)*dReQ6n;
1740   Double_t reQ4nQ2nQ6nstar = dReQ6n*dReQ4n*dReQ2n-dReQ6n*dImQ4n*dImQ2n+dImQ6n*dReQ4n*dImQ2n
1741                            + dImQ6n*dImQ4n*dReQ2n;
1742   Double_t reQ4nQ1nstarQ3nstar = dReQ4n*(dReQ3n*dReQ1n-dImQ3n*dImQ1n)+dImQ4n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
1743   Double_t reQ2nQ2nQ4nstar = pow(dReQ2n,2.)*dReQ4n+2.*dReQ2n*dImQ2n*dImQ4n-pow(dImQ2n,2.)*dReQ4n;
1744   Double_t reQ2nQ1nQ3nstar = dReQ3n*dReQ2n*dReQ1n-dReQ3n*dImQ2n*dImQ1n+dImQ3n*dReQ2n*dImQ1n
1745                            + dImQ3n*dImQ2n*dReQ1n;
1746   Double_t reQ2nQ1nstarQ1nstar = pow(dReQ1n,2.)*dReQ2n + 2.*dReQ1n*dImQ1n*dImQ2n - pow(dImQ1n,2.)*dReQ2n;
1747   // Analytic expression for 5-particle azimuthal correlator:
1748   five2n2n2n3n3n = (reQ2nQ2nQ2nQ3nstarQ3nstar-reQ2nQ2nQ2nQ6nstar-3.*reQ4nQ2nQ3nstarQ3nstar
1749                  - 6.*reQ2nQ2nQ1nstarQ3nstar+2.*reQ6nQ3nstarQ3nstar+3.*reQ4nQ2nQ6nstar
1750                  + 6.*reQ4nQ1nstarQ3nstar+6.*reQ2nQ2nQ4nstar
1751                  + 12.*reQ2nQ1nQ3nstar+6.*reQ2nQ1nstarQ1nstar
1752                  - 2.*((pow(dReQ6n,2.)+pow(dImQ6n,2.))
1753                  + 3.*(pow(dReQ4n,2.)+pow(dImQ4n,2.))
1754                  + 6.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
1755                  + 9.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
1756                  + 6.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))-12.*dMult))
1757                  /(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
1758   // Fill all-events profile:
1759   f5pCorrelatorPro->Fill(0.5,five2n2n2n3n3n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
1760  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1761
1762  // b) Calculate 5-p correlator with using particle weights:
1763  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1764  {
1765   // ...
1766  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1767
1768 } // end of void AliFlowAnalysisWithMixedHarmonics::Calculate5pCorrelator()
1769
1770 //================================================================================================================
1771
1772 void AliFlowAnalysisWithMixedHarmonics::CalculateNonIsotropicTerms()
1773 {
1774  // Calculate non-isotropic terms which appear in the decomposition of 3-p correlator <cos[n(phi1+phi2-2phi3)]>.
1775
1776  // a) Calculate without using particle weights;
1777  // b) Calculate using particle weights.
1778
1779  // For detector with uniform acceptance all these terms vanish. These non-isotropic terms are stored in fNonIsotropicTermsPro.
1780  // Binning of fNonIsotropicTermsPro is organized as follows:
1781  //  1st bin: <<cos(n*phi1)>>
1782  //  2nd bin: <<sin(n*phi1)>>
1783  //  3rd bin: <<cos(2n*phi1)>>
1784  //  4th bin: <<sin(2n*phi1)>>
1785  //  5th bin: <<cos(n*(phi1+phi2)>>
1786  //  6th bin: <<sin(n*(phi1+phi2)>>
1787  //  7th bin: <<cos(n*(2phi1-phi2)>>
1788  //  8th bin: <<sin(n*(2phi1-phi2)>>
1789  //  9th bin: <<cos(n*(phi1-phi2-phi3)>> // not needed
1790  // 10th bin: <<sin(n*(phi1-phi2-phi3)>> // not needed
1791
1792  // a) Calculate without using particle weights:
1793  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1794  {
1795   // Multiplicity (number of RPs):
1796   Double_t dMult = (*fSpk)(0,0);
1797   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n and 2n:
1798   Double_t dReQ1n = (*fReQnk)(0,0);
1799   Double_t dReQ2n = (*fReQnk)(1,0);
1800   Double_t dImQ1n = (*fImQnk)(0,0);
1801   Double_t dImQ2n = (*fImQnk)(1,0);
1802   // 1-particle terms:
1803   Double_t cosP1n = 0.; // <cos(n*(phi1))>
1804   Double_t sinP1n = 0.; // <sin(n*(phi1))>
1805   Double_t cosP2n = 0.; // <cos(2n*(phi1))>
1806   Double_t sinP2n = 0.; // <sin(2n*(phi1))>
1807   if(dMult>0)
1808   {
1809    cosP1n = dReQ1n/dMult;
1810    sinP1n = dImQ1n/dMult;
1811    cosP2n = dReQ2n/dMult;
1812    sinP2n = dImQ2n/dMult;
1813    // All-event avarages:
1814    fNonIsotropicTermsPro->Fill(0.5,cosP1n,dMult); // <<cos(n*(phi1))>>
1815    fNonIsotropicTermsPro->Fill(1.5,sinP1n,dMult); // <<sin(n*(phi1))>>
1816    fNonIsotropicTermsPro->Fill(2.5,cosP2n,dMult); // <<cos(2n*(phi1))>>
1817    fNonIsotropicTermsPro->Fill(3.5,sinP2n,dMult); // <<sin(2n*(phi1))>>
1818    // All-event avarages vs M:
1819    if(fCalculateVsM)
1820    {
1821     if(dMult<fMinMultiplicity)
1822     {
1823      fNonIsotropicTermsVsMPro->Fill(0.5,0.5,cosP1n,dMult); // <<cos(n*(phi1))>>
1824      fNonIsotropicTermsVsMPro->Fill(1.5,0.5,sinP1n,dMult); // <<sin(n*(phi1))>>
1825      fNonIsotropicTermsVsMPro->Fill(2.5,0.5,cosP2n,dMult); // <<cos(2n*(phi1))>>
1826      fNonIsotropicTermsVsMPro->Fill(3.5,0.5,sinP2n,dMult); // <<sin(2n*(phi1))>>
1827     } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
1828       {
1829        fNonIsotropicTermsVsMPro->Fill(0.5,0.5+fNoOfMultipicityBins+1,cosP1n,dMult); // <<cos(n*(phi1))>>
1830        fNonIsotropicTermsVsMPro->Fill(1.5,0.5+fNoOfMultipicityBins+1,sinP1n,dMult); // <<sin(n*(phi1))>>
1831        fNonIsotropicTermsVsMPro->Fill(2.5,0.5+fNoOfMultipicityBins+1,cosP2n,dMult); // <<cos(2n*(phi1))>>
1832        fNonIsotropicTermsVsMPro->Fill(3.5,0.5+fNoOfMultipicityBins+1,sinP2n,dMult); // <<sin(2n*(phi1))>>
1833       } else
1834         {
1835          fNonIsotropicTermsVsMPro->Fill(0.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP1n,dMult); // <<cos(n*(phi1))>>
1836          fNonIsotropicTermsVsMPro->Fill(1.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP1n,dMult); // <<sin(n*(phi1))>>
1837          fNonIsotropicTermsVsMPro->Fill(2.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP2n,dMult); // <<cos(2n*(phi1))>>
1838          fNonIsotropicTermsVsMPro->Fill(3.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP2n,dMult); // <<sin(2n*(phi1))>>
1839         }
1840    } // end of if(fCalculateVsM)
1841   } // end of if(dMult>0)
1842   // 2-particle terms:
1843   Double_t cosP1nP1n = 0.; // <cos(n*(phi1+phi2))>
1844   Double_t sinP1nP1n = 0.; // <sin(n*(phi1+phi2))>
1845   Double_t cosP2nM1n = 0.; // <cos(n*(2phi1-phi2))>
1846   Double_t sinP2nM1n = 0.; // <sin(n*(2phi1-phi2))>
1847   if(dMult>1)
1848   {
1849    cosP1nP1n = (pow(dReQ1n,2)-pow(dImQ1n,2)-dReQ2n)/(dMult*(dMult-1));
1850    sinP1nP1n = (2.*dReQ1n*dImQ1n-dImQ2n)/(dMult*(dMult-1));
1851    cosP2nM1n = (dReQ2n*dReQ1n+dImQ2n*dImQ1n-dReQ1n)/(dMult*(dMult-1));
1852    sinP2nM1n = (dImQ2n*dReQ1n-dReQ2n*dImQ1n-dImQ1n)/(dMult*(dMult-1));
1853    // All-event avarages:
1854    fNonIsotropicTermsPro->Fill(4.5,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
1855    fNonIsotropicTermsPro->Fill(5.5,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
1856    fNonIsotropicTermsPro->Fill(6.5,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
1857    fNonIsotropicTermsPro->Fill(7.5,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
1858    // All-event avarages vs M:
1859    if(fCalculateVsM)
1860    {
1861     if(dMult<fMinMultiplicity)
1862     {
1863      fNonIsotropicTermsVsMPro->Fill(4.5,0.5,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
1864      fNonIsotropicTermsVsMPro->Fill(5.5,0.5,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
1865      fNonIsotropicTermsVsMPro->Fill(6.5,0.5,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
1866      fNonIsotropicTermsVsMPro->Fill(7.5,0.5,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
1867     } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
1868       {
1869        fNonIsotropicTermsVsMPro->Fill(4.5,0.5+fNoOfMultipicityBins+1,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
1870        fNonIsotropicTermsVsMPro->Fill(5.5,0.5+fNoOfMultipicityBins+1,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
1871        fNonIsotropicTermsVsMPro->Fill(6.5,0.5+fNoOfMultipicityBins+1,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
1872        fNonIsotropicTermsVsMPro->Fill(7.5,0.5+fNoOfMultipicityBins+1,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
1873       } else
1874         {
1875          fNonIsotropicTermsVsMPro->Fill(4.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
1876          fNonIsotropicTermsVsMPro->Fill(5.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
1877          fNonIsotropicTermsVsMPro->Fill(6.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
1878          fNonIsotropicTermsVsMPro->Fill(7.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
1879         }
1880    } // end of if(fCalculateVsM)
1881   } // end of if(dMult>1)
1882   // 3-particle: correct and ready but not needed, hence commented out.
1883   /*
1884   Double_t cosP1nM1nM1n = 0.; // <cos(n*(phi1-phi2-phi3))>
1885   Double_t sinP1nM1nM1n = 0.; // <sin(n*(phi1-phi2-phi3))>
1886   if(dMult>2)
1887   {
1888    cosP1nM1nM1n = (dReQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))-dReQ1n*dReQ2n-dImQ1n*dImQ2n-2.*(dMult-1)*dReQ1n)
1889                 / (dMult*(dMult-1)*(dMult-2));
1890    sinP1nM1nM1n = (-dImQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))+dReQ1n*dImQ2n-dImQ1n*dReQ2n+2.*(dMult-1)*dImQ1n)
1891                 / (dMult*(dMult-1)*(dMult-2));
1892    // All-events avarages:
1893    fNonIsotropicTermsPro->Fill(8.5,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
1894    fNonIsotropicTermsPro->Fill(9.5,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
1895    // All-events avarages vs M:
1896    if(fCalculateVsM)
1897    {
1898     if(dMult<fMinMultiplicity)
1899     {
1900      fNonIsotropicTermsVsMPro->Fill(8.5,0.5,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
1901      fNonIsotropicTermsVsMPro->Fill(9.5,0.5,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
1902     } else if(dMult>=fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)
1903       {
1904        fNonIsotropicTermsVsMPro->Fill(8.5,0.5+fNoOfMultipicityBins+1,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
1905        fNonIsotropicTermsVsMPro->Fill(9.5,0.5+fNoOfMultipicityBins+1,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
1906       } else
1907         {
1908          fNonIsotropicTermsVsMPro->Fill(8.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),
1909                                         cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
1910          fNonIsotropicTermsVsMPro->Fill(9.5,1.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),
1911                                         sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
1912         }
1913    } // end of if(fCalculateVsM)
1914   } // end of if(dMult>2)
1915   */
1916  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1917
1918  // b) Calculate using particle weights:
1919  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1920  {
1921   // ...
1922  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1923
1924 } // end of void AliFlowAnalysisWithMixedHarmonics::CalculateNonIsotropicTerms()
1925
1926 //================================================================================================================
1927
1928 void AliFlowAnalysisWithMixedHarmonics::CalculateDifferential3pCorrelator(Double_t &gIntegratedValue)
1929 {
1930  // Calculate differential 3-p azimuthal correlator cos[n(psi1+psi2-2phi3)] in terms of Q_{2n}, p_{n}, q1_{n} and q2_{n}.
1931
1932  // a) Calculate differential 3-p correlator without using particle weights;
1933  // b) Calculate differential 3-p correlator with using particle weights.
1934
1935  // a) Calculate differential 3-p correlator without using particle weights:
1936  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1937  {
1938    Int_t iBinCounter = 0;
1939    Double_t gSumBinContentTimesWeight = 0., gSumWeight = 0.;
1940    Double_t gSumBinContentTimesWeightSquared = 0.;
1941
1942   // Multiplicity (number of RPs):
1943   Double_t dMult = (*fSpk)(0,0);
1944   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonic 2n:
1945   Double_t dReQ2n = (*fReQnk)(1,0);
1946   Double_t dImQ2n = (*fImQnk)(1,0);
1947   for(Int_t sd=0;sd<2;sd++)
1948   {
1949     // [(p1+p2)/2,|p1-p2|]
1950    // looping over all bins and calculating reduced correlations:
1951    for(Int_t b=1;b<=fnBinsPt;b++)
1952    {
1953     // real and imaginary parts of p_{n}:
1954     Double_t p1nRe = fRePEBE[sd]->GetBinContent(b)*fRePEBE[sd]->GetBinEntries(b);
1955     Double_t p1nIm = fImPEBE[sd]->GetBinContent(b)*fImPEBE[sd]->GetBinEntries(b);
1956     // overlap 1: to be improved (terminology)
1957     Double_t overlap1 = fOverlapEBE[0][sd]->GetBinContent(b)*fOverlapEBE[0][sd]->GetBinEntries(b);
1958     // overlap 2: to be improved (terminology)
1959     Double_t overlap2 = fOverlapEBE[1][sd]->GetBinContent(b)*fOverlapEBE[1][sd]->GetBinEntries(b);
1960     // number of pairs of POIs in particular (p1+p2)/2 or |p1-p2| bin:
1961     Double_t mp = fRePEBE[sd]->GetBinEntries(b);
1962     // number of pairs of POI1/RP and POI2 in particular (p1+p2)/2 or |p1-p2| bin:
1963     Double_t mOverlap1 = fOverlapEBE[0][sd]->GetBinEntries(b);
1964     // number of pairs of POI2/RP and POI1 in particular (p1+p2)/2 or |p1-p2| bin:
1965     Double_t mOverlap2 = fOverlapEBE[1][sd]->GetBinEntries(b);
1966     // e-b-e weight for cos[n(psi1+psi2-2phi3)]:
1967     Double_t weight = mp*dMult-mOverlap1-mOverlap2;
1968
1969     Double_t cosP2nphi1M1npsi2M1npsi2 = 0; // cos[n(psi1+psi2-2phi3)]
1970     if(weight>0.)
1971     {
1972      cosP2nphi1M1npsi2M1npsi2 = (p1nRe*dReQ2n+p1nIm*dImQ2n-overlap1-overlap2)/(weight);
1973     }
1974     f3pCorrelatorVsPtSumDiffPro[sd]->Fill(fPtMin+(b-1)*fPtBinWidth,cosP2nphi1M1npsi2M1npsi2,weight);
1975     if(sd == 0) {
1976       iBinCounter += 1;
1977
1978       gSumBinContentTimesWeight += f3pCorrelatorVsPtSumDiffPro[sd]->GetBinContent(b)*f3pCorrelatorVsPtSumDiffPro[sd]->GetBinEntries(b);
1979       gSumWeight += f3pCorrelatorVsPtSumDiffPro[sd]->GetBinEntries(b);
1980       gSumBinContentTimesWeightSquared += TMath::Power(gSumBinContentTimesWeight,2);
1981     }
1982    } // end of for(Int_t b=1;b<=fnBinsPt;b++)
1983
1984     // [(eta1+eta2)/2,|eta1-eta2|]
1985    // looping over all bins and calculating reduced correlations:
1986    for(Int_t k=1;k<=fnBinsEta;k++)
1987    {
1988     // real and imaginary parts of p_{n}:
1989     Double_t p1nRe = fReEtaEBE[sd]->GetBinContent(k)*fReEtaEBE[sd]->GetBinEntries(k);
1990     Double_t p1nIm = fImEtaEBE[sd]->GetBinContent(k)*fImEtaEBE[sd]->GetBinEntries(k);
1991     // overlap 1: to be improved (terminology)
1992     Double_t overlap1 = fOverlapEBE2[0][sd]->GetBinContent(k)*fOverlapEBE2[0][sd]->GetBinEntries(k);
1993     // overlap 2: to be improved (terminology)
1994     Double_t overlap2 = fOverlapEBE2[1][sd]->GetBinContent(k)*fOverlapEBE2[1][sd]->GetBinEntries(k);
1995     // number of pairs of POIs in particular (eta1+eta2)/2 or |eta1-eta2| bin:
1996     Double_t mp = fReEtaEBE[sd]->GetBinEntries(k);
1997     // number of pairs of POI1/RP and POI2 in particular (eta1+eta2)/2 or |eta1-eta2| bin:
1998     Double_t mOverlap1 = fOverlapEBE2[0][sd]->GetBinEntries(k);
1999     // number of pairs of POI2/RP and POI1 in particular (eta1+eta2)/2 or |eta1-eta2| bin:
2000     Double_t mOverlap2 = fOverlapEBE2[1][sd]->GetBinEntries(k);
2001     // e-b-e weight for cos[n(psi1+psi2-2phi3)]:
2002     Double_t weight = mp*dMult-mOverlap1-mOverlap2;
2003
2004     Double_t cosP2nphi1M1npsi2M1npsi2 = 0; // cos[n(psi1+psi2-2phi3)]
2005     if(weight>0.)
2006     {
2007      cosP2nphi1M1npsi2M1npsi2 = (p1nRe*dReQ2n+p1nIm*dImQ2n-overlap1-overlap2)/(weight);
2008     }
2009     f3pCorrelatorVsEtaSumDiffPro[sd]->Fill(fEtaMin+(k-1)*fEtaBinWidth,cosP2nphi1M1npsi2M1npsi2,weight);
2010    } // end of for(Int_t k=1;k<=fnBinsEta;k++)
2011   } // end of for(Int_t sd=0;sd<2;sd++)
2012
2013   gIntegratedValue = -1000.;
2014   if((gSumWeight)&&(iBinCounter))
2015     gIntegratedValue = gSumBinContentTimesWeight/gSumWeight;
2016  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
2017
2018  // b) Calculate differential 3-p correlator by using particle weights:
2019  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2020  {
2021   // ...
2022  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
2023
2024 } // end of void AliFlowAnalysisWithMixedHarmonics::CalculateDifferential3pCorrelator()
2025
2026 //================================================================================================================
2027
2028 void AliFlowAnalysisWithMixedHarmonics::GetCorrelatorAndError(TProfile *g3pCorrelatorVsPt, Double_t &g3pCorrelatorValue, Double_t &g3pCorrelatorError) {
2029   //Retrieves the 3p correlator <<cos[n(psi1+psi2-2phi3)]>>
2030   //and its error
2031   Double_t gSumXi = 0.;
2032   Double_t gSumYi = 0.;
2033   Double_t gSumXiYi = 0.;
2034   Double_t gSumXiYi2 = 0.;
2035   Double_t gSumXi2Yi2 = 0.;
2036   Double_t gSumDeltaXi2 = 0.;
2037   Double_t gSumYi2DeltaXi2 = 0.;
2038
2039   for(Int_t iBin = 1; iBin <= g3pCorrelatorVsPt->GetNbinsX(); iBin++) {
2040     gSumXi += g3pCorrelatorVsPt->GetBinEntries(iBin);
2041     gSumYi += g3pCorrelatorVsPt->GetBinContent(iBin);
2042     gSumXiYi += g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin);
2043     gSumXiYi2 += g3pCorrelatorVsPt->GetBinEntries(iBin)*TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2);
2044     gSumXi2Yi2 += TMath::Power(g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin),2);
2045     gSumDeltaXi2 += TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);
2046     gSumYi2DeltaXi2 += TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2) + TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);
2047   }
2048
2049   g3pCorrelatorValue = -1000.;
2050   g3pCorrelatorError = 1000.;
2051
2052   if(gSumXi != 0.)
2053     g3pCorrelatorValue = gSumXiYi/gSumXi;
2054   if((gSumXi != 0.)&&(gSumXiYi != 0.))
2055     g3pCorrelatorError = TMath::Abs((gSumXiYi/gSumXi))*TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumYi2DeltaXi2)/gSumXiYi),2) + TMath::Power((gSumDeltaXi2/gSumXi),2));
2056 }
2057 //================================================================================================================
2058