AliFlowAnalysisWithQCumulants.cxx

Go to the documentation of this file.

/*************************************************************************
* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  * 
**************************************************************************/

/********************************** 
 * flow analysis with Q-cumulants * 
 *                                * 
 * author: Ante Bilandzic         * 
 *        (abilandzic@gmail.com)  *
 *********************************/ 

#define AliFlowAnalysisWithQCumulants_cxx

#include "Riostream.h"
#include "AliFlowCommonConstants.h"
#include "AliFlowCommonHist.h"
#include "AliFlowCommonHistResults.h"
#include "TChain.h"

#include "TFile.h"
#include "TList.h"
#include "TGraph.h"
#include "TParticle.h"
#include "TRandom3.h"
#include "TStyle.h"
#include "TProfile.h"
#include "TProfile2D.h" 
#include "TMath.h"
#include "TArrow.h"
#include "TPaveLabel.h"
#include "TCanvas.h"
#include "AliFlowEventSimple.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowAnalysisWithQCumulants.h"
#include "TArrayD.h"
#include "TRandom.h"
#include "TF1.h"

class TH1;
class TH2;
class TGraph;
class TPave;
class TLatex;
class TMarker;
class TRandom3;
class TObjArray;
class TList;
class TCanvas;
class TSystem;
class TROOT;
class AliFlowVector;
class TVector;

//================================================================================================================

using std::endl;
using std::cout;
using std::flush;
ClassImp(AliFlowAnalysisWithQCumulants)

AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants(): 
 // 0.) base:
 fHistList(NULL),
 // 1.) common:
 fBookOnlyBasicCCH(kTRUE),
 fCommonHists(NULL),
 fCommonHists2nd(NULL), 
 fCommonHists4th(NULL),
 fCommonHists6th(NULL),
 fCommonHists8th(NULL),
 fCommonHistsResults2nd(NULL),
 fCommonHistsResults4th(NULL),
 fCommonHistsResults6th(NULL),
 fCommonHistsResults8th(NULL),
 fnBinsPhi(0),
 fPhiMin(0),
 fPhiMax(0),
 fPhiBinWidth(0),
 fnBinsPt(0),
 fPtMin(0),
 fPtMax(0),
 fPtBinWidth(0),
 fnBinsEta(0),
 fEtaMin(0),
 fEtaMax(0),
 fEtaBinWidth(0),
 fCommonConstants(NULL),
 fFillMultipleControlHistograms(kFALSE),
 fHarmonic(2),
 fAnalysisLabel(NULL),
 fMaxCommonResultsHistogram(8),
 // 2a.) particle weights:
 fWeightsList(NULL),
 fUsePhiWeights(kFALSE),
 fUsePtWeights(kFALSE),
 fUseEtaWeights(kFALSE),
 fUseTrackWeights(kFALSE),
 fUseParticleWeights(NULL),
 fPhiWeights(NULL),
 fPtWeights(NULL),
 fEtaWeights(NULL),
 // 2b.) event weights:
 fMultiplicityWeight(NULL),
 fMultiplicityIs(AliFlowCommonConstants::kRP),
 // 3.) integrated flow:
 fIntFlowList(NULL), 
 fIntFlowProfiles(NULL),
 fIntFlowResults(NULL),
 fIntFlowAllCorrelationsVsM(NULL),
 fIntFlowFlags(NULL),
 fApplyCorrectionForNUA(kFALSE),  
 fApplyCorrectionForNUAVsM(kFALSE),
 fnBinsMult(10000),
 fMinMult(0.),  
 fMaxMult(10000.), 
 fPropagateErrorAlsoFromNIT(kFALSE), 
 fCalculateCumulantsVsM(kFALSE),
 fCalculateAllCorrelationsVsM(kFALSE), 
 fMinimumBiasReferenceFlow(kTRUE), 
 fForgetAboutCovariances(kFALSE), 
 fStorePhiDistributionForOneEvent(kFALSE),
 fExactNoRPs(0),
 fUse2DHistograms(kFALSE),
 fFillProfilesVsMUsingWeights(kTRUE),
 fUseQvectorTerms(kFALSE),
 fReQ(NULL),
 fImQ(NULL),
 fSpk(NULL),
 fIntFlowCorrelationsEBE(NULL),
 fIntFlowEventWeightsForCorrelationsEBE(NULL),
 fIntFlowCorrelationsAllEBE(NULL),
 fNumberOfRPsEBE(0.),
 fNumberOfPOIsEBE(0.),
 fReferenceMultiplicityEBE(0.),  
 fAvMultiplicity(NULL),
 fIntFlowCorrelationsPro(NULL),
 fIntFlowSquaredCorrelationsPro(NULL),
 fIntFlowCorrelationsAllPro(NULL),
 fIntFlowExtraCorrelationsPro(NULL),
 fIntFlowProductOfCorrelationsPro(NULL),
 fIntFlowProductOfCorrectionTermsForNUAPro(NULL),
 fIntFlowCorrelationsHist(NULL),
 fIntFlowCorrelationsAllHist(NULL),
 fIntFlowCovariances(NULL),
 fIntFlowSumOfProductOfEventWeights(NULL),
 fIntFlowCovariancesNUA(NULL),
 fIntFlowSumOfProductOfEventWeightsNUA(NULL),
 fIntFlowQcumulants(NULL),
 fIntFlowQcumulantsRebinnedInM(NULL), 
 fIntFlowQcumulantsErrorSquaredRatio(NULL), 
 fIntFlow(NULL),
 fIntFlowRebinnedInM(NULL),
 fIntFlowDetectorBias(NULL),
 // 4.) differential flow:
 fDiffFlowList(NULL),
 fDiffFlowProfiles(NULL),
 fDiffFlowResults(NULL),
 fDiffFlow2D(NULL),
 fDiffFlowFlags(NULL),
 fCalculateDiffFlow(kTRUE),
 fCalculate2DDiffFlow(kFALSE),
 fCalculateDiffFlowVsEta(kTRUE),
 // 5.) other differential correlators:
 fOtherDiffCorrelatorsList(NULL),
 // 6.) distributions:
 fDistributionsList(NULL),
 fDistributionsFlags(NULL),
 fStoreDistributions(kFALSE),
 fnBinsForCorrelations(10000),
 // 7.) various:
 fVariousList(NULL),
 fPhiDistributionForOneEvent(NULL),
 // 8.) debugging and cross-checking:
 fNestedLoopsList(NULL),
 fEvaluateIntFlowNestedLoops(kFALSE),
 fEvaluateDiffFlowNestedLoops(kFALSE),
 fMaxAllowedMultiplicity(10),
 fEvaluateNestedLoops(NULL),
 fIntFlowDirectCorrelations(NULL),
 fIntFlowExtraDirectCorrelations(NULL),
 fCrossCheckInPtBinNo(10),
 fCrossCheckInEtaBinNo(20),
 fNoOfParticlesInBin(NULL),
 fMixedHarmonicsNestedLoops(NULL),
 // 9.) mixed harmonics: 
 fMixedHarmonicsList(NULL),
 fMixedHarmonicsProfiles(NULL),
 fMixedHarmonicsResults(NULL),
 fMixedHarmonicsErrorPropagation(NULL),
 fMixedHarmonicsFlags(NULL),
 fCalculateMixedHarmonics(kFALSE),
 fCalculateMixedHarmonicsVsM(kFALSE),
 f2pCorrelations(NULL),
 f3pCorrelations(NULL),
 f4pCorrelations(NULL),
 f5pCorrelations(NULL),
 f6pCorrelations(NULL),
 f7pCorrelations(NULL),
 f8pCorrelations(NULL),
 f2pCumulants(NULL),
 f3pCumulants(NULL),
 f4pCumulants(NULL),
 f5pCumulants(NULL),
 f6pCumulants(NULL),
 f7pCumulants(NULL),
 f8pCumulants(NULL),
 fMixedHarmonicProductOfEventWeights(NULL),
 fMixedHarmonicProductOfCorrelations(NULL),
 // 10.) Control histograms:
 fControlHistogramsList(NULL), 
 fControlHistogramsFlags(NULL),
 fStoreControlHistograms(kFALSE),
 fCorrelationNoRPsVsRefMult(NULL), 
 fCorrelationNoPOIsVsRefMult(NULL),
 fCorrelationNoRPsVsNoPOIs(NULL),
 // 11.) Bootstrap:
 fBootstrapList(NULL),
 fBootstrapProfilesList(NULL),
 fBootstrapResultsList(NULL),
 fBootstrapFlags(NULL),
 fUseBootstrap(kFALSE),
 fUseBootstrapVsM(kFALSE),
 fnSubsamples(10),
 fRandom(NULL),
 fBootstrapCorrelations(NULL),
 fBootstrapCumulants(NULL)
 {
  // constructor  
  
  // base list to hold all output objects:
  fHistList = new TList();
  fHistList->SetName("cobjQC");
  fHistList->SetOwner(kTRUE);
  
  // list to hold histograms with phi, pt and eta weights:      
  fWeightsList = new TList();
  
  // multiplicity weight:
  fMultiplicityWeight = new TString("combinations");

  // analysis label;
  fAnalysisLabel = new TString();
      
  // initialize all arrays:  
  this->InitializeArraysForIntFlow();
  this->InitializeArraysForDiffFlow();
  this->InitializeArraysForDistributions();
  this->InitializeArraysForVarious();
  this->InitializeArraysForNestedLoops();
  this->InitializeArraysForMixedHarmonics();
  this->InitializeArraysForControlHistograms();
  this->InitializeArraysForBootstrap();
  
 } // end of constructor
 
//================================================================================================================  

AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()
{
 // destructor
 
 delete fHistList;

} // end of AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()

//================================================================================================================

void AliFlowAnalysisWithQCumulants::Init()
{
 // a) Cross check if the settings make sense before starting the QC adventure;
 // b) Access all common constants;
 // c) Book all objects;
 // d) Store flags for integrated and differential flow;
 // e) Store flags for distributions of corelations;
 // f) Store harmonic which will be estimated;
 // g) Store flags for mixed harmonics;
 // h) Store flags for control histograms;
 // i) Store bootstrap flags.
  
 //save old value and prevent histograms from being added to directory
 //to avoid name clashes in case multiple analaysis objects are used
 //in an analysis
 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 TH1::AddDirectory(kFALSE);
 
 // a) Cross check if the settings make sense before starting the QC adventure; 
 this->CrossCheckSettings();
 // b) Access all common constants and book a profile to hold them:
 this->CommonConstants("Init");
 // c) Book all objects:
 this->BookAndFillWeightsHistograms(); 
 this->BookAndNestAllLists();
 this->BookCommonHistograms();
 this->BookEverythingForIntegratedFlow(); 
 this->BookEverythingForDifferentialFlow(); 
 this->BookEverythingFor2DDifferentialFlow(); 
 this->BookEverythingForDistributions();
 this->BookEverythingForVarious();
 this->BookEverythingForNestedLoops();
 this->BookEverythingForMixedHarmonics();
 this->BookEverythingForControlHistograms();
 this->BookEverythingForBootstrap();

 // d) Store flags for integrated and differential flow:
 this->StoreIntFlowFlags();
 this->StoreDiffFlowFlags();
 // e) Store flags for distributions of corelations:
 this->StoreFlagsForDistributions();
 // f) Store harmonic which will be estimated:
 this->StoreHarmonic();
 // g) Store flags for mixed harmonics:
 this->StoreMixedHarmonicsFlags();
 // h) Store flags for control histograms:
 this->StoreControlHistogramsFlags();
 // i) Store bootstrap flags:
 this->StoreBootstrapFlags();

 TH1::AddDirectory(oldHistAddStatus);

} // end of void AliFlowAnalysisWithQCumulants::Init()

//================================================================================================================

void AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
{
 // Running over data only in this method.
 
 // a) Check all pointers used in this method;
 // b) Define local variables;
 // c) Fill the common control histograms and call the method to fill fAvMultiplicity;
 // d) Loop over data and calculate e-b-e quantities Q_{n,k}, S_{p,k} and s_{p,k};
 // e) Calculate the final expressions for S_{p,k} and s_{p,k} (important !!!!); 
 // f) Call the methods which calculate correlations for reference flow;
 // g) Call the methods which calculate correlations for differential flow;
 // h) Call the methods which calculate correlations for 2D differential flow;
 // i) Call the methods which calculate other differential correlators;
 // j) Distributions of correlations;
 // k) Store phi distribution for one event to illustrate flow;
 // l) Cross-check with nested loops correlators for reference flow;
 // m) Cross-check with nested loops correlators for differential flow;
 // n) Reset all event-by-event quantities (very important !!!!). 

 // a) Check all pointers used in this method:
 this->CheckPointersUsedInMake();
 
 // b) Define local variables:
 Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 Double_t dPt  = 0.; // transverse momentum
 Double_t dEta = 0.; // pseudorapidity
 Double_t wPhi = 1.; // phi weight
 Double_t wPt  = 1.; // pt weight
 Double_t wEta = 1.; // eta weight
 Double_t wTrack = 1.; // track weight
 Int_t nCounterNoRPs = 0; // needed only for shuffling
 fNumberOfRPsEBE = anEvent->GetNumberOfRPs(); // number of RPs (i.e. number of reference particles)
 if(fExactNoRPs > 0 && fNumberOfRPsEBE<fExactNoRPs){return;}
 fNumberOfPOIsEBE = anEvent->GetNumberOfPOIs(); // number of POIs (i.e. number of particles of interest)
 fReferenceMultiplicityEBE = anEvent->GetReferenceMultiplicity(); // reference multiplicity for current event
 //Printf("Reference multiplicity (QC): %.1f",fReferenceMultiplicityEBE);
 Double_t ptEta[2] = {0.,0.}; // 0 = dPt, 1 = dEta
  
 // c) Fill the common control histograms and call the method to fill fAvMultiplicity:
 this->FillCommonControlHistograms(anEvent);                                                               
 this->FillAverageMultiplicities((Int_t)(fNumberOfRPsEBE)); 
 if(fStoreControlHistograms){this->FillControlHistograms(anEvent);}                                                              
                                                                                                                                                                                                                                                                                        
 // d) Loop over data and calculate e-b-e quantities Q_{n,k}, S_{p,k} and s_{p,k}:
 Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks
 AliFlowTrackSimple *aftsTrack = NULL;
 Int_t n = fHarmonic; // shortcut for the harmonic 
 for(Int_t i=0;i<nPrim;i++) 
 { 
  if(fExactNoRPs > 0 && nCounterNoRPs>fExactNoRPs){continue;}
  aftsTrack=anEvent->GetTrack(i);
  if(aftsTrack)
  {
   if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())){continue;} // safety measure: consider only tracks which are RPs or POIs
   if(aftsTrack->InRPSelection()) // RP condition:
   {    
    nCounterNoRPs++;
    dPhi = aftsTrack->Phi();
    dPt  = aftsTrack->Pt();
    dEta = aftsTrack->Eta();
    if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi weight for this particle:
    {
     wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
    }
    if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt weight for this particle:
    {
     wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); 
    }              
    if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta weight for this particle: 
    {
     wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); 
    }      
    // Access track weight:
    if(fUseTrackWeights)
    {
     wTrack = aftsTrack->Weight(); 
    }
    // Calculate Re[Q_{m*n,k}] and Im[Q_{m*n,k}] for this event (m = 1,2,...,12, k = 0,1,...,8):
    for(Int_t m=0;m<12;m++) // to be improved - hardwired 6 
    {
     for(Int_t k=0;k<9;k++) // to be improved - hardwired 9
     {
      (*fReQ)(m,k)+=pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1)*n*dPhi); 
      (*fImQ)(m,k)+=pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1)*n*dPhi); 
     } 
    }
    // Calculate S_{p,k} for this event (Remark: final calculation of S_{p,k} follows after the loop over data bellow):
    for(Int_t p=0;p<8;p++)
    {
     for(Int_t k=0;k<9;k++)
     {     
      (*fSpk)(p,k)+=pow(wPhi*wPt*wEta*wTrack,k);
     }
    } 
    // Differential flow:
    if(fCalculateDiffFlow || fCalculate2DDiffFlow)
    {
     ptEta[0] = dPt; 
     ptEta[1] = dEta; 
     // Calculate r_{m*n,k} and s_{p,k} (r_{m,k} is 'p-vector' for RPs): 
     for(Int_t k=0;k<9;k++) // to be improved - hardwired 9
     {
      for(Int_t m=0;m<4;m++) // to be improved - hardwired 4
      {
       if(fCalculateDiffFlow)
       {
        for(Int_t pe=0;pe<1+(Int_t)fCalculateDiffFlowVsEta;pe++) // pt or eta
        {
         fReRPQ1dEBE[0][pe][m][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
         fImRPQ1dEBE[0][pe][m][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1.)*n*dPhi),1.);          
         if(m==0) // s_{p,k} does not depend on index m
         {
          fs1dEBE[0][pe][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k),1.);
         } // end of if(m==0) // s_{p,k} does not depend on index m
        } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
       } // end of if(fCalculateDiffFlow) 
       if(fCalculate2DDiffFlow)
       {
        fReRPQ2dEBE[0][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
        fImRPQ2dEBE[0][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1.)*n*dPhi),1.);      
        if(m==0) // s_{p,k} does not depend on index m
        {
         fs2dEBE[0][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k),1.);
        } // end of if(m==0) // s_{p,k} does not depend on index m
       } // end of if(fCalculate2DDiffFlow)
      } // end of for(Int_t m=0;m<4;m++) // to be improved - hardwired 4
     } // end of for(Int_t k=0;k<9;k++) // to be improved - hardwired 9
     // Checking if RP particle is also POI particle:      
     if(aftsTrack->InPOISelection())
     {
      // Calculate q_{m*n,k} and s_{p,k} ('q-vector' and 's' for RPs && POIs): 
      for(Int_t k=0;k<9;k++) // to be improved - hardwired 9
      {
       for(Int_t m=0;m<4;m++) // to be improved - hardwired 4
       {
        if(fCalculateDiffFlow)
        {
         for(Int_t pe=0;pe<1+(Int_t)fCalculateDiffFlowVsEta;pe++) // pt or eta
         {
          fReRPQ1dEBE[2][pe][m][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
          fImRPQ1dEBE[2][pe][m][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1.)*n*dPhi),1.);          
          if(m==0) // s_{p,k} does not depend on index m
          {
           fs1dEBE[2][pe][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k),1.);
          } // end of if(m==0) // s_{p,k} does not depend on index m
         } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
        } // end of if(fCalculateDiffFlow) 
        if(fCalculate2DDiffFlow)
        {
         fReRPQ2dEBE[2][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
         fImRPQ2dEBE[2][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1.)*n*dPhi),1.);      
         if(m==0) // s_{p,k} does not depend on index m
         {
          fs2dEBE[2][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k),1.);
         } // end of if(m==0) // s_{p,k} does not depend on index m
        } // end of if(fCalculate2DDiffFlow)
       } // end of for(Int_t m=0;m<4;m++) // to be improved - hardwired 4
      } // end of for(Int_t k=0;k<9;k++) // to be improved - hardwired 9    
     } // end of if(aftsTrack->InPOISelection())  
    } // end of if(fCalculateDiffFlow || fCalculate2DDiffFlow)         
   } // end of if(pTrack->InRPSelection())
   if(aftsTrack->InPOISelection())
   {
    dPhi = aftsTrack->Phi();
    dPt  = aftsTrack->Pt();
    dEta = aftsTrack->Eta();
    wPhi = 1.;
    wPt  = 1.;
    wEta = 1.;
    wTrack = 1.;
    if(fUsePhiWeights && fPhiWeights && fnBinsPhi && aftsTrack->InRPSelection()) // determine phi weight for POI && RP particle:
    {
     wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
    }
    if(fUsePtWeights && fPtWeights && fnBinsPt && aftsTrack->InRPSelection()) // determine pt weight for POI && RP particle:
    {
     wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); 
    }              
    if(fUseEtaWeights && fEtaWeights && fEtaBinWidth && aftsTrack->InRPSelection()) // determine eta weight for POI && RP particle: 
    {
     wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); 
    }      
    // Access track weight for POI && RP particle:
    if(aftsTrack->InRPSelection() && fUseTrackWeights)
    {
     wTrack = aftsTrack->Weight(); 
    }
    ptEta[0] = dPt;
    ptEta[1] = dEta;
    // Calculate p_{m*n,k} ('p-vector' for POIs): 
    for(Int_t k=0;k<9;k++) // to be improved - hardwired 9
    {
     for(Int_t m=0;m<4;m++) // to be improved - hardwired 4
     {
      if(fCalculateDiffFlow)
      {
       for(Int_t pe=0;pe<1+(Int_t)fCalculateDiffFlowVsEta;pe++) // pt or eta
       {
        fReRPQ1dEBE[1][pe][m][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
        fImRPQ1dEBE[1][pe][m][k]->Fill(ptEta[pe],pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1.)*n*dPhi),1.);          
       } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
      } // end of if(fCalculateDiffFlow) 
      if(fCalculate2DDiffFlow)
      {
       fReRPQ2dEBE[1][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
       fImRPQ2dEBE[1][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta*wTrack,k)*TMath::Sin((m+1.)*n*dPhi),1.);      
      } // end of if(fCalculate2DDiffFlow)
     } // end of for(Int_t m=0;m<4;m++) // to be improved - hardwired 4
    } // end of for(Int_t k=0;k<9;k++) // to be improved - hardwired 9    
   } // end of if(pTrack->InPOISelection())    
  } else // to if(aftsTrack)
    {
     printf("\n WARNING (QC): No particle (i.e. aftsTrack is a NULL pointer in AFAWQC::Make())!!!!\n\n");
    }
 } // end of for(Int_t i=0;i<nPrim;i++) 

 // e) Calculate the final expressions for S_{p,k} and s_{p,k} (important !!!!):
 for(Int_t p=0;p<8;p++)
 {
  for(Int_t k=0;k<9;k++)
  {
   (*fSpk)(p,k)=pow((*fSpk)(p,k),p+1);
   // ... for the time being s_{p,k} dosn't need higher powers, so no need to finalize it here ...
  } // end of for(Int_t k=0;k<9;k++)  
 } // end of for(Int_t p=0;p<8;p++)
 
 // f) Call the methods which calculate correlations for reference flow:
 if(!fEvaluateIntFlowNestedLoops)
 {
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   if(fNumberOfRPsEBE>1){this->CalculateIntFlowCorrelations();} // without using particle weights
  } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
    {
     if(fNumberOfRPsEBE>1){this->CalculateIntFlowCorrelationsUsingParticleWeights();} // with using particle weights   
    }        
  // Whether or not using particle weights the following is calculated in the same way:  
  if(fNumberOfRPsEBE>3){this->CalculateIntFlowProductOfCorrelations();}
  if(fNumberOfRPsEBE>1){this->CalculateIntFlowSumOfEventWeights();}
  if(fNumberOfRPsEBE>1){this->CalculateIntFlowSumOfProductOfEventWeights();}  
  // Non-isotropic terms:
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   if(fNumberOfRPsEBE>0){this->CalculateIntFlowCorrectionsForNUASinTerms();}
   if(fNumberOfRPsEBE>0){this->CalculateIntFlowCorrectionsForNUACosTerms();}
  } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
    {
     if(fNumberOfRPsEBE>0){this->CalculateIntFlowCorrectionsForNUASinTermsUsingParticleWeights();}
     if(fNumberOfRPsEBE>0){this->CalculateIntFlowCorrectionsForNUACosTermsUsingParticleWeights();}     
    }      
  // Whether or not using particle weights the following is calculated in the same way:  
  if(fNumberOfRPsEBE>0){this->CalculateIntFlowProductOfCorrectionTermsForNUA();}     
  if(fNumberOfRPsEBE>0){this->CalculateIntFlowSumOfEventWeightsNUA();}     
  if(fNumberOfRPsEBE>0){this->CalculateIntFlowSumOfProductOfEventWeightsNUA();}     
  // Mixed harmonics:
  if(fCalculateMixedHarmonics){this->CalculateMixedHarmonics();}
 } // end of if(!fEvaluateIntFlowNestedLoops)

 // g) Call the methods which calculate correlations for differential flow:
 if(!fEvaluateDiffFlowNestedLoops && fCalculateDiffFlow)
 {
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   // Without using particle weights:
   this->CalculateDiffFlowCorrelations("RP","Pt"); 
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrelations("RP","Eta");}
   this->CalculateDiffFlowCorrelations("POI","Pt");
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrelations("POI","Eta");}
   // Non-isotropic terms:
   this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Pt");
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Eta");}
   this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Pt");
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Eta");}
   this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Pt");
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Eta");}
   this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Pt");
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Eta");}   
  } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
    {
     // With using particle weights:   
     this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Pt"); 
     if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Eta");} 
     this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Pt"); 
     if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Eta");} 
     // Non-isotropic terms:
     this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("RP","Pt");
     if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("RP","Eta");}
     this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("POI","Pt");
     if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("POI","Eta");}
     this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("RP","Pt");
     if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("RP","Eta");}
     this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("POI","Pt");
     if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("POI","Eta");}   
    }     
  // Whether or not using particle weights the following is calculated in the same way:  
  this->CalculateDiffFlowProductOfCorrelations("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowProductOfCorrelations("RP","Eta");}
  this->CalculateDiffFlowProductOfCorrelations("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowProductOfCorrelations("POI","Eta");}
  this->CalculateDiffFlowSumOfEventWeights("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowSumOfEventWeights("RP","Eta");}
  this->CalculateDiffFlowSumOfEventWeights("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowSumOfEventWeights("POI","Eta");}
  this->CalculateDiffFlowSumOfProductOfEventWeights("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowSumOfProductOfEventWeights("RP","Eta");}
  this->CalculateDiffFlowSumOfProductOfEventWeights("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowSumOfProductOfEventWeights("POI","Eta");}   
 } // end of if(!fEvaluateDiffFlowNestedLoops && fCalculateDiffFlow)

 // h) Call the methods which calculate correlations for 2D differential flow:
 if(!fEvaluateDiffFlowNestedLoops && fCalculate2DDiffFlow)
 {
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   // Without using particle weights:
   this->Calculate2DDiffFlowCorrelations("RP"); 
   this->Calculate2DDiffFlowCorrelations("POI");
   // Non-isotropic terms:
   // ... to be ctd ...
  } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
    {
     // With using particle weights:   
     // ... to be ctd ...  
     // Non-isotropic terms:
     // ... to be ctd ...
    }     
  // Whether or not using particle weights the following is calculated in the same way:  
  // ... to be ctd ...   
 } // end of if(!fEvaluateDiffFlowNestedLoops && fCalculate2DDiffFlow)
 
 // i) Call the methods which calculate other differential correlators:
 if(!fEvaluateDiffFlowNestedLoops && fCalculateDiffFlow)
 {
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   // Without using particle weights:
   this->CalculateOtherDiffCorrelators("RP","Pt"); 
   if(fCalculateDiffFlowVsEta){this->CalculateOtherDiffCorrelators("RP","Eta");}
   this->CalculateOtherDiffCorrelators("POI","Pt"); 
   if(fCalculateDiffFlowVsEta){this->CalculateOtherDiffCorrelators("POI","Eta");}     
  } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
    {
     // With using particle weights:   
     // ... to be ctd ...  
    }     
  // Whether or not using particle weights the following is calculated in the same way:  
  // ... to be ctd ...   
 } // end of if(!fEvaluateDiffFlowNestedLoops)
 
 // j) Distributions of correlations:
 if(fStoreDistributions){this->StoreDistributionsOfCorrelations();}
 
 // k) Store phi distribution for one event to illustrate flow: 
 if(fStorePhiDistributionForOneEvent){this->StorePhiDistributionForOneEvent(anEvent);}
   
 // l) Cross-check with nested loops correlators for reference flow:
 if(fEvaluateIntFlowNestedLoops){this->EvaluateIntFlowNestedLoops(anEvent);} 

 // m) Cross-check with nested loops correlators for differential flow:
 if(fEvaluateDiffFlowNestedLoops){this->EvaluateDiffFlowNestedLoops(anEvent);} 
 
 // n) Reset all event-by-event quantities (very important !!!!):
 this->ResetEventByEventQuantities();
 
} // end of AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::Finish()
{
 // Calculate the final results.

 // a) Check all pointers used in this method;
 // b) Access the constants;
 // c) Access the flags;
 // d) Calculate reference cumulants (not corrected for detector effects);
 // e) Correct reference cumulants for detector effects;
 // f) Calculate reference flow;
 // g) Store results for reference flow in AliFlowCommonHistResults and print them on the screen;
 // h) Calculate the final results for differential flow (without/with weights);
 // i) Correct the results for differential flow (without/with weights) for effects of non-uniform acceptance (NUA);
 // j) Calculate the final results for integrated flow (RP/POI) and store in AliFlowCommonHistResults;
 // k) Store results for differential flow in AliFlowCommonHistResults;
 // l) Print the final results for integrated flow (RP/POI) on the screen; 
 // m) Cross-checking: Results from Q-vectors vs results from nested loops;
 // i) Calculate cumulants for mixed harmonics;
 // j) Calculate cumulants for bootstrap.

 // a) Check all pointers used in this method:
 this->CheckPointersUsedInFinish();
  
 // b) Access the constants:
 this->CommonConstants("Finish");          
 
 if(fCommonHists && fCommonHists->GetHarmonic()) // to be improved (moved somewhere else)
 {
  fHarmonic = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1);
 } 
 
 // c) Access the flags: // tbi (implement a method for this? should I store again the flags because they can get modified with redoFinish?)
 fUsePhiWeights = (Bool_t)fUseParticleWeights->GetBinContent(1); 
 fUsePtWeights = (Bool_t)fUseParticleWeights->GetBinContent(2); 
 fUseEtaWeights = (Bool_t)fUseParticleWeights->GetBinContent(3);  
 fUseTrackWeights = (Bool_t)fUseParticleWeights->GetBinContent(4);  
 fApplyCorrectionForNUA = (Bool_t)fIntFlowFlags->GetBinContent(3); 
 fPrintFinalResults[0] = (Bool_t)fIntFlowFlags->GetBinContent(4);
 fPrintFinalResults[1] = (Bool_t)fIntFlowFlags->GetBinContent(5);
 fPrintFinalResults[2] = (Bool_t)fIntFlowFlags->GetBinContent(6);
 fPrintFinalResults[3] = (Bool_t)fIntFlowFlags->GetBinContent(7);
 fApplyCorrectionForNUAVsM = (Bool_t)fIntFlowFlags->GetBinContent(8);  
 fPropagateErrorAlsoFromNIT = (Bool_t)fIntFlowFlags->GetBinContent(9);  
 fCalculateCumulantsVsM = (Bool_t)fIntFlowFlags->GetBinContent(10); 
 fMinimumBiasReferenceFlow = (Bool_t)fIntFlowFlags->GetBinContent(11); 
 fForgetAboutCovariances = (Bool_t)fIntFlowFlags->GetBinContent(12);
 fStorePhiDistributionForOneEvent = (Bool_t)fIntFlowFlags->GetBinContent(13);
 fFillMultipleControlHistograms = (Bool_t)fIntFlowFlags->GetBinContent(14); 
 fCalculateAllCorrelationsVsM = (Bool_t)fIntFlowFlags->GetBinContent(15);
 fUse2DHistograms = (Bool_t)fIntFlowFlags->GetBinContent(18);
 fFillProfilesVsMUsingWeights = (Bool_t)fIntFlowFlags->GetBinContent(19);
 fUseQvectorTerms = (Bool_t)fIntFlowFlags->GetBinContent(20);
 fMaxCommonResultsHistogram = (Bool_t)fIntFlowFlags->GetBinContent(21);
 fEvaluateIntFlowNestedLoops = (Bool_t)fEvaluateNestedLoops->GetBinContent(1);
 fEvaluateDiffFlowNestedLoops = (Bool_t)fEvaluateNestedLoops->GetBinContent(2); 
 fCrossCheckInPtBinNo = (Int_t)fEvaluateNestedLoops->GetBinContent(3);
 fCrossCheckInEtaBinNo = (Int_t)fEvaluateNestedLoops->GetBinContent(4); 
 fCalculateMixedHarmonics = (Bool_t)fMixedHarmonicsFlags->GetBinContent(1);
 //fHarmonic = (Int_t)fMixedHarmonicsFlags->GetBinContent(2); // TBI should I add inpdependent generic harmonic here?
 fCalculateMixedHarmonicsVsM = (Bool_t)fMixedHarmonicsFlags->GetBinContent(3);
 fUseBootstrap = (Bool_t)fBootstrapFlags->GetBinContent(1); 
 fUseBootstrapVsM = (Bool_t)fBootstrapFlags->GetBinContent(2); 
 fnSubsamples = (Int_t)fBootstrapFlags->GetBinContent(3); 

 // d) Calculate reference cumulants (not corrected for detector effects):
 this->FinalizeCorrelationsIntFlow();
 this->CalculateCovariancesIntFlow();
 this->CalculateCumulantsIntFlow();

 // e) Correct reference cumulants for detector effects:
 this->FinalizeCorrectionTermsForNUAIntFlow();
 this->CalculateCovariancesNUAIntFlow(); 
 this->CalculateQcumulantsCorrectedForNUAIntFlow();  

 // f) Calculate reference flow:
 this->CalculateReferenceFlow(); 
  
 // g) Store results for reference flow in AliFlowCommonHistResults and print them on the screen:
 this->FillCommonHistResultsIntFlow();  
 if(fPrintFinalResults[0]){this->PrintFinalResultsForIntegratedFlow("RF");}
 if(fPrintFinalResults[3] && fCalculateCumulantsVsM){this->PrintFinalResultsForIntegratedFlow("RF, rebinned in M");}
 
 // h) Calculate the final results for differential flow (without/with weights):
 if(fCalculateDiffFlow)
 {
  this->FinalizeReducedCorrelations("RP","Pt"); 
  if(fCalculateDiffFlowVsEta){this->FinalizeReducedCorrelations("RP","Eta");} 
  this->FinalizeReducedCorrelations("POI","Pt"); 
  if(fCalculateDiffFlowVsEta){this->FinalizeReducedCorrelations("POI","Eta");}
  this->CalculateDiffFlowCovariances("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCovariances("RP","Eta");}
  this->CalculateDiffFlowCovariances("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCovariances("POI","Eta");}
  this->CalculateDiffFlowCumulants("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCumulants("RP","Eta");}
  this->CalculateDiffFlowCumulants("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCumulants("POI","Eta");}
  this->CalculateDiffFlow("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlow("RP","Eta");}
  this->CalculateDiffFlow("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlow("POI","Eta");}
 } // if(fCalculateDiffFlow)
 
 // i) Correct the results for differential flow (without/with weights) for effects of non-uniform acceptance (NUA):
 if(fCalculateDiffFlow)
 {
  this->FinalizeCorrectionTermsForNUADiffFlow("RP","Pt");
  if(fCalculateDiffFlowVsEta){this->FinalizeCorrectionTermsForNUADiffFlow("RP","Eta");}
  this->FinalizeCorrectionTermsForNUADiffFlow("POI","Pt");
  if(fCalculateDiffFlowVsEta){this->FinalizeCorrectionTermsForNUADiffFlow("POI","Eta");}      
  this->CalculateDiffFlowCumulantsCorrectedForNUA("RP","Pt");   
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCumulantsCorrectedForNUA("RP","Eta");}   
  this->CalculateDiffFlowCumulantsCorrectedForNUA("POI","Pt");   
  if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCumulantsCorrectedForNUA("POI","Eta");}  
  if(fApplyCorrectionForNUA)
  {
   this->CalculateDiffFlowCorrectedForNUA("RP","Pt"); 
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectedForNUA("RP","Eta");} 
   this->CalculateDiffFlowCorrectedForNUA("POI","Pt"); 
   if(fCalculateDiffFlowVsEta){this->CalculateDiffFlowCorrectedForNUA("POI","Eta");} 
  }
 } // end of if(fCalculateDiffFlow && fApplyCorrectionForNUA)
 
 // i) Calcualate final results for 2D differential flow: 
 if(fCalculate2DDiffFlow)
 {
  this->Calculate2DDiffFlowCumulants("RP");
  this->Calculate2DDiffFlowCumulants("POI");
  this->Calculate2DDiffFlow("RP");  
  this->Calculate2DDiffFlow("POI");  
 } // end of if(fCalculate2DDiffFlow)
    
 // j) Calculate the final results for integrated flow (RP/POI) and store in AliFlowCommonHistResults:
 if(fCalculateDiffFlow)
 {
  this->CalculateFinalResultsForRPandPOIIntegratedFlow("RP");
  this->CalculateFinalResultsForRPandPOIIntegratedFlow("POI");
 }
 
 // k) Store results for differential flow in AliFlowCommonHistResults:
 if(fCalculateDiffFlow)
 {
  this->FillCommonHistResultsDiffFlow("RP");
  this->FillCommonHistResultsDiffFlow("POI");
 }
 
 // l) Print the final results for integrated flow (RP/POI) on the screen:
 if(fPrintFinalResults[1] && fCalculateDiffFlow){this->PrintFinalResultsForIntegratedFlow("RP");} 
 if(fPrintFinalResults[2] && fCalculateDiffFlow){this->PrintFinalResultsForIntegratedFlow("POI");}
    
 // m) Cross-checking: Results from Q-vectors vs results from nested loops:
 //  m1) Reference flow:
 if(fEvaluateIntFlowNestedLoops)
 {
  this->CrossCheckIntFlowCorrelations();
  this->CrossCheckIntFlowCorrectionTermsForNUA(); 
  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights){this->CrossCheckIntFlowExtraCorrelations();}
  if(fCalculateMixedHarmonics){this->CrossCheckIntFlowCorrelations();}     
 } // end of if(fEvaluateIntFlowNestedLoops)  
 //  m2) Differential flow: 
 if(fEvaluateDiffFlowNestedLoops && fCalculateDiffFlow) 
 {
  // Correlations:
  this->PrintNumberOfParticlesInSelectedBin();
  this->CrossCheckDiffFlowCorrelations("RP","Pt");  
  if(fCalculateDiffFlowVsEta){this->CrossCheckDiffFlowCorrelations("RP","Eta");} 
  this->CrossCheckDiffFlowCorrelations("POI","Pt");  
  if(fCalculateDiffFlowVsEta){this->CrossCheckDiffFlowCorrelations("POI","Eta");}
  // Correction terms for non-uniform acceptance:
  this->CrossCheckDiffFlowCorrectionTermsForNUA("RP","Pt");      
  if(fCalculateDiffFlowVsEta){this->CrossCheckDiffFlowCorrectionTermsForNUA("RP","Eta");}       
  this->CrossCheckDiffFlowCorrectionTermsForNUA("POI","Pt");      
  if(fCalculateDiffFlowVsEta){this->CrossCheckDiffFlowCorrectionTermsForNUA("POI","Eta");}
  // Other differential correlators:       
  this->CrossCheckOtherDiffCorrelators("RP","Pt");  
  if(fCalculateDiffFlowVsEta){this->CrossCheckOtherDiffCorrelators("RP","Eta");} 
  this->CrossCheckOtherDiffCorrelators("POI","Pt");  
  if(fCalculateDiffFlowVsEta){this->CrossCheckOtherDiffCorrelators("POI","Eta");}
 } // end of if(fEvaluateDiffFlowNestedLoops)
                                                                                
 // i) Calculate cumulants for mixed harmonics: 
 if(fCalculateMixedHarmonics){this->CalculateCumulantsMixedHarmonics();} 

 // j) Calculate cumulants for bootstrap:
 if(fUseBootstrap||fUseBootstrapVsM){this->CalculateCumulantsForBootstrap();} 

} // end of AliFlowAnalysisWithQCumulants::Finish()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::EvaluateIntFlowNestedLoops(AliFlowEventSimple* anEvent)
{
 // Evaluate all correlators for reference flow with nested loops.
 
 Int_t nPrim = anEvent->NumberOfTracks(); // number of primaries
 if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10 
 {
  // Without using particle weights:
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   // Correlations:
   this->CalculateIntFlowCorrelations(); // from Q-vectors
   this->EvaluateIntFlowCorrelationsWithNestedLoops(anEvent); // from nested loops (TBI: do I have to pass here anEvent or not?)
   // Correction for non-uniform acceptance:
   this->CalculateIntFlowCorrectionsForNUASinTerms(); // from Q-vectors (sin terms)
   this->CalculateIntFlowCorrectionsForNUACosTerms(); // from Q-vectors (cos terms)
   this->EvaluateIntFlowCorrectionsForNUAWithNestedLoops(anEvent); // from nested loops (both sin and cos terms)
   // Mixed harmonics:
   if(fCalculateMixedHarmonics)
   {
    this->CalculateMixedHarmonics(); // from Q-vectors  
    this->EvaluateMixedHarmonicsWithNestedLoops(anEvent); // from nested loops (TBI: do I have to pass here anEvent or not?)
   } // end of if(fCalculateMixedHarmonics)
  }
  // Using particle weights:
  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights)
  {
   // Correlations
   this->CalculateIntFlowCorrelationsUsingParticleWeights(); // from Q-vectors
   this->EvaluateIntFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent); // from nested loops (to be improved: do I have to pass here anEvent or not?)
   // Correction for non-uniform acceptance:
   this->CalculateIntFlowCorrectionsForNUASinTermsUsingParticleWeights(); // from Q-vectors (sin terms)
   this->CalculateIntFlowCorrectionsForNUACosTermsUsingParticleWeights(); // from Q-vectors (cos terms)
   this->EvaluateIntFlowCorrectionsForNUAWithNestedLoopsUsingParticleWeights(anEvent); // from nested loops (both sin and cos terms)   
  }
 } else if(nPrim>fMaxAllowedMultiplicity) // to if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity)
   {
    cout<<endl;
    cout<<"Skipping the event because multiplicity is "<<nPrim<<". Too high to evaluate nested loops!"<<endl;
   } else
     {
      cout<<endl;
      cout<<"Skipping the event because multiplicity is "<<nPrim<<"."<<endl;      
     } 

} // end of void AliFlowAnalysisWithQCumulants::EvaluateIntFlowNestedLoops(AliFlowEventSimple* anEvent)

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::EvaluateDiffFlowNestedLoops(AliFlowEventSimple* anEvent)
{
 // Evalauted all correlators for differential flow with nested loops.

 if(!fCalculateDiffFlow){return;}

 Int_t nPrim = anEvent->NumberOfTracks(); // number of primaries
 if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10
 {
  // Without using particle weights:
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   // 1.) Reduced correlations:
   //  Q-vectors:
   this->CalculateDiffFlowCorrelations("RP","Pt");
   this->CalculateDiffFlowCorrelations("RP","Eta");
   this->CalculateDiffFlowCorrelations("POI","Pt");
   this->CalculateDiffFlowCorrelations("POI","Eta");
   //  Nested loops:
   this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"RP","Pt"); 
   this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"RP","Eta"); 
   this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"POI","Pt"); 
   this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"POI","Eta"); 
   // 2.) Reduced corrections for non-uniform acceptance:
   //  Q-vectors:
   this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Pt");
   this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Eta");
   this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Pt");
   this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Eta");
   this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Pt");
   this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Eta");
   this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Pt");
   this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Eta");
   //  Nested loops:
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"RP","Pt");
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"RP","Eta");
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"POI","Pt"); 
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"POI","Eta"); 
   // 3.) Other differential correlators:
   //  Q-vectors:
   this->CalculateOtherDiffCorrelators("RP","Pt");
   this->CalculateOtherDiffCorrelators("RP","Eta");
   this->CalculateOtherDiffCorrelators("POI","Pt");
   this->CalculateOtherDiffCorrelators("POI","Eta");   
   //  Nested loops:
   this->EvaluateOtherDiffCorrelatorsWithNestedLoops(anEvent,"RP","Pt");
   this->EvaluateOtherDiffCorrelatorsWithNestedLoops(anEvent,"RP","Eta");
   this->EvaluateOtherDiffCorrelatorsWithNestedLoops(anEvent,"POI","Pt");
   this->EvaluateOtherDiffCorrelatorsWithNestedLoops(anEvent,"POI","Eta");   
  } // end of if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  // Using particle weights:
  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights)
  {
   this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Pt"); 
   this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Eta"); 
   this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Pt"); 
   this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Eta"); 
   this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("RP","Pt");
   this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("RP","Eta");
   this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("POI","Pt");
   this->CalculateDiffFlowCorrectionsForNUASinTermsUsingParticleWeights("POI","Eta");
   this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("RP","Pt");
   this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("RP","Eta");
   this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("POI","Pt");
   this->CalculateDiffFlowCorrectionsForNUACosTermsUsingParticleWeights("POI","Eta");
   this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"RP","Pt"); 
   this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"RP","Eta");
   this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"POI","Pt"); 
   this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"POI","Eta");   
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoopsUsingParticleWeights(anEvent,"RP","Pt"); 
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoopsUsingParticleWeights(anEvent,"RP","Eta"); 
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoopsUsingParticleWeights(anEvent,"POI","Pt"); 
   this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoopsUsingParticleWeights(anEvent,"POI","Eta"); 
  } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights)
 } // end of if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10

} // end of void AliFlowAnalysisWithQCumulants::EvaluateDiffFlowNestedLoops(AliFlowEventSimple* anEvent)

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUACosTerms()
{
 // Calculate correction terms for non-uniform acceptance of the detector for reference flow (cos terms).
 
 // multiplicity:
 Double_t dMult = (*fSpk)(0,0);
 
 // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
 Double_t dReQ1n = (*fReQ)(0,0);
 Double_t dReQ2n = (*fReQ)(1,0);
 //Double_t dReQ3n = (*fReQ)(2,0);
 //Double_t dReQ4n = (*fReQ)(3,0);
 Double_t dImQ1n = (*fImQ)(0,0);
 Double_t dImQ2n = (*fImQ)(1,0);
 //Double_t dImQ3n = (*fImQ)(2,0);
 //Double_t dImQ4n = (*fImQ)(3,0);

 // Multiplicity bin of an event (relevant for all histos vs M): 
 Double_t dMultiplicityBin = 0.;
 if(fMultiplicityIs==AliFlowCommonConstants::kRP)
 {
  dMultiplicityBin = fNumberOfRPsEBE+0.5;
 } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
   {
    dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
     {
      dMultiplicityBin = fNumberOfPOIsEBE+0.5;
     } 

 //                                  *************************************************************
 //                                  **** corrections for non-uniform acceptance (cos terms): ****
 //                                  *************************************************************
 //
 // Remark 1: corrections for non-uniform acceptance (cos terms) calculated with non-weighted Q-vectors 
 //           are stored in 1D profile fQCorrectionsCos.
 // Remark 2: binning of fIntFlowCorrectionTermsForNUAPro[1] is organized as follows:
 // --------------------------------------------------------------------------------------------------------------------
 // 1st bin: <<cos(n*(phi1))>> = cosP1n
 // 2nd bin: <<cos(n*(phi1+phi2))>> = cosP1nP1n
 // 3rd bin: <<cos(n*(phi1-phi2-phi3))>> = cosP1nM1nM1n
 // 4th bin: <<cos(n*(2phi1-phi2))>> = cosP2nM1n
 // --------------------------------------------------------------------------------------------------------------------
  
 // 1-particle:
 Double_t cosP1n = 0.; // <<cos(n*(phi1))>>
   
 if(dMult>0)
 {
  cosP1n = dReQ1n/dMult; 
  
  // average non-weighted 1-particle correction (cos terms) for non-uniform acceptance for single event:
  fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(1,cosP1n);
  // event weights for NUA terms:
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[1]->SetBinContent(1,dMult);
  
  // final average non-weighted 1-particle correction (cos terms) for non-uniform acceptance for all events:
  fIntFlowCorrectionTermsForNUAPro[1]->Fill(0.5,cosP1n,dMult);  
  if(fCalculateCumulantsVsM){fIntFlowCorrectionTermsForNUAVsMPro[1][0]->Fill(dMultiplicityBin,cosP1n,dMult);}    
 } 
 
 // 2-particle:
 Double_t cosP1nP1n = 0.; // <<cos(n*(phi1+phi2))>>
 Double_t cosP2nM1n = 0.; // <<cos(n*(2phi1-phi2))>>
 
 if(dMult>1)
 {
  cosP1nP1n = (pow(dReQ1n,2)-pow(dImQ1n,2)-dReQ2n)/(dMult*(dMult-1)); 
  cosP2nM1n = (dReQ2n*dReQ1n+dImQ2n*dImQ1n-dReQ1n)/(dMult*(dMult-1)); 
  
  // average non-weighted 2-particle correction (cos terms) for non-uniform acceptance for single event:
  fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(2,cosP1nP1n);
  fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(4,cosP2nM1n);
  // event weights for NUA terms:
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[1]->SetBinContent(2,dMult*(dMult-1));
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[1]->SetBinContent(4,dMult*(dMult-1));
      
  // final average non-weighted 2-particle correction (cos terms) for non-uniform acceptance for all events:
  fIntFlowCorrectionTermsForNUAPro[1]->Fill(1.5,cosP1nP1n,dMult*(dMult-1));  
  fIntFlowCorrectionTermsForNUAPro[1]->Fill(3.5,cosP2nM1n,dMult*(dMult-1));
  if(fCalculateCumulantsVsM)
  {
   fIntFlowCorrectionTermsForNUAVsMPro[1][1]->Fill(dMultiplicityBin,cosP1nP1n,dMult*(dMult-1));  
   fIntFlowCorrectionTermsForNUAVsMPro[1][3]->Fill(dMultiplicityBin,cosP2nM1n,dMult*(dMult-1));
  }
 } 
 
 // 3-particle:
 Double_t cosP1nM1nM1n = 0.; // <<cos(n*(phi1-phi2-phi3))>>
 
 if(dMult>2)
 {
  cosP1nM1nM1n = (dReQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))-dReQ1n*dReQ2n-dImQ1n*dImQ2n-2.*(dMult-1)*dReQ1n)
               / (dMult*(dMult-1)*(dMult-2)); 
  
  // average non-weighted 3-particle correction (cos terms) for non-uniform acceptance for single event:
  fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(3,cosP1nM1nM1n);
  // event weights for NUA terms:
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[1]->SetBinContent(3,dMult*(dMult-1)*(dMult-2));
  
  // final average non-weighted 3-particle correction (cos terms) for non-uniform acceptance for all events:
  fIntFlowCorrectionTermsForNUAPro[1]->Fill(2.5,cosP1nM1nM1n,dMult*(dMult-1)*(dMult-2));
  if(fCalculateCumulantsVsM){fIntFlowCorrectionTermsForNUAVsMPro[1][2]->Fill(dMultiplicityBin,cosP1nM1nM1n,dMult*(dMult-1)*(dMult-2));}  
 } 
 
} // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUACosTerms()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUASinTerms()
{
 // calculate corrections for non-uniform acceptance of the detector for no-name integrated flow (sin terms)
 
 // multiplicity:
 Double_t dMult = (*fSpk)(0,0);
 
 // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
 Double_t dReQ1n = (*fReQ)(0,0);
 Double_t dReQ2n = (*fReQ)(1,0);
 //Double_t dReQ3n = (*fReQ)(2,0);
 //Double_t dReQ4n = (*fReQ)(3,0);
 Double_t dImQ1n = (*fImQ)(0,0);
 Double_t dImQ2n = (*fImQ)(1,0);
 //Double_t dImQ3n = (*fImQ)(2,0);
 //Double_t dImQ4n = (*fImQ)(3,0);
        
 // Multiplicity bin of an event (relevant for all histos vs M): 
 Double_t dMultiplicityBin = 0.;
 if(fMultiplicityIs==AliFlowCommonConstants::kRP)
 {
  dMultiplicityBin = fNumberOfRPsEBE+0.5;
 } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
   {
    dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
     {
      dMultiplicityBin = fNumberOfPOIsEBE+0.5;
     } 

 //                                  *************************************************************
 //                                  **** corrections for non-uniform acceptance (sin terms): ****
 //                                  *************************************************************
 //
 // Remark 1: corrections for non-uniform acceptance (sin terms) calculated with non-weighted Q-vectors 
 //           are stored in 1D profile fQCorrectionsSin.
 // Remark 2: binning of fIntFlowCorrectionTermsForNUAPro[0] is organized as follows:
 // --------------------------------------------------------------------------------------------------------------------
 // 1st bin: <<sin(n*(phi1))>> = sinP1n
 // 2nd bin: <<sin(n*(phi1+phi2))>> = sinP1nP1n
 // 3rd bin: <<sin(n*(phi1-phi2-phi3))>> = sinP1nM1nM1n
 // 4th bin: <<sin(n*(2phi1-phi2))>> = sinP2nM1n
 // --------------------------------------------------------------------------------------------------------------------
 
 // 1-particle:
 Double_t sinP1n = 0.; // <sin(n*(phi1))>
 
 if(dMult>0)
 {
  sinP1n = dImQ1n/dMult; 
     
  // average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for single event:
  fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(1,sinP1n);  
  // event weights for NUA terms:
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[0]->SetBinContent(1,dMult);
  
  // final average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for all events:   
  fIntFlowCorrectionTermsForNUAPro[0]->Fill(0.5,sinP1n,dMult);  
  if(fCalculateCumulantsVsM){fIntFlowCorrectionTermsForNUAVsMPro[0][0]->Fill(dMultiplicityBin,sinP1n,dMult);} 
 } 
 
 // 2-particle:
 Double_t sinP1nP1n = 0.; // <<sin(n*(phi1+phi2))>>
 Double_t sinP2nM1n = 0.; // <<sin(n*(2phi1-phi2))>>
 if(dMult>1)
 {
  sinP1nP1n = (2.*dReQ1n*dImQ1n-dImQ2n)/(dMult*(dMult-1)); 
  sinP2nM1n = (dImQ2n*dReQ1n-dReQ2n*dImQ1n-dImQ1n)/(dMult*(dMult-1)); 
     
  // average non-weighted 2-particle correction (sin terms) for non-uniform acceptance for single event:
  fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(2,sinP1nP1n);
  fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(4,sinP2nM1n);
  // event weights for NUA terms:
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[0]->SetBinContent(2,dMult*(dMult-1));
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[0]->SetBinContent(4,dMult*(dMult-1));
  
  // final average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for all events:      
  fIntFlowCorrectionTermsForNUAPro[0]->Fill(1.5,sinP1nP1n,dMult*(dMult-1));  
  fIntFlowCorrectionTermsForNUAPro[0]->Fill(3.5,sinP2nM1n,dMult*(dMult-1));  
  if(fCalculateCumulantsVsM)
  {
   fIntFlowCorrectionTermsForNUAVsMPro[0][1]->Fill(dMultiplicityBin,sinP1nP1n,dMult*(dMult-1));  
   fIntFlowCorrectionTermsForNUAVsMPro[0][3]->Fill(dMultiplicityBin,sinP2nM1n,dMult*(dMult-1));    
  }
 } 
 
 // 3-particle:
 Double_t sinP1nM1nM1n = 0.; // <<sin(n*(phi1-phi2-phi3))>>
 
 if(dMult>2)
 {
  sinP1nM1nM1n = (-dImQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))+dReQ1n*dImQ2n-dImQ1n*dReQ2n+2.*(dMult-1)*dImQ1n)
               / (dMult*(dMult-1)*(dMult-2)); 
  
  // average non-weighted 3-particle correction (sin terms) for non-uniform acceptance for single event:
  fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(3,sinP1nM1nM1n);
  // event weights for NUA terms:
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[0]->SetBinContent(3,dMult*(dMult-1)*(dMult-2));
  
  // final average non-weighted 3-particle correction (sin terms) for non-uniform acceptance for all events:  
  fIntFlowCorrectionTermsForNUAPro[0]->Fill(2.5,sinP1nM1nM1n,dMult*(dMult-1)*(dMult-2));
  if(fCalculateCumulantsVsM){fIntFlowCorrectionTermsForNUAVsMPro[0][2]->Fill(dMultiplicityBin,sinP1nM1nM1n,dMult*(dMult-1)*(dMult-2));}  
 } 
 
} // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUASinTerms()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)
{
 // a) Get pointers for common control and common result histograms;
 // b) Get pointers for histograms holding particle weights;
 // c) Get pointers for reference flow histograms;
 // d) Get pointers for differential flow histograms;
 // e) Get pointers for 2D differential flow histograms;
 // f) Get pointers for other differential correlators;
 // g) Get pointers for mixed harmonics histograms;
 // h) Get pointers for nested loops' histograms;
 // i) Get pointers for control histograms;
 // j) Get pointers for bootstrap.
 
 if(outputListHistos)
 {  
  this->SetHistList(outputListHistos);
  if(!fHistList)
  {
   printf("\n WARNING (QC): fHistList is NULL in AFAWQC::GOH() !!!!\n\n");
   exit(0);
  }
  this->GetPointersForCommonHistograms(); 
  this->GetPointersForParticleWeightsHistograms(); 
  this->GetPointersForIntFlowHistograms();
  this->GetPointersForDiffFlowHistograms(); 
  this->GetPointersFor2DDiffFlowHistograms(); 
  this->GetPointersForOtherDiffCorrelators();  
  this->GetPointersForMixedHarmonicsHistograms(); 
  this->GetPointersForNestedLoopsHistograms(); 
  this->GetPointersForControlHistograms();
  this->GetPointersForBootstrap();
 } else 
   {
    printf("\n WARNING (QC): outputListHistos is NULL in AFAWQC::GOH() !!!!\n\n");
    exit(0);
   }
   
} // end of void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)

//=======================================================================================================================

TProfile* AliFlowAnalysisWithQCumulants::MakePtProjection(TProfile2D *profilePtEta) const
{
 // project 2D profile onto pt axis to get 1D profile
 
 Int_t nBinsPt   = profilePtEta->GetNbinsX();
 Double_t dPtMin = (profilePtEta->GetXaxis())->GetXmin();
 Double_t dPtMax = (profilePtEta->GetXaxis())->GetXmax();
 
 Int_t nBinsEta   = profilePtEta->GetNbinsY();
 
 TProfile *profilePt = new TProfile("","",nBinsPt,dPtMin,dPtMax); 
 
 for(Int_t p=1;p<=nBinsPt;p++)
 {
  Double_t contentPt = 0.;
  Double_t entryPt = 0.;
  Double_t spreadPt = 0.;
  Double_t sum1 = 0.;
  Double_t sum2 = 0.;
  Double_t sum3 = 0.;
  for(Int_t e=1;e<=nBinsEta;e++)
  {
   contentPt += (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)))
              * (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
   entryPt   += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
   
   sum1 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)))
         * (pow(profilePtEta->GetBinError(profilePtEta->GetBin(p,e)),2.)
            + pow(profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)),2.)); 
   sum2 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
   sum3 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)))
         * (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)));            
  }
  if(sum2>0. && sum1/sum2-pow(sum3/sum2,2.) > 0.)
  {
   spreadPt = pow(sum1/sum2-pow(sum3/sum2,2.),0.5);
  }
  profilePt->SetBinContent(p,contentPt);
  profilePt->SetBinEntries(p,entryPt);
  {
   profilePt->SetBinError(p,spreadPt);
  }
  
 }
 
 return profilePt;
 
} // end of TProfile* AliFlowAnalysisWithQCumulants::MakePtProjection(TProfile2D *profilePtEta)


//=======================================================================================================================


TProfile* AliFlowAnalysisWithQCumulants::MakeEtaProjection(TProfile2D *profilePtEta) const
{
 // project 2D profile onto eta axis to get 1D profile
 
 Int_t nBinsEta   = profilePtEta->GetNbinsY();
 Double_t dEtaMin = (profilePtEta->GetYaxis())->GetXmin();
 Double_t dEtaMax = (profilePtEta->GetYaxis())->GetXmax();
 
 Int_t nBinsPt = profilePtEta->GetNbinsX();
 
 TProfile *profileEta = new TProfile("","",nBinsEta,dEtaMin,dEtaMax); 
 
 for(Int_t e=1;e<=nBinsEta;e++)
 {
  Double_t contentEta = 0.;
  Double_t entryEta = 0.;
  for(Int_t p=1;p<=nBinsPt;p++)
  {
   contentEta += (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)))
              * (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
   entryEta   += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
  }
  profileEta->SetBinContent(e,contentEta);
  profileEta->SetBinEntries(e,entryEta);
 }
 
 return profileEta;
 
} // end of TProfile* AliFlowAnalysisWithQCumulants::MakeEtaProjection(TProfile2D *profilePtEta)

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::PrintFinalResultsForIntegratedFlow(TString type)
{
 // Printing on the screen the final results for integrated flow (RF, POI and RP). 
 
 Int_t n = fHarmonic; 
 
 Double_t dVn[4] = {0.}; // array to hold Vn{2}, Vn{4}, Vn{6} and Vn{8}   
 Double_t dVnErr[4] = {0.}; // array to hold errors of Vn{2}, Vn{4}, Vn{6} and Vn{8}   
 
 if(type == "RF")
 {
  for(Int_t b=0;b<4;b++)
  {
   if(fCommonHistsResults2nd)
   {
    dVn[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinContent(1);
    dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinError(1);
   }
   if(fCommonHistsResults4th)
   {
    dVn[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinContent(1);
    dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinError(1);
   }
   if(fCommonHistsResults6th)
   {
    dVn[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinContent(1);
    dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinError(1);
   }
   if(fCommonHistsResults8th)
   {
    dVn[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinContent(1);
    dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinError(1);
   }
  } // for(Int_t b=0;b<4;b++)
 } else if(type == "RP")
   {
    if(fCommonHistsResults2nd)
    {
     dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinContent(1);
     dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinError(1);
    }
    if(fCommonHistsResults4th)
    {
     dVn[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinContent(1);
     dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinError(1);
    }
    if(fCommonHistsResults6th)
    {
     dVn[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinContent(1);
     dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinError(1);
    }
    if(fCommonHistsResults8th)
    {
     dVn[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinContent(1);
     dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinError(1);
    }
   } else if(type == "POI")
     {
      if(fCommonHistsResults2nd)
      {
       dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinContent(1);
       dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinError(1);
      }
      if(fCommonHistsResults4th)
      {
       dVn[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinContent(1);
       dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinError(1);
      }
      if(fCommonHistsResults6th)
      {
       dVn[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinContent(1);
       dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinError(1);
      }
      if(fCommonHistsResults8th)
      {
       dVn[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinContent(1);
       dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinError(1);
      }
     } else if(type == "RF, rebinned in M" && fCalculateCumulantsVsM)
       {
        for(Int_t b=0;b<4;b++)
        {
         dVn[b] = fIntFlowRebinnedInM->GetBinContent(b+1); 
         dVnErr[b] = fIntFlowRebinnedInM->GetBinError(b+1);
        }  
       }
 
 TString title = " flow estimates from Q-cumulants"; 
 TString subtitle = "    ("; 
 TString subtitle2 = "       (rebinned in M)"; 
 
 if(type != "RF, rebinned in M")
 {
  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
  {
   subtitle.Append(type);
   subtitle.Append(", without weights)");
  } else  
    {
     subtitle.Append(type);
     subtitle.Append(", with weights)");
    }
 } else
   {
    if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
    {
     subtitle.Append("RF");
     subtitle.Append(", without weights)");
    } else  
      {
       subtitle.Append("RF");
       subtitle.Append(", with weights)");      
      }
   } 
   
 cout<<endl;
 cout<<"*************************************"<<endl;
 cout<<"*************************************"<<endl;
 cout<<title.Data()<<endl; 
 cout<<subtitle.Data()<<endl; 
 if(type == "RF, rebinned in M"){cout<<subtitle2.Data()<<endl;}
 cout<<endl;
  
 for(Int_t i=0;i<4;i++)
 {
  cout<<"  v_"<<n<<"{"<<2*(i+1)<<"} = "<<dVn[i]<<" +/- "<<dVnErr[i]<<endl;
 }
 
 cout<<endl;
 if(type == "RF")
 {
  if(fApplyCorrectionForNUA)
  {
   cout<<" detector bias (corrected for): "<<endl;
  } else
    {
     cout<<" detector bias (not corrected for):"<<endl;  
    }
  cout<<"  to QC{2}: "<<fIntFlowDetectorBias->GetBinContent(1)<<" +/- "<<fIntFlowDetectorBias->GetBinError(1)<<endl;
  cout<<"  to QC{4}: "<<fIntFlowDetectorBias->GetBinContent(2)<<" +/- "<<fIntFlowDetectorBias->GetBinError(2)<<endl;
  cout<<endl;
 }
 if(type == "RF" || type == "RF, rebinned in M")
 {
  cout<<"     nEvts = "<<(Int_t)fCommonHists->GetHistMultRP()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()<<endl; 
 }
 else if (type == "RP")
 {
  cout<<"     nEvts = "<<(Int_t)fCommonHists->GetHistMultRP()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()<<endl;  
 } 
 else if (type == "POI")
 {
  cout<<"     nEvts = "<<(Int_t)fCommonHists->GetHistMultPOI()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultPOI()->GetMean()<<endl;
 }  
 
 cout<<"*************************************"<<endl;
 cout<<"*************************************"<<endl;
 cout<<endl; 
  
}// end of AliFlowAnalysisWithQCumulants::PrintFinalResultsForIntegratedFlow(TString type="RF");

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::WriteHistograms(TString outputFileName)
{
 //store the final results in output .root file
 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 //output->WriteObject(fHistList, "cobjQC","SingleKey");
 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 delete output;
}


//=======================================================================================================================


void AliFlowAnalysisWithQCumulants::WriteHistograms(TDirectoryFile *outputFileName)
{
 //store the final results in output .root file
 fHistList->SetName("cobjQC");
 fHistList->SetOwner(kTRUE);
 outputFileName->Add(fHistList);
 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
}

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::BookCommonHistograms()
{
 // Book common control histograms and common histograms for final results.
 //  a) Book common control histograms;
 //  b) Book common result histograms.
 
 // a) Book common control histograms: 
 //  Common control histograms (all events):
 TString commonHistsName = "AliFlowCommonHistQC";
 commonHistsName += fAnalysisLabel->Data();
 fCommonHists = new AliFlowCommonHist(commonHistsName.Data(),commonHistsName.Data(),fBookOnlyBasicCCH);
 fHistList->Add(fCommonHists);  
 //  Common control histograms (selected events):
 if(fFillMultipleControlHistograms)
 {
  // Common control histogram filled for events with 2 and more reference particles:
  TString commonHists2ndOrderName = "AliFlowCommonHist2ndOrderQC";
  commonHists2ndOrderName += fAnalysisLabel->Data();
  fCommonHists2nd = new AliFlowCommonHist(commonHists2ndOrderName.Data(),commonHists2ndOrderName.Data(),fBookOnlyBasicCCH);
  fHistList->Add(fCommonHists2nd);  
  // Common control histogram filled for events with 2 and more reference particles:
  TString commonHists4thOrderName = "AliFlowCommonHist4thOrderQC";
  commonHists4thOrderName += fAnalysisLabel->Data();
  fCommonHists4th = new AliFlowCommonHist(commonHists4thOrderName.Data(),commonHists4thOrderName.Data(),fBookOnlyBasicCCH);
  fHistList->Add(fCommonHists4th);  
  // Common control histogram filled for events with 6 and more reference particles:
  TString commonHists6thOrderName = "AliFlowCommonHist6thOrderQC";
  commonHists6thOrderName += fAnalysisLabel->Data();
  fCommonHists6th = new AliFlowCommonHist(commonHists6thOrderName.Data(),commonHists6thOrderName.Data(),fBookOnlyBasicCCH);
  fHistList->Add(fCommonHists6th);  
  // Common control histogram filled for events with 8 and more reference particles:
  TString commonHists8thOrderName = "AliFlowCommonHist8thOrderQC";
  commonHists8thOrderName += fAnalysisLabel->Data();
  fCommonHists8th = new AliFlowCommonHist(commonHists8thOrderName.Data(),commonHists8thOrderName.Data(),fBookOnlyBasicCCH);
  fHistList->Add(fCommonHists8th);    
 } // end of if(fFillMultipleControlHistograms)
 
 // b) Book common result histograms: 
 if(fMaxCommonResultsHistogram >= 2)
 {
  //  Common result histograms for QC{2}:
  TString commonHistResults2ndOrderName = "AliFlowCommonHistResults2ndOrderQC";
  commonHistResults2ndOrderName += fAnalysisLabel->Data();
  fCommonHistsResults2nd = new AliFlowCommonHistResults(commonHistResults2ndOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults2nd);
 }
 if(fMaxCommonResultsHistogram >= 4)
 {
  //  Common result histograms for QC{4}:
  TString commonHistResults4thOrderName = "AliFlowCommonHistResults4thOrderQC";
  commonHistResults4thOrderName += fAnalysisLabel->Data();
  fCommonHistsResults4th = new AliFlowCommonHistResults(commonHistResults4thOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults4th);
 }
 if(fMaxCommonResultsHistogram >= 6)
 {
  //  Common result histograms for QC{6}:
  TString commonHistResults6thOrderName = "AliFlowCommonHistResults6thOrderQC";
  commonHistResults6thOrderName += fAnalysisLabel->Data();
  fCommonHistsResults6th = new AliFlowCommonHistResults(commonHistResults6thOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults6th);
 }
 if(fMaxCommonResultsHistogram >= 8)
 {
  //  Common result histograms for QC{8}:
  TString commonHistResults8thOrderName = "AliFlowCommonHistResults8thOrderQC";
  commonHistResults8thOrderName += fAnalysisLabel->Data();
  fCommonHistsResults8th = new AliFlowCommonHistResults(commonHistResults8thOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults8th);
 }

} // end of void AliFlowAnalysisWithQCumulants::BookCommonHistograms()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::BookAndFillWeightsHistograms()
{
 // Book and fill histograms which hold phi, pt and eta weights.

 if(!fWeightsList)
 {
  printf("\n WARNING (QC): fWeightsList is NULL in AFAWQC::BAFWH() !!!! \n\n");
  exit(0);  
 }
    
 TString fUseParticleWeightsName = "fUseParticleWeightsQC";
 fUseParticleWeightsName += fAnalysisLabel->Data();
 fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",4,0,4);
 fUseParticleWeights->SetLabelSize(0.06);
 fUseParticleWeights->SetStats(kFALSE);
 (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
 (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
 (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
 (fUseParticleWeights->GetXaxis())->SetBinLabel(4,"w_{track}");
 fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
 fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
 fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
 fUseParticleWeights->Fill(3.5,(Int_t)fUseTrackWeights);
 fWeightsList->Add(fUseParticleWeights); 
  
 if(fUsePhiWeights)
 {
  if(fWeightsList->FindObject("phi_weights"))
  {
   fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
   if(!fPhiWeights)
   {
    printf("\n WARNING (QC): fPhiWeights is NULL in AFAWQC::BAFWH() !!!!\n\n");
    exit(0);
   }
   if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
   {
    cout<<endl;
    cout<<"WARNING (QC): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
    cout<<endl;
    //exit(0);
   }
  } else 
    {
     cout<<"WARNING: fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
     exit(0);
    }
 } // end of if(fUsePhiWeights)
 
 if(fUsePtWeights) 
 {
  if(fWeightsList->FindObject("pt_weights"))
  {
   fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
   if(!fPtWeights)
   {
    printf("\n WARNING (QC): fPtWeights is NULL in AFAWQC::BAFWH() !!!!\n\n");
    exit(0);
   }
   if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
   {
    cout<<endl;
    cout<<"WARNING (QC): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
    cout<<endl;
    //exit(0);
   }
  } else 
    {
     cout<<"WARNING: fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
     exit(0);
    }
 } // end of if(fUsePtWeights)    

 if(fUseEtaWeights) 
 {
  if(fWeightsList->FindObject("eta_weights"))
  {
   fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
   if(!fEtaWeights)
   {
    printf("\n WARNING (QC): fEtaWeights is NULL in AFAWQC::BAFWH() !!!!\n\n");
    exit(0);
   }
   if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
   {
    cout<<endl;
    cout<<"WARNING (QC): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
    cout<<endl;
    //exit(0);
   }
  } else 
    {
     cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
     exit(0);
    }
 } // end of if(fUseEtaWeights)
 
} // end of AliFlowAnalysisWithQCumulants::BookAndFillWeightsHistograms()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
{
 // Book all objects for integrated flow:
 //  a) Book profile to hold all flags for integrated flow;
 //  b) Book event-by-event quantities;
 //  c) Book profiles; // to be improved (comment)
 //  d) Book histograms holding the final results.
 
 TString sinCosFlag[2] = {"sin","cos"}; // to be improved (should I promote this to data members?)
 TString powerFlag[2] = {"linear","quadratic"}; // to be improved (should I promote this to data members?)
 
 // a) Book profile to hold all flags for integrated flow:
 TString intFlowFlagsName = "fIntFlowFlags";
 intFlowFlagsName += fAnalysisLabel->Data();
 fIntFlowFlags = new TProfile(intFlowFlagsName.Data(),"Flags for Integrated Flow",21,0.,21.);
 fIntFlowFlags->SetTickLength(-0.01,"Y");
 fIntFlowFlags->SetMarkerStyle(25);
 fIntFlowFlags->SetLabelSize(0.04);
 fIntFlowFlags->SetLabelOffset(0.02,"Y");
 fIntFlowFlags->SetStats(kFALSE);
 fIntFlowFlags->GetXaxis()->SetBinLabel(1,"Particle Weights");
 fIntFlowFlags->GetXaxis()->SetBinLabel(2,"Event Weights");
 fIntFlowFlags->GetXaxis()->SetBinLabel(3,"Corrected for NUA?");
 fIntFlowFlags->GetXaxis()->SetBinLabel(4,"Print RF results");
 fIntFlowFlags->GetXaxis()->SetBinLabel(5,"Print RP results");
 fIntFlowFlags->GetXaxis()->SetBinLabel(6,"Print POI results");
 fIntFlowFlags->GetXaxis()->SetBinLabel(7,"Print RF (rebinned in M) results");
 fIntFlowFlags->GetXaxis()->SetBinLabel(8,"Corrected for NUA vs M?");
 fIntFlowFlags->GetXaxis()->SetBinLabel(9,"Propagate errors to v_{n} from correlations?");
 fIntFlowFlags->GetXaxis()->SetBinLabel(10,"Calculate cumulants vs M");
 fIntFlowFlags->GetXaxis()->SetBinLabel(11,"fMinimumBiasReferenceFlow");
 fIntFlowFlags->GetXaxis()->SetBinLabel(12,"fForgetAboutCovariances");
 fIntFlowFlags->GetXaxis()->SetBinLabel(13,"fStorePhiDistributionForOneEvent");
 fIntFlowFlags->GetXaxis()->SetBinLabel(14,"fFillMultipleControlHistograms");
 fIntFlowFlags->GetXaxis()->SetBinLabel(15,"Calculate all correlations vs M");
 fIntFlowFlags->GetXaxis()->SetBinLabel(16,"fMultiplicityIs");
 fIntFlowFlags->GetXaxis()->SetBinLabel(17,"fExactNoRPs");
 fIntFlowFlags->GetXaxis()->SetBinLabel(18,"fUse2DHistograms");
 fIntFlowFlags->GetXaxis()->SetBinLabel(19,"fFillProfilesVsMUsingWeights");
 fIntFlowFlags->GetXaxis()->SetBinLabel(20,"fUseQvectorTerms");
 fIntFlowFlags->GetXaxis()->SetBinLabel(21,"fMaxCommonResultsHistogram");
 fIntFlowList->Add(fIntFlowFlags);

 // b) Book event-by-event quantities:
 // Re[Q_{m*n,k}], Im[Q_{m*n,k}] and S_{p,k}^M: 
 fReQ = new TMatrixD(12,9);
 fImQ = new TMatrixD(12,9);
 fSpk = new TMatrixD(8,9);
 // average correlations <2>, <4>, <6> and <8> for single event (bining is the same as in fIntFlowCorrelationsPro and fIntFlowCorrelationsHist):
 TString intFlowCorrelationsEBEName = "fIntFlowCorrelationsEBE";
 intFlowCorrelationsEBEName += fAnalysisLabel->Data();
 fIntFlowCorrelationsEBE = new TH1D(intFlowCorrelationsEBEName.Data(),intFlowCorrelationsEBEName.Data(),4,0,4);
 // weights for average correlations <2>, <4>, <6> and <8> for single event:
 TString intFlowEventWeightsForCorrelationsEBEName = "fIntFlowEventWeightsForCorrelationsEBE";
 intFlowEventWeightsForCorrelationsEBEName += fAnalysisLabel->Data();
 fIntFlowEventWeightsForCorrelationsEBE = new TH1D(intFlowEventWeightsForCorrelationsEBEName.Data(),intFlowEventWeightsForCorrelationsEBEName.Data(),4,0,4);
 // average all correlations for single event (bining is the same as in fIntFlowCorrelationsAllPro and fIntFlowCorrelationsAllHist):
 TString intFlowCorrelationsAllEBEName = "fIntFlowCorrelationsAllEBE";
 intFlowCorrelationsAllEBEName += fAnalysisLabel->Data();
 fIntFlowCorrelationsAllEBE = new TH1D(intFlowCorrelationsAllEBEName.Data(),intFlowCorrelationsAllEBEName.Data(),64,0,64);
 // average correction terms for non-uniform acceptance for single event 
 // (binning is the same as in fIntFlowCorrectionTermsForNUAPro[2] and fIntFlowCorrectionTermsForNUAHist[2]):
 TString fIntFlowCorrectionTermsForNUAEBEName = "fIntFlowCorrectionTermsForNUAEBE";
 fIntFlowCorrectionTermsForNUAEBEName += fAnalysisLabel->Data();
 for(Int_t sc=0;sc<2;sc++) // sin or cos terms
 {
  fIntFlowCorrectionTermsForNUAEBE[sc] = new TH1D(Form("%s: %s terms",fIntFlowCorrectionTermsForNUAEBEName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),4,0,4);  
 }
 // event weights for terms for non-uniform acceptance: 
 TString fIntFlowEventWeightForCorrectionTermsForNUAEBEName = "fIntFlowEventWeightForCorrectionTermsForNUAEBE";
 fIntFlowEventWeightForCorrectionTermsForNUAEBEName += fAnalysisLabel->Data();
 for(Int_t sc=0;sc<2;sc++) // sin or cos terms
 {
  fIntFlowEventWeightForCorrectionTermsForNUAEBE[sc] = new TH1D(Form("%s: %s terms",fIntFlowEventWeightForCorrectionTermsForNUAEBEName.Data(),sinCosFlag[sc].Data()),Form("Event weights for terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),4,0,4); // to be improved - 4  
 }
 // c) Book profiles: // to be improved (comment)
 // profile to hold average multiplicities and number of events for events with nRP>=0, nRP>=1, ... , and nRP>=8:
 TString avMultiplicityName = "fAvMultiplicity";
 avMultiplicityName += fAnalysisLabel->Data();
 fAvMultiplicity = new TProfile(avMultiplicityName.Data(),"Average multiplicities of reference particles (RPs)",9,0,9);
 fAvMultiplicity->SetTickLength(-0.01,"Y");
 fAvMultiplicity->SetMarkerStyle(25);
 fAvMultiplicity->SetLabelSize(0.05);
 fAvMultiplicity->SetLabelOffset(0.02,"Y");
 fAvMultiplicity->SetYTitle("Average multiplicity");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(1,"all evts");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(2,"n_{RP} #geq 1");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(3,"n_{RP} #geq 2");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(4,"n_{RP} #geq 3");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(5,"n_{RP} #geq 4");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(6,"n_{RP} #geq 5");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(7,"n_{RP} #geq 6");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(8,"n_{RP} #geq 7");
 (fAvMultiplicity->GetXaxis())->SetBinLabel(9,"n_{RP} #geq 8");
 fIntFlowProfiles->Add(fAvMultiplicity);
 // Average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with wrong errors!):
 TString correlationFlag[4] = {"#LT#LT2#GT#GT","#LT#LT4#GT#GT","#LT#LT6#GT#GT","#LT#LT8#GT#GT"};
 TString intFlowCorrelationsProName = "fIntFlowCorrelationsPro";
 intFlowCorrelationsProName += fAnalysisLabel->Data();
 fIntFlowCorrelationsPro = new TProfile(intFlowCorrelationsProName.Data(),"Average correlations for all events",4,0,4,"s");
 fIntFlowCorrelationsPro->Sumw2();
 fIntFlowCorrelationsPro->SetTickLength(-0.01,"Y");
 fIntFlowCorrelationsPro->SetMarkerStyle(25);
 fIntFlowCorrelationsPro->SetLabelSize(0.06);
 fIntFlowCorrelationsPro->SetLabelOffset(0.01,"Y");
 for(Int_t b=0;b<4;b++)
 {
  (fIntFlowCorrelationsPro->GetXaxis())->SetBinLabel(b+1,correlationFlag[b].Data());
 }
 fIntFlowProfiles->Add(fIntFlowCorrelationsPro);
 // Average correlations squared <<2>^2>, <<4>^2>, <<6>^2> and <<8>^2> for all events:
 TString squaredCorrelationFlag[4] = {"#LT#LT2#GT^{2}#GT","#LT#LT4#GT^{2}#GT","#LT#LT6#GT^{2}#GT","#LT#LT8#GT^{2}#GT"};
 TString intFlowSquaredCorrelationsProName = "fIntFlowSquaredCorrelationsPro";
 intFlowSquaredCorrelationsProName += fAnalysisLabel->Data();
 fIntFlowSquaredCorrelationsPro = new TProfile(intFlowSquaredCorrelationsProName.Data(),"Average squared correlations for all events",4,0,4,"s");
 fIntFlowSquaredCorrelationsPro->Sumw2();
 fIntFlowSquaredCorrelationsPro->SetTickLength(-0.01,"Y");
 fIntFlowSquaredCorrelationsPro->SetMarkerStyle(25);
 fIntFlowSquaredCorrelationsPro->SetLabelSize(0.06);
 fIntFlowSquaredCorrelationsPro->SetLabelOffset(0.01,"Y");
 for(Int_t b=0;b<4;b++)
 {
  (fIntFlowSquaredCorrelationsPro->GetXaxis())->SetBinLabel(b+1,squaredCorrelationFlag[b].Data());
 }
 fIntFlowProfiles->Add(fIntFlowSquaredCorrelationsPro);
 if(fCalculateCumulantsVsM)
 {
  for(Int_t ci=0;ci<4;ci++) // correlation index
  {
   // average correlations <<2>>, <<4>>, <<6>> and <<8>> versus multiplicity for all events (with wrong errors):
   TString intFlowCorrelationsVsMProName = "fIntFlowCorrelationsVsMPro";
   intFlowCorrelationsVsMProName += fAnalysisLabel->Data();
   fIntFlowCorrelationsVsMPro[ci] = new TProfile(Form("%s, %s",intFlowCorrelationsVsMProName.Data(),correlationFlag[ci].Data()),
                                                 Form("%s vs multiplicity",correlationFlag[ci].Data()),
                                                 fnBinsMult,fMinMult,fMaxMult,"s");   
   fIntFlowCorrelationsVsMPro[ci]->Sumw2();                                                                                       
   fIntFlowCorrelationsVsMPro[ci]->GetYaxis()->SetTitle(correlationFlag[ci].Data());
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("# RPs"); 
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("# POIs"); 
       } 
   fIntFlowProfiles->Add(fIntFlowCorrelationsVsMPro[ci]);
   // average squared correlations <<2>^2>, <<4>^2>, <<6>^2> and <<8>^2> versus multiplicity for all events:  
   TString intFlowSquaredCorrelationsVsMProName = "fIntFlowSquaredCorrelationsVsMPro";
   intFlowSquaredCorrelationsVsMProName += fAnalysisLabel->Data();
   fIntFlowSquaredCorrelationsVsMPro[ci] = new TProfile(Form("%s, %s",intFlowSquaredCorrelationsVsMProName.Data(),squaredCorrelationFlag[ci].Data()),
                                                        Form("%s vs multiplicity",squaredCorrelationFlag[ci].Data()),
                                                        fnBinsMult,fMinMult,fMaxMult,"s");   
   fIntFlowSquaredCorrelationsVsMPro[ci]->Sumw2();                                                                                              
   fIntFlowSquaredCorrelationsVsMPro[ci]->GetYaxis()->SetTitle(squaredCorrelationFlag[ci].Data());
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowSquaredCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("# RPs");
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowSquaredCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowSquaredCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowProfiles->Add(fIntFlowSquaredCorrelationsVsMPro[ci]);
  } // end of for(Int_t ci=0;ci<4;ci++) // correlation index  
 } // end of if(fCalculateCumulantsVsM)
 // averaged all correlations for all events (with wrong errors!):
 TString intFlowCorrelationsAllProName = "fIntFlowCorrelationsAllPro";
 intFlowCorrelationsAllProName += fAnalysisLabel->Data();
 fIntFlowCorrelationsAllPro = new TProfile(intFlowCorrelationsAllProName.Data(),"Average all correlations for all events",64,0,64);
 fIntFlowCorrelationsAllPro->Sumw2();
 fIntFlowCorrelationsAllPro->SetTickLength(-0.01,"Y");
 fIntFlowCorrelationsAllPro->SetMarkerStyle(25);
 fIntFlowCorrelationsAllPro->SetLabelSize(0.03);
 fIntFlowCorrelationsAllPro->SetLabelOffset(0.01,"Y");
 // 2-p correlations:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(1,"#LT#LT2#GT#GT_{n|n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(2,"#LT#LT2#GT#GT_{2n|2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(3,"#LT#LT2#GT#GT_{3n|3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(4,"#LT#LT2#GT#GT_{4n|4n}");
 // 3-p correlations:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(6,"#LT#LT3#GT#GT_{2n|n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(7,"#LT#LT3#GT#GT_{3n|2n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(8,"#LT#LT3#GT#GT_{4n|2n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(9,"#LT#LT3#GT#GT_{4n|3n,n}");
 // 4-p correlations:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(11,"#LT#LT4#GT#GT_{n,n|n,n}"); 
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(12,"#LT#LT4#GT#GT_{2n,n|2n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(13,"#LT#LT4#GT#GT_{2n,2n|2n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(14,"#LT#LT4#GT#GT_{3n|n,n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(15,"#LT#LT4#GT#GT_{3n,n|3n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(16,"#LT#LT4#GT#GT_{3n,n|2n,2n}"); 
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(17,"#LT#LT4#GT#GT_{4n|2n,n,n}");
 // 5-p correlations:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(19,"#LT#LT5#GT#GT_{2n,n|n,n,n}"); 
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(20,"#LT#LT5#GT#GT_{2n,2n|2n,n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(21,"#LT#LT5#GT#GT_{3n,n|2n,n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(22,"#LT#LT5#GT#GT_{4n|n,n,n,n}");
 // 6-p correlations:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(24,"#LT#LT6#GT#GT_{n,n,n|n,n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(25,"#LT#LT6#GT#GT_{2n,n,n|2n,n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(26,"#LT#LT6#GT#GT_{2n,2n|n,n,n,n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(27,"#LT#LT6#GT#GT_{3n,n|n,n,n,n}");
 // 7-p correlations:  
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(29,"#LT#LT7#GT#GT_{2n,n,n|n,n,n,n}");
 // 8-p correlations:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(31,"#LT#LT8#GT#GT_{n,n,n,n|n,n,n,n}");
 //  EXTRA correlations for v3{5} study:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(33,"#LT#LT4#GT#GT_{4n,2n|3n,3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(34,"#LT#LT5#GT#GT_{3n,3n|2n,2n,2n}");
 //  EXTRA correlations for Teaney-Yan study:
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(35,"#LT#LT2#GT#GT_{5n|5n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(36,"#LT#LT2#GT#GT_{6n|6n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(37,"#LT#LT3#GT#GT_{5n|3n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(38,"#LT#LT3#GT#GT_{5n|4n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(39,"#LT#LT3#GT#GT_{6n|3n,3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(40,"#LT#LT3#GT#GT_{6n|4n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(41,"#LT#LT3#GT#GT_{6n|5n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(42,"#LT#LT4#GT#GT_{6n|3n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(43,"#LT#LT4#GT#GT_{3n,2n|3n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(44,"#LT#LT4#GT#GT_{4n,1n|3n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(45,"#LT#LT4#GT#GT_{3n,3n|3n,3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(46,"#LT#LT4#GT#GT_{4n,2n|3n,3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(47,"#LT#LT4#GT#GT_{5n,1n|3n,3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(48,"#LT#LT4#GT#GT_{4n,2n|4n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(49,"#LT#LT4#GT#GT_{5n,1n|4n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(50,"#LT#LT4#GT#GT_{5n|3n,1n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(51,"#LT#LT4#GT#GT_{5n|2n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(52,"#LT#LT4#GT#GT_{5n,1n|5n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(53,"#LT#LT5#GT#GT_{3n,3n|3n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(54,"#LT#LT5#GT#GT_{4n,2n|3n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(55,"#LT#LT5#GT#GT_{3n,2n|3n,1n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(56,"#LT#LT5#GT#GT_{3n,2n|2n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(57,"#LT#LT5#GT#GT_{5n,1n|3n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(58,"#LT#LT6#GT#GT_{3n,2n,1n|3n,2n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(59,"#LT#LT4#GT#GT_{6n|4n,1n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(60,"#LT#LT4#GT#GT_{6n|2n,2n,2n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(61,"#LT#LT5#GT#GT_{6n|2n,2n,1n,1n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(62,"#LT#LT5#GT#GT_{4n,1n,1n|3n,3n}");
 (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(63,"#LT#LT6#GT#GT_{3n,3n|2n,2n,1n,1n}");
 fIntFlowProfiles->Add(fIntFlowCorrelationsAllPro);
 // average all correlations versus multiplicity (errors via Sumw2 - to be improved):
 if(fCalculateAllCorrelationsVsM)
 {
  // 2-p correlations vs M:  
  fIntFlowCorrelationsAllVsMPro[0] = new TProfile("two1n1n","#LT#LT2#GT#GT_{n|n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[1] = new TProfile("two2n2n","#LT#LT2#GT#GT_{2n|2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[2] = new TProfile("two3n3n","#LT#LT2#GT#GT_{3n|3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[3] = new TProfile("two4n4n","#LT#LT2#GT#GT_{4n|4n}",fnBinsMult,fMinMult,fMaxMult);
  // 3-p correlations vs M:
  fIntFlowCorrelationsAllVsMPro[5] = new TProfile("three2n1n1n","#LT#LT3#GT#GT_{2n|n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[6] = new TProfile("three3n2n1n","#LT#LT3#GT#GT_{3n|2n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[7] = new TProfile("three4n2n2n","#LT#LT3#GT#GT_{4n|2n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[8] = new TProfile("three4n3n1n","#LT#LT3#GT#GT_{4n|3n,n}",fnBinsMult,fMinMult,fMaxMult);
  // 4-p correlations vs M:
  fIntFlowCorrelationsAllVsMPro[10] = new TProfile("four1n1n1n1n","#LT#LT4#GT#GT_{n,n|n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[11] = new TProfile("four2n1n2n1n","#LT#LT4#GT#GT_{2n,n|2n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[12] = new TProfile("four2n2n2n2n","#LT#LT4#GT#GT_{2n,2n|2n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[13] = new TProfile("four3n1n1n1n","#LT#LT4#GT#GT_{3n|n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[14] = new TProfile("four3n1n3n1n","#LT#LT4#GT#GT_{3n,n|3n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[15] = new TProfile("four3n1n2n2n","#LT#LT4#GT#GT_{3n,n|2n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[16] = new TProfile("four4n2n1n1n","#LT#LT4#GT#GT_{4n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  // 5-p correlations vs M:
  fIntFlowCorrelationsAllVsMPro[18] = new TProfile("five2n1n1n1n1n","#LT#LT5#GT#GT_{2n,n|n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[19] = new TProfile("five2n2n2n1n1n","#LT#LT5#GT#GT_{2n,2n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[20] = new TProfile("five3n1n2n1n1n","#LT#LT5#GT#GT_{3n,n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[21] = new TProfile("five4n1n1n1n1n","#LT#LT5#GT#GT_{4n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  // 6-p correlations vs M:
  fIntFlowCorrelationsAllVsMPro[23] = new TProfile("six1n1n1n1n1n1n","#LT#LT6#GT#GT_{n,n,n|n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[24] = new TProfile("six2n1n1n2n1n1n","#LT#LT6#GT#GT_{2n,n,n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[25] = new TProfile("six2n2n1n1n1n1n","#LT#LT6#GT#GT_{2n,2n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[26] = new TProfile("six3n1n1n1n1n1n","#LT#LT6#GT#GT_{3n,n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  // 7-p correlations vs M:
  fIntFlowCorrelationsAllVsMPro[28] = new TProfile("seven2n1n1n1n1n1n1n","#LT#LT7#GT#GT_{2n,n,n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  // 8-p correlations vs M:
  fIntFlowCorrelationsAllVsMPro[30] = new TProfile("eight1n1n1n1n1n1n1n1n","#LT#LT8#GT#GT_{n,n,n,n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
  // EXTRA correlations vs M for v3{5} study (to be improved - put them in a right order somewhere):
  fIntFlowCorrelationsAllVsMPro[32] = new TProfile("four4n2n3n3n","#LT#LT4#GT#GT_{4n,2n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[33] = new TProfile("five3n3n2n2n2n","#LT#LT5#GT#GT_{3n,3n|2n,2n,2n}",fnBinsMult,fMinMult,fMaxMult);
  // EXTRA correlations vs M for Teaney-Yan study (to be improved - put them in a right order somewhere):
  fIntFlowCorrelationsAllVsMPro[34] = new TProfile("two5n5n","#LT#LT2#GT#GT_{5n|5n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[35] = new TProfile("two6n6n","#LT#LT2#GT#GT_{6n|6n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[36] = new TProfile("three5n3n2n","#LT#LT3#GT#GT_{5n|3n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[37] = new TProfile("three5n4n1n","#LT#LT3#GT#GT_{5n|4n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[38] = new TProfile("three6n3n3n","#LT#LT3#GT#GT_{6n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[39] = new TProfile("three6n4n2n","#LT#LT3#GT#GT_{6n|4n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[40] = new TProfile("three6n5n1n","#LT#LT3#GT#GT_{6n|5n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[41] = new TProfile("four6n3n2n1n","#LT#LT4#GT#GT_{6n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[42] = new TProfile("four3n2n3n2n","#LT#LT4#GT#GT_{3n,2n|3n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[43] = new TProfile("four4n1n3n2n","#LT#LT4#GT#GT_{4n,1n|3n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[44] = new TProfile("four3n3n3n3n","#LT#LT4#GT#GT_{3n,3n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[45] = new TProfile("four4n2n3n3n","#LT#LT4#GT#GT_{4n,2n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[46] = new TProfile("four5n1n3n3n","#LT#LT4#GT#GT_{5n,1n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[47] = new TProfile("four4n2n4n2n","#LT#LT4#GT#GT_{4n,2n|4n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[48] = new TProfile("four5n1n4n2n","#LT#LT4#GT#GT_{5n,1n|4n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[49] = new TProfile("four5n3n1n1n","#LT#LT4#GT#GT_{5n|3n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[50] = new TProfile("four5n2n2n1n","#LT#LT4#GT#GT_{5n|2n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[51] = new TProfile("four5n1n5n1n","#LT#LT4#GT#GT_{5n,1n|5n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[52] = new TProfile("five3n3n3n2n1n","#LT#LT5#GT#GT_{3n,3n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[53] = new TProfile("five4n2n3n2n1n","#LT#LT5#GT#GT_{4n,2n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[54] = new TProfile("five3n2n3n1n1n","#LT#LT5#GT#GT_{3n,2n|3n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[55] = new TProfile("five3n2n2n2n1n","#LT#LT5#GT#GT_{3n,2n|2n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[56] = new TProfile("five5n1n3n2n1n","#LT#LT5#GT#GT_{5n,1n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[57] = new TProfile("six3n2n1n3n2n1n","#LT#LT6#GT#GT_{3n,2n,1n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[58] = new TProfile("four6n4n1n1n","#LT#LT4#GT#GT_{6n|4n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[59] = new TProfile("four6n2n2n2n","#LT#LT4#GT#GT_{6n|2n,2n,2n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[60] = new TProfile("five6n2n2n1n1n","#LT#LT5#GT#GT_{6n|2n,2n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[61] = new TProfile("five4n1n1n3n3n","#LT#LT5#GT#GT_{4n,1n,1n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
  fIntFlowCorrelationsAllVsMPro[62] = new TProfile("six3n3n2n2n1n1n","#LT#LT6#GT#GT_{3n,3n|2n,2n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
  for(Int_t n=0;n<63;n++)
  {
   if(fIntFlowCorrelationsAllVsMPro[n])
   {
    fIntFlowCorrelationsAllVsMPro[n]->Sumw2();
    if(fMultiplicityIs==AliFlowCommonConstants::kRP)
    {
     fIntFlowCorrelationsAllVsMPro[n]->GetXaxis()->SetTitle("# RPs"); 
    } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
      {
       fIntFlowCorrelationsAllVsMPro[n]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
      } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
        {
         fIntFlowCorrelationsAllVsMPro[n]->GetXaxis()->SetTitle("# POIs");
        }
    fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[n]);
   } // end of if(fIntFlowCorrelationsAllVsMPro[n])
  } // end of for(Int_t n=0;n<63;n++)
 } // end of if(fCalculateAllCorrelationsVsM)
 // when particle weights are used some extra correlations appear:
 if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights) 
 {
  TString intFlowExtraCorrelationsProName = "fIntFlowExtraCorrelationsPro";
  intFlowExtraCorrelationsProName += fAnalysisLabel->Data();
  fIntFlowExtraCorrelationsPro = new TProfile(intFlowExtraCorrelationsProName.Data(),"Average extra correlations for all events",100,0,100,"s");
  fIntFlowExtraCorrelationsPro->SetTickLength(-0.01,"Y");
  fIntFlowExtraCorrelationsPro->SetMarkerStyle(25);
  fIntFlowExtraCorrelationsPro->SetLabelSize(0.03);
  fIntFlowExtraCorrelationsPro->SetLabelOffset(0.01,"Y");
  // extra 2-p correlations:
  (fIntFlowExtraCorrelationsPro->GetXaxis())->SetBinLabel(1,"<<w1^3 w2 cos(n*(phi1-phi2))>>");
  (fIntFlowExtraCorrelationsPro->GetXaxis())->SetBinLabel(2,"<<w1 w2 w3^2 cos(n*(phi1-phi2))>>");
  fIntFlowProfiles->Add(fIntFlowExtraCorrelationsPro);
 } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights)
 // average product of correlations <2>, <4>, <6> and <8>:  
 TString productFlag[6] = {"#LT#LT2#GT#LT4#GT#GT","#LT#LT2#GT#LT6#GT#GT","#LT#LT2#GT#LT8#GT#GT",
                           "#LT#LT4#GT#LT6#GT#GT","#LT#LT4#GT#LT8#GT#GT","#LT#LT6#GT#LT8#GT#GT"};
 TString intFlowProductOfCorrelationsProName = "fIntFlowProductOfCorrelationsPro";
 intFlowProductOfCorrelationsProName += fAnalysisLabel->Data();
 fIntFlowProductOfCorrelationsPro = new TProfile(intFlowProductOfCorrelationsProName.Data(),"Average products of correlations",6,0,6);
 fIntFlowProductOfCorrelationsPro->SetTickLength(-0.01,"Y");
 fIntFlowProductOfCorrelationsPro->SetMarkerStyle(25); 
 fIntFlowProductOfCorrelationsPro->SetLabelSize(0.05);
 fIntFlowProductOfCorrelationsPro->SetLabelOffset(0.01,"Y");
 for(Int_t b=0;b<6;b++)
 {
  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(b+1,productFlag[b].Data());
 }
 fIntFlowProfiles->Add(fIntFlowProductOfCorrelationsPro); 
 // average product of correlations <2>, <4>, <6> and <8> versus multiplicity
 // [0=<<2><4>>,1=<<2><6>>,2=<<2><8>>,3=<<4><6>>,4=<<4><8>>,5=<<6><8>>]  
 if(fCalculateCumulantsVsM)
 {
  TString intFlowProductOfCorrelationsVsMProName = "fIntFlowProductOfCorrelationsVsMPro";
  intFlowProductOfCorrelationsVsMProName += fAnalysisLabel->Data();
  for(Int_t pi=0;pi<6;pi++)
  { 
   fIntFlowProductOfCorrelationsVsMPro[pi] = new TProfile(Form("%s, %s",intFlowProductOfCorrelationsVsMProName.Data(),productFlag[pi].Data()),
                                                          Form("%s versus multiplicity",productFlag[pi].Data()),
                                                          fnBinsMult,fMinMult,fMaxMult);             
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowProductOfCorrelationsVsMPro[pi]->GetXaxis()->SetTitle("# RPs");
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowProductOfCorrelationsVsMPro[pi]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowProductOfCorrelationsVsMPro[pi]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowProfiles->Add(fIntFlowProductOfCorrelationsVsMPro[pi]);
  } // end of for(Int_t pi=0;pi<6;pi++)
 } // end of if(fCalculateCumulantsVsM) 
 // average product of correction terms for NUA:  
 TString intFlowProductOfCorrectionTermsForNUAProName = "fIntFlowProductOfCorrectionTermsForNUAPro";
 intFlowProductOfCorrectionTermsForNUAProName += fAnalysisLabel->Data();
 fIntFlowProductOfCorrectionTermsForNUAPro = new TProfile(intFlowProductOfCorrectionTermsForNUAProName.Data(),"Average products of correction terms for NUA",27,0,27);
 fIntFlowProductOfCorrectionTermsForNUAPro->SetTickLength(-0.01,"Y");
 fIntFlowProductOfCorrectionTermsForNUAPro->SetMarkerStyle(25); 
 fIntFlowProductOfCorrectionTermsForNUAPro->SetLabelSize(0.03);
 fIntFlowProductOfCorrectionTermsForNUAPro->SetLabelOffset(0.01,"Y");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(1,"<<2><cos(#phi)>>");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(2,"<<2><sin(#phi)>>");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(3,"<<cos(#phi)><sin(#phi)>>");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(4,"Cov(<2>,<cos(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(5,"Cov(<2>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(6,"Cov(<2>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(7,"Cov(<2>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(8,"Cov(<4>,<cos(#phi)>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(9,"Cov(<4>,<sin(#phi)>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(10,"Cov(<4>,<cos(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(11,"Cov(<4>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(12,"Cov(<4>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(13,"Cov(<4>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(14,"Cov(<cos(#phi)>,<cos(#phi_{1}+#phi_{2})>)"); 
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(15,"Cov(<cos(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(16,"Cov(<cos(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(17,"Cov(<cos(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(18,"Cov(<sin(#phi)>,<cos(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(19,"Cov(<sin(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(20,"Cov(<sin(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(21,"Cov(<sin(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(22,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(23,"Cov(<cos(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(24,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(25,"Cov(<sin(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(26,"Cov(<sin(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(27,"Cov(<cos(#phi_{1}-#phi_{2}-#phi_{3}>,<sin(#phi_{1}-#phi_{2}-#phi_{3}>)");
 fIntFlowProfiles->Add(fIntFlowProductOfCorrectionTermsForNUAPro);
 // average correction terms for non-uniform acceptance (with wrong errors!):
 for(Int_t sc=0;sc<2;sc++) // sin or cos terms
 {
  TString intFlowCorrectionTermsForNUAProName = "fIntFlowCorrectionTermsForNUAPro";
  intFlowCorrectionTermsForNUAProName += fAnalysisLabel->Data();
  fIntFlowCorrectionTermsForNUAPro[sc] = new TProfile(Form("%s: %s terms",intFlowCorrectionTermsForNUAProName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),4,0,4,"s");
  fIntFlowCorrectionTermsForNUAPro[sc]->SetTickLength(-0.01,"Y");
  fIntFlowCorrectionTermsForNUAPro[sc]->SetMarkerStyle(25);
  fIntFlowCorrectionTermsForNUAPro[sc]->SetLabelSize(0.05);
  fIntFlowCorrectionTermsForNUAPro[sc]->SetLabelOffset(0.01,"Y");
  (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(1,Form("#LT#LT%s(n(#phi_{1}))#GT#GT",sinCosFlag[sc].Data()));
  (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(2,Form("#LT#LT%s(n(#phi_{1}+#phi_{2}))#GT#GT",sinCosFlag[sc].Data()));  
  (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(3,Form("#LT#LT%s(n(#phi_{1}-#phi_{2}-#phi_{3}))#GT#GT",sinCosFlag[sc].Data()));  
  (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(4,Form("#LT#LT%s(n(2#phi_{1}-#phi_{2}))#GT#GT",sinCosFlag[sc].Data()));  
  fIntFlowProfiles->Add(fIntFlowCorrectionTermsForNUAPro[sc]);
  // versus multiplicity:
  if(fCalculateCumulantsVsM)
  {
   TString correctionTermFlag[4] = {"(n(phi1))","(n(phi1+phi2))","(n(phi1-phi2-phi3))","(n(2phi1-phi2))"}; // to be improved - hardwired 4
   for(Int_t ci=0;ci<4;ci++) // correction term index (to be improved - hardwired 4)
   {
    TString intFlowCorrectionTermsForNUAVsMProName = "fIntFlowCorrectionTermsForNUAVsMPro";
    intFlowCorrectionTermsForNUAVsMProName += fAnalysisLabel->Data();
    fIntFlowCorrectionTermsForNUAVsMPro[sc][ci] = new TProfile(Form("%s: #LT#LT%s%s#GT#GT",intFlowCorrectionTermsForNUAVsMProName.Data(),sinCosFlag[sc].Data(),correctionTermFlag[ci].Data()),Form("#LT#LT%s%s#GT#GT vs M",sinCosFlag[sc].Data(),correctionTermFlag[ci].Data()),fnBinsMult,fMinMult,fMaxMult,"s");
    fIntFlowProfiles->Add(fIntFlowCorrectionTermsForNUAVsMPro[sc][ci]);
   }
  } // end of if(fCalculateCumulantsVsM)
 } // end of for(Int_t sc=0;sc<2;sc++) 
 
 // d) Book histograms holding the final results:
 // average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with correct errors!):
 TString intFlowCorrelationsHistName = "fIntFlowCorrelationsHist";
 intFlowCorrelationsHistName += fAnalysisLabel->Data();
 fIntFlowCorrelationsHist = new TH1D(intFlowCorrelationsHistName.Data(),"Average correlations for all events",4,0,4);
 fIntFlowCorrelationsHist->SetTickLength(-0.01,"Y");
 fIntFlowCorrelationsHist->SetMarkerStyle(25);
 fIntFlowCorrelationsHist->SetLabelSize(0.06);
 fIntFlowCorrelationsHist->SetLabelOffset(0.01,"Y");
 (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(1,"#LT#LT2#GT#GT");
 (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(2,"#LT#LT4#GT#GT");
 (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(3,"#LT#LT6#GT#GT");
 (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(4,"#LT#LT8#GT#GT");
 fIntFlowResults->Add(fIntFlowCorrelationsHist);
 // average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with correct errors!) vs M:
 if(fCalculateCumulantsVsM)
 {
  for(Int_t ci=0;ci<4;ci++) // correlation index
  {
   TString intFlowCorrelationsVsMHistName = "fIntFlowCorrelationsVsMHist";
   intFlowCorrelationsVsMHistName += fAnalysisLabel->Data();
   fIntFlowCorrelationsVsMHist[ci] = new TH1D(Form("%s, %s",intFlowCorrelationsVsMHistName.Data(),correlationFlag[ci].Data()),
                                              Form("%s vs multiplicity",correlationFlag[ci].Data()),
                                              fnBinsMult,fMinMult,fMaxMult);                                            
   fIntFlowCorrelationsVsMHist[ci]->GetYaxis()->SetTitle(correlationFlag[ci].Data());
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowCorrelationsVsMHist[ci]->GetXaxis()->SetTitle("# RPs");
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowCorrelationsVsMHist[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowCorrelationsVsMHist[ci]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowResults->Add(fIntFlowCorrelationsVsMHist[ci]);
  } // end of for(Int_t ci=0;ci<4;ci++) // correlation index   
 } // end of if(fCalculateCumulantsVsM) 
 // average all correlations for all events (with correct errors!):
 TString intFlowCorrelationsAllHistName = "fIntFlowCorrelationsAllHist";
 intFlowCorrelationsAllHistName += fAnalysisLabel->Data();
 fIntFlowCorrelationsAllHist = new TH1D(intFlowCorrelationsAllHistName.Data(),"Average correlations for all events",34,0,34);
 fIntFlowCorrelationsAllHist->SetTickLength(-0.01,"Y");
 fIntFlowCorrelationsAllHist->SetMarkerStyle(25);
 fIntFlowCorrelationsAllHist->SetLabelSize(0.03);
 fIntFlowCorrelationsAllHist->SetLabelOffset(0.01,"Y");
 // 2-p correlations:
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
 // 3-p correlations:
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
 // 4-p correlations:
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
 // 5-p correlations:
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
 // 6-p correlations:
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
 // 7-p correlations:  
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
 // 8-p correlations:
 (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
 fIntFlowResults->Add(fIntFlowCorrelationsAllHist);
 // average correction terms for non-uniform acceptance (with correct errors!):
 for(Int_t sc=0;sc<2;sc++) // sin or cos terms
 {
  TString intFlowCorrectionTermsForNUAHistName = "fIntFlowCorrectionTermsForNUAHist";
  intFlowCorrectionTermsForNUAHistName += fAnalysisLabel->Data();
  fIntFlowCorrectionTermsForNUAHist[sc] = new TH1D(Form("%s: %s terms",intFlowCorrectionTermsForNUAHistName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),4,0,4);
  fIntFlowCorrectionTermsForNUAHist[sc]->SetTickLength(-0.01,"Y");
  fIntFlowCorrectionTermsForNUAHist[sc]->SetMarkerStyle(25);
  fIntFlowCorrectionTermsForNUAHist[sc]->SetLabelSize(0.05);
  fIntFlowCorrectionTermsForNUAHist[sc]->SetLabelOffset(0.01,"Y");
  (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(1,Form("#LT#LT%s(n(#phi_{1}))#GT#GT",sinCosFlag[sc].Data()));
  (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(2,Form("#LT#LT%s(n(#phi_{1}+#phi_{2}))#GT#GT",sinCosFlag[sc].Data()));  
  (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(3,Form("#LT#LT%s(n(#phi_{1}-#phi_{2}-#phi_{3}))#GT#GT",sinCosFlag[sc].Data()));  
  (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(4,Form("#LT#LT%s(n(2#phi_{1}-#phi_{2}))#GT#GT",sinCosFlag[sc].Data()));   
  fIntFlowResults->Add(fIntFlowCorrectionTermsForNUAHist[sc]);
 } // end of for(Int_t sc=0;sc<2;sc++) 
 // covariances (multiplied with weight dependent prefactor):
 TString intFlowCovariancesName = "fIntFlowCovariances";
 intFlowCovariancesName += fAnalysisLabel->Data();
 fIntFlowCovariances = new TH1D(intFlowCovariancesName.Data(),"Covariances (multiplied with weight dependent prefactor)",6,0,6);
 fIntFlowCovariances->SetLabelSize(0.04);
 fIntFlowCovariances->SetMarkerStyle(25);
 (fIntFlowCovariances->GetXaxis())->SetBinLabel(1,"Cov(#LT2#GT,#LT4#GT)");
 (fIntFlowCovariances->GetXaxis())->SetBinLabel(2,"Cov(#LT2#GT,#LT6#GT)");
 (fIntFlowCovariances->GetXaxis())->SetBinLabel(3,"Cov(#LT2#GT,#LT8#GT)");
 (fIntFlowCovariances->GetXaxis())->SetBinLabel(4,"Cov(#LT4#GT,#LT6#GT)");
 (fIntFlowCovariances->GetXaxis())->SetBinLabel(5,"Cov(#LT4#GT,#LT8#GT)");
 (fIntFlowCovariances->GetXaxis())->SetBinLabel(6,"Cov(#LT6#GT,#LT8#GT)");  
 fIntFlowResults->Add(fIntFlowCovariances);
 // sum of linear and quadratic event weights for <2>, <4>, <6> and <8>:
 TString intFlowSumOfEventWeightsName = "fIntFlowSumOfEventWeights";
 intFlowSumOfEventWeightsName += fAnalysisLabel->Data();
 for(Int_t power=0;power<2;power++)
 {
  fIntFlowSumOfEventWeights[power] = new TH1D(Form("%s: %s",intFlowSumOfEventWeightsName.Data(),powerFlag[power].Data()),Form("Sum of %s event weights for correlations",powerFlag[power].Data()),4,0,4);
  fIntFlowSumOfEventWeights[power]->SetLabelSize(0.04);
  fIntFlowSumOfEventWeights[power]->SetMarkerStyle(25);
  if(power == 0)
  {
   (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{#LT2#GT}");
   (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{#LT4#GT}");
   (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{#LT6#GT}");
   (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{#LT8#GT}");
  } else if (power == 1) 
    {
     (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{#LT2#GT}^{2}");
     (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{#LT4#GT}^{2}");
     (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{#LT6#GT}^{2}");
     (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{#LT8#GT}^{2}");
    }
  fIntFlowResults->Add(fIntFlowSumOfEventWeights[power]);
 } 
 // sum of products of event weights for correlations <2>, <4>, <6> and <8>:  
 TString intFlowSumOfProductOfEventWeightsName = "fIntFlowSumOfProductOfEventWeights";
 intFlowSumOfProductOfEventWeightsName += fAnalysisLabel->Data();
 fIntFlowSumOfProductOfEventWeights = new TH1D(intFlowSumOfProductOfEventWeightsName.Data(),"Sum of product of event weights for correlations",6,0,6);
 fIntFlowSumOfProductOfEventWeights->SetLabelSize(0.04);
 fIntFlowSumOfProductOfEventWeights->SetMarkerStyle(25);
 (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LT4#GT}");
 (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LT6#GT}");
 (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LT8#GT}");
 (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LT6#GT}");
 (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(5,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LT8#GT}");
 (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(6,"#sum_{i=1}^{N} w_{#LT6#GT} w_{#LT8#GT}");
 fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeights);
 // final result for covariances of correlations (multiplied with weight dependent prefactor) versus M
 // [0=Cov(2,4),1=Cov(2,6),2=Cov(2,8),3=Cov(4,6),4=Cov(4,8),5=Cov(6,8)]:
 if(fCalculateCumulantsVsM)
 {
  TString intFlowCovariancesVsMName = "fIntFlowCovariancesVsM";
  intFlowCovariancesVsMName += fAnalysisLabel->Data();
  TString covarianceFlag[6] = {"Cov(<2>,<4>)","Cov(<2>,<6>)","Cov(<2>,<8>)","Cov(<4>,<6>)","Cov(<4>,<8>)","Cov(<6>,<8>)"};
  for(Int_t ci=0;ci<6;ci++)
  {
   fIntFlowCovariancesVsM[ci] = new TH1D(Form("%s, %s",intFlowCovariancesVsMName.Data(),covarianceFlag[ci].Data()),
                                         Form("%s vs multiplicity",covarianceFlag[ci].Data()),
                                         fnBinsMult,fMinMult,fMaxMult);
   fIntFlowCovariancesVsM[ci]->GetYaxis()->SetTitle(covarianceFlag[ci].Data());
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowCovariancesVsM[ci]->GetXaxis()->SetTitle("# RPs");
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowCovariancesVsM[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowCovariancesVsM[ci]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowResults->Add(fIntFlowCovariancesVsM[ci]);
  }
 } // end of if(fCalculateCumulantsVsM) 
 // sum of linear and quadratic event weights for <2>, <4>, <6> and <8> versus multiplicity
 // [0=sum{w_{<2>}},1=sum{w_{<4>}},2=sum{w_{<6>}},3=sum{w_{<8>}}][0=linear 1,1=quadratic]:
 if(fCalculateCumulantsVsM)
 {
  TString intFlowSumOfEventWeightsVsMName = "fIntFlowSumOfEventWeightsVsM";
  intFlowSumOfEventWeightsVsMName += fAnalysisLabel->Data();
  TString sumFlag[2][4] = {{"#sum_{i=1}^{N} w_{<2>}","#sum_{i=1}^{N} w_{<4>}","#sum_{i=1}^{N} w_{<6>}","#sum_{i=1}^{N} w_{<8>}"},
                           {"#sum_{i=1}^{N} w_{<2>}^{2}","#sum_{i=1}^{N} w_{<4>}^{2}","#sum_{i=1}^{N} w_{<6>}^{2}","#sum_{i=1}^{N} w_{<8>}^{2}"}};
  for(Int_t si=0;si<4;si++)
  {
   for(Int_t power=0;power<2;power++)
   {
    fIntFlowSumOfEventWeightsVsM[si][power] = new TH1D(Form("%s, %s",intFlowSumOfEventWeightsVsMName.Data(),sumFlag[power][si].Data()),
                                                       Form("%s vs multiplicity",sumFlag[power][si].Data()),
                                                       fnBinsMult,fMinMult,fMaxMult);    
    fIntFlowSumOfEventWeightsVsM[si][power]->GetYaxis()->SetTitle(sumFlag[power][si].Data());  
    if(fMultiplicityIs==AliFlowCommonConstants::kRP)
    {
     fIntFlowSumOfEventWeightsVsM[si][power]->GetXaxis()->SetTitle("# RPs");
    } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
      {
       fIntFlowSumOfEventWeightsVsM[si][power]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
      } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
        {
         fIntFlowSumOfEventWeightsVsM[si][power]->GetXaxis()->SetTitle("# POIs");
        }
    fIntFlowResults->Add(fIntFlowSumOfEventWeightsVsM[si][power]);
   } // end of for(Int_t power=0;power<2;power++)
  } // end of for(Int_t si=0;si<4;si++)   
 } // end of if(fCalculateCumulantsVsM)
 // sum of products of event weights for correlations <2>, <4>, <6> and <8> vs M
 // [0=sum{w_{<2>}w_{<4>}},1=sum{w_{<2>}w_{<6>}},2=sum{w_{<2>}w_{<8>}},
 //  3=sum{w_{<4>}w_{<6>}},4=sum{w_{<4>}w_{<8>}},5=sum{w_{<6>}w_{<8>}}]:  
 if(fCalculateCumulantsVsM)
 {
  TString intFlowSumOfProductOfEventWeightsVsMName = "fIntFlowSumOfProductOfEventWeightsVsM";
  intFlowSumOfProductOfEventWeightsVsMName += fAnalysisLabel->Data();
  TString sopowFlag[6] = {"#sum_{i=1}^{N} w_{<2>} w_{<4>}","#sum_{i=1}^{N} w_{<2>} w_{<6>}","#sum_{i=1}^{N} w_{<2>} w_{<8>}",
                          "#sum_{i=1}^{N} w_{<4>} w_{<6>}","#sum_{i=1}^{N} w_{<4>} w_{<8>}","#sum_{i=1}^{N} w_{<6>} w_{<8>}"}; 
  for(Int_t pi=0;pi<6;pi++)
  {
   fIntFlowSumOfProductOfEventWeightsVsM[pi] = new TH1D(Form("%s, %s",intFlowSumOfProductOfEventWeightsVsMName.Data(),sopowFlag[pi].Data()),
                                                        Form("%s versus multiplicity",sopowFlag[pi].Data()),
                                                        fnBinsMult,fMinMult,fMaxMult); 
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetXaxis()->SetTitle("# RPs"); 
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetYaxis()->SetTitle(sopowFlag[pi].Data()); 
   fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeightsVsM[pi]);
  } // end of for(Int_t pi=0;pi<6;pi++) 
 } // end of if(fCalculateCumulantsVsM)
 // covariances of NUA terms (multiplied with weight dependent prefactor):
 TString intFlowCovariancesNUAName = "fIntFlowCovariancesNUA";
 intFlowCovariancesNUAName += fAnalysisLabel->Data();
 fIntFlowCovariancesNUA = new TH1D(intFlowCovariancesNUAName.Data(),"Covariances for NUA (multiplied with weight dependent prefactor)",27,0,27);
 fIntFlowCovariancesNUA->SetLabelSize(0.04);
 fIntFlowCovariancesNUA->SetMarkerStyle(25);
 fIntFlowCovariancesNUA->GetXaxis()->SetLabelSize(0.02);
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(1,"Cov(<2>,<cos(#phi)>");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(2,"Cov(<2>,<sin(#phi)>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(3,"Cov(<cos(#phi)>,<sin(#phi)>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(4,"Cov(<2>,<cos(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(5,"Cov(<2>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(6,"Cov(<2>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(7,"Cov(<2>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(8,"Cov(<4>,<cos(#phi)>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(9,"Cov(<4>,<sin(#phi)>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(10,"Cov(<4>,<cos(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(11,"Cov(<4>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(12,"Cov(<4>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(13,"Cov(<4>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(14,"Cov(<cos(#phi)>,<cos(#phi_{1}+#phi_{2})>)"); 
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(15,"Cov(<cos(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(16,"Cov(<cos(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(17,"Cov(<cos(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(18,"Cov(<sin(#phi)>,<cos(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(19,"Cov(<sin(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(20,"Cov(<sin(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(21,"Cov(<sin(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(22,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}+#phi_{2})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(23,"Cov(<cos(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(24,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(25,"Cov(<sin(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(26,"Cov(<sin(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
 (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(27,"Cov(<cos(#phi_{1}-#phi_{2}-#phi_{3}>,<sin(#phi_{1}-#phi_{2}-#phi_{3}>)");
 fIntFlowResults->Add(fIntFlowCovariancesNUA);
 // sum of linear and quadratic event weights for NUA terms:
 TString intFlowSumOfEventWeightsNUAName = "fIntFlowSumOfEventWeightsNUA";
 intFlowSumOfEventWeightsNUAName += fAnalysisLabel->Data();
 for(Int_t sc=0;sc<2;sc++)
 {
  for(Int_t power=0;power<2;power++)
  {
   fIntFlowSumOfEventWeightsNUA[sc][power] = new TH1D(Form("%s: %s, %s",intFlowSumOfEventWeightsNUAName.Data(),powerFlag[power].Data(),sinCosFlag[sc].Data()),Form("Sum of %s event weights for NUA %s terms",powerFlag[power].Data(),sinCosFlag[sc].Data()),4,0,4); // to be improved - 4
   fIntFlowSumOfEventWeightsNUA[sc][power]->SetLabelSize(0.05);
   fIntFlowSumOfEventWeightsNUA[sc][power]->SetMarkerStyle(25);
   if(power == 0)
   {
    (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(1,Form("#sum_{i=1}^{N} w_{<%s(#phi)>}",sinCosFlag[sc].Data()));
    (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(2,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}+#phi_{2})>}",sinCosFlag[sc].Data()));
    (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(3,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}-#phi_{2}-#phi_{3})>}",sinCosFlag[sc].Data()));   
    (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(4,Form("#sum_{i=1}^{N} w_{<%s(2#phi_{1}-#phi_{2})>}",sinCosFlag[sc].Data()));
   } else if(power == 1) 
     {
      (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(1,Form("#sum_{i=1}^{N} w_{<%s(#phi)>}^{2}",sinCosFlag[sc].Data()));
      (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(2,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}+#phi_{2})>}^{2}",sinCosFlag[sc].Data()));
      (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(3,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}-#phi_{2}-#phi_{3})>}^{2}",sinCosFlag[sc].Data()));
      (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(4,Form("#sum_{i=1}^{N} w_{<%s(2#phi_{1}-#phi_{2})>}^{2}",sinCosFlag[sc].Data()));
     }
   fIntFlowResults->Add(fIntFlowSumOfEventWeightsNUA[sc][power]);
  }
 }  
 // sum of products of event weights for NUA terms:  
 TString intFlowSumOfProductOfEventWeightsNUAName = "fIntFlowSumOfProductOfEventWeightsNUA";
 intFlowSumOfProductOfEventWeightsNUAName += fAnalysisLabel->Data();
 fIntFlowSumOfProductOfEventWeightsNUA = new TH1D(intFlowSumOfProductOfEventWeightsNUAName.Data(),"Sum of product of event weights for NUA terms",27,0,27);
 fIntFlowSumOfProductOfEventWeightsNUA->SetLabelSize(0.02);
 fIntFlowSumOfProductOfEventWeightsNUA->SetMarkerStyle(25);
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LTcos(#phi)#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LTsin(#phi)#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{#LTcos(#phi)#GT} w_{#LTsin(#phi)#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LTcos(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(5,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LTsin(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(6,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(7,"#sum_{i=1}^{N} w_{#LT2#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(8,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LTcos(#phi)#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(9,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LTsin(#phi)#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(10,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LTcos(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(11,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LTsin(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(12,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(13,"#sum_{i=1}^{N} w_{#LT4#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(14,"#sum_{i=1}^{N} w_{#LTcos(#phi)#GT} w_{#LTcos(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(15,"#sum_{i=1}^{N} w_{#LTcos(#phi)#GT} w_{#LTsin(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(16,"#sum_{i=1}^{N} w_{#LTcos(#phi)#GT} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(17,"#sum_{i=1}^{N} w_{#LTcos(#phi)#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(18,"#sum_{i=1}^{N} w_{#LTsin(#phi)#GT} w_{#LTcos(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(19,"#sum_{i=1}^{N} w_{#LTsin(#phi)#GT} w_{#LTsin(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(20,"#sum_{i=1}^{N} w_{#LTsin(#phi)#GT} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(21,"#sum_{i=1}^{N} w_{#LTsin(#phi)#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(22,"#sum_{i=1}^{N} w_{#LTcos(#phi_{1}+#phi_{2})#GT} w_{#LTsin(#phi_{1}+#phi_{2})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(23,"#sum_{i=1}^{N} w_{#LTcos(#phi_{1}+#phi_{2})#GT} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(24,"#sum_{i=1}^{N} w_{#LTcos(#phi_{1}+#phi_{2})#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(25,"#sum_{i=1}^{N} w_{#LTsin(#phi_{1}+#phi_{2})#GT} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(26,"#sum_{i=1}^{N} w_{#LTsin(#phi_{1}+#phi_{2})#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(27,"#sum_{i=1}^{N} w_{#LTcos(#phi_{1}-#phi_{2}-#phi_{3})#GT} w_{#LTsin(#phi_{1}-#phi_{2}-#phi_{3})#GT}");
 fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeightsNUA);
 // Final results for reference Q-cumulants:
 TString cumulantFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
 TString intFlowQcumulantsName = "fIntFlowQcumulants";
 intFlowQcumulantsName += fAnalysisLabel->Data();
 fIntFlowQcumulants = new TH1D(intFlowQcumulantsName.Data(),"Reference Q-cumulants",4,0,4);
 if(fPropagateErrorAlsoFromNIT)
 {
  fIntFlowQcumulants->SetTitle("Reference Q-cumulants (error from non-isotropic terms also propagated)");
 }
 fIntFlowQcumulants->SetLabelSize(0.05);
 fIntFlowQcumulants->SetMarkerStyle(25);
 for(Int_t b=0;b<4;b++)
 {
  (fIntFlowQcumulants->GetXaxis())->SetBinLabel(b+1,cumulantFlag[b].Data());
 } 
 fIntFlowResults->Add(fIntFlowQcumulants);
 // Final results for reference Q-cumulants rebinned in M: 
 if(fCalculateCumulantsVsM)
 {
  TString intFlowQcumulantsRebinnedInMName = "fIntFlowQcumulantsRebinnedInM";
  intFlowQcumulantsRebinnedInMName += fAnalysisLabel->Data();
  fIntFlowQcumulantsRebinnedInM = new TH1D(intFlowQcumulantsRebinnedInMName.Data(),"Reference Q-cumulants rebinned in M",4,0,4);
  fIntFlowQcumulantsRebinnedInM->SetLabelSize(0.05);
  fIntFlowQcumulantsRebinnedInM->SetMarkerStyle(25);
  for(Int_t b=0;b<4;b++)
  {
   (fIntFlowQcumulantsRebinnedInM->GetXaxis())->SetBinLabel(b+1,cumulantFlag[b].Data());
  } 
  fIntFlowResults->Add(fIntFlowQcumulantsRebinnedInM);
 } // end of if(fCalculateCumulantsVsM) 
 // Ratio between error squared: with/without non-isotropic terms:
 TString intFlowQcumulantsErrorSquaredRatioName = "fIntFlowQcumulantsErrorSquaredRatio";
 intFlowQcumulantsErrorSquaredRatioName += fAnalysisLabel->Data();
 fIntFlowQcumulantsErrorSquaredRatio = new TH1D(intFlowQcumulantsErrorSquaredRatioName.Data(),"Error squared of reference Q-cumulants: #frac{with NUA terms}{without NUA terms}",4,0,4);
 fIntFlowQcumulantsErrorSquaredRatio->SetLabelSize(0.05);
 fIntFlowQcumulantsErrorSquaredRatio->SetMarkerStyle(25);
 for(Int_t b=0;b<4;b++)
 {
  (fIntFlowQcumulantsErrorSquaredRatio->GetXaxis())->SetBinLabel(b+1,cumulantFlag[b].Data());
 } 
 fIntFlowResults->Add(fIntFlowQcumulantsErrorSquaredRatio);
 // final results for integrated Q-cumulants versus multiplicity:
 if(fCalculateCumulantsVsM)
 {
  TString intFlowQcumulantsVsMName = "fIntFlowQcumulantsVsM";
  intFlowQcumulantsVsMName += fAnalysisLabel->Data();
  for(Int_t co=0;co<4;co++) // cumulant order
  {
   fIntFlowQcumulantsVsM[co] = new TH1D(Form("%s, %s",intFlowQcumulantsVsMName.Data(),cumulantFlag[co].Data()),
                                        Form("%s vs multiplicity",cumulantFlag[co].Data()),
                                        fnBinsMult,fMinMult,fMaxMult);
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowQcumulantsVsM[co]->GetXaxis()->SetTitle("# RPs");
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowQcumulantsVsM[co]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowQcumulantsVsM[co]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowQcumulantsVsM[co]->GetYaxis()->SetTitle(cumulantFlag[co].Data());  
   fIntFlowResults->Add(fIntFlowQcumulantsVsM[co]);                                    
  } // end of for(Int_t co=0;co<4;co++) // cumulant order
 } // end of if(fCalculateCumulantsVsM)
 // final integrated flow estimates from Q-cumulants:
 TString flowFlag[4] = {Form("v_{%d}{2,QC}",fHarmonic),Form("v_{%d}{4,QC}",fHarmonic),Form("v_{%d}{6,QC}",fHarmonic),Form("v_{%d}{8,QC}",fHarmonic)};
 TString intFlowName = "fIntFlow";
 intFlowName += fAnalysisLabel->Data();  
 // integrated flow from Q-cumulants:
 fIntFlow = new TH1D(intFlowName.Data(),"Reference flow estimates from Q-cumulants",4,0,4);
 fIntFlow->SetLabelSize(0.05);
 fIntFlow->SetMarkerStyle(25);
 for(Int_t b=0;b<4;b++)
 {
  (fIntFlow->GetXaxis())->SetBinLabel(b+1,flowFlag[b].Data()); 
 }
 fIntFlowResults->Add(fIntFlow); 
 // Reference flow vs M rebinned in one huge bin:
 if(fCalculateCumulantsVsM)
 { 
  TString intFlowRebinnedInMName = "fIntFlowRebinnedInM";
  intFlowRebinnedInMName += fAnalysisLabel->Data();  
  fIntFlowRebinnedInM = new TH1D(intFlowRebinnedInMName.Data(),"Reference flow estimates from Q-cumulants (rebinned in M)",4,0,4);
  fIntFlowRebinnedInM->SetLabelSize(0.05);
  fIntFlowRebinnedInM->SetMarkerStyle(25);
  for(Int_t b=0;b<4;b++)
  {
   (fIntFlowRebinnedInM->GetXaxis())->SetBinLabel(b+1,flowFlag[b].Data()); 
  }
  fIntFlowResults->Add(fIntFlowRebinnedInM); 
 } 
 // integrated flow from Q-cumulants: versus multiplicity:
 if(fCalculateCumulantsVsM)
 {
  TString intFlowVsMName = "fIntFlowVsM";
  intFlowVsMName += fAnalysisLabel->Data();
  for(Int_t co=0;co<4;co++) // cumulant order
  {
   fIntFlowVsM[co] = new TH1D(Form("%s, %s",intFlowVsMName.Data(),flowFlag[co].Data()),
                              Form("%s vs multiplicity",flowFlag[co].Data()),
                              fnBinsMult,fMinMult,fMaxMult);
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowVsM[co]->GetXaxis()->SetTitle("# RPs");
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowVsM[co]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowVsM[co]->GetXaxis()->SetTitle("# POIs"); 
       }
   fIntFlowVsM[co]->GetYaxis()->SetTitle(flowFlag[co].Data());  
   fIntFlowResults->Add(fIntFlowVsM[co]);                                    
  } // end of for(Int_t co=0;co<4;co++) // cumulant order
 } // end of if(fCalculateCumulantsVsM)
 // quantifying detector effects effects to correlations:
 TString intFlowDetectorBiasName = "fIntFlowDetectorBias";
 intFlowDetectorBiasName += fAnalysisLabel->Data();  
 fIntFlowDetectorBias = new TH1D(intFlowDetectorBiasName.Data(),"Quantifying detector bias",4,0,4);
 fIntFlowDetectorBias->SetLabelSize(0.05);
 fIntFlowDetectorBias->SetMarkerStyle(25);
 for(Int_t ci=0;ci<4;ci++)
 {  
  (fIntFlowDetectorBias->GetXaxis())->SetBinLabel(ci+1,Form("#frac{corrected}{measured} %s",cumulantFlag[ci].Data()));
 }
 fIntFlowResults->Add(fIntFlowDetectorBias); 
 // quantifying detector effects to correlations versus multiplicity:
 if(fCalculateCumulantsVsM)
 {
  TString intFlowDetectorBiasVsMName = "fIntFlowDetectorBiasVsM";
  intFlowDetectorBiasVsMName += fAnalysisLabel->Data();
  for(Int_t ci=0;ci<4;ci++) // correlation index
  {
   fIntFlowDetectorBiasVsM[ci] = new TH1D(Form("%s for %s",intFlowDetectorBiasVsMName.Data(),cumulantFlag[ci].Data()),
                                          Form("Quantifying detector bias for %s vs multiplicity",cumulantFlag[ci].Data()),
                                          fnBinsMult,fMinMult,fMaxMult);
   if(fMultiplicityIs==AliFlowCommonConstants::kRP)
   {
    fIntFlowDetectorBiasVsM[ci]->GetXaxis()->SetTitle("# RPs"); 
   } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
     {
      fIntFlowDetectorBiasVsM[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
     } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
       {
        fIntFlowDetectorBiasVsM[ci]->GetXaxis()->SetTitle("# POIs");
       }
   fIntFlowDetectorBiasVsM[ci]->GetYaxis()->SetTitle("#frac{corrected}{measured}");  
   fIntFlowResults->Add(fIntFlowDetectorBiasVsM[ci]);                                    
  } // end of for(Int_t co=0;co<4;co++) // cumulant order
 } // end of if(fCalculateCumulantsVsM)
   
} // end of AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::BookEverythingForControlHistograms()
{
 // Book all objects for control histograms.

 // a) Book profile to hold all flags for control histograms;
 // b) Book all control histograms.

 // a) Book profile to hold all flags for control histograms:
 TString controlHistogramsFlagsName = "fControlHistogramsFlags";
 controlHistogramsFlagsName += fAnalysisLabel->Data();
 fControlHistogramsFlags = new TProfile(controlHistogramsFlagsName.Data(),"Flags for Control Histograms",2,0,2);
 fControlHistogramsFlags->SetTickLength(-0.01,"Y");
 fControlHistogramsFlags->SetMarkerStyle(25);
 fControlHistogramsFlags->SetLabelSize(0.04);
 fControlHistogramsFlags->SetLabelOffset(0.02,"Y");
 fControlHistogramsFlags->SetStats(kFALSE);
 fControlHistogramsFlags->GetXaxis()->SetBinLabel(1,"fStoreControlHistograms");
 fControlHistogramsFlags->GetXaxis()->SetBinLabel(2,"fUseQvectorTerms");
 fControlHistogramsList->Add(fControlHistogramsFlags);

 if(!fStoreControlHistograms){return;}

 // b) Book all control histograms:
 //  b1) Correlation between # RPs and ref. mult. determined centrally:
 TString sCorrelationNoRPsVsRefMultName = "fCorrelationNoRPsVsRefMult";
 sCorrelationNoRPsVsRefMultName += fAnalysisLabel->Data();
 fCorrelationNoRPsVsRefMult = new TH2D(sCorrelationNoRPsVsRefMultName.Data(),"# RPs vs. Reference Multiplicity",fnBinsMult,fMinMult,fMaxMult,fnBinsMult,fMinMult,fMaxMult);
 fCorrelationNoRPsVsRefMult->SetTickLength(-0.01,"Y");
 fCorrelationNoRPsVsRefMult->SetLabelSize(0.04);
 fCorrelationNoRPsVsRefMult->SetLabelOffset(0.02,"Y");
 fCorrelationNoRPsVsRefMult->SetStats(kTRUE);
 fCorrelationNoRPsVsRefMult->GetXaxis()->SetTitle("# RPs");
 fCorrelationNoRPsVsRefMult->GetYaxis()->SetTitle("Reference Multiplicity");
 fControlHistogramsList->Add(fCorrelationNoRPsVsRefMult);
 //  b2) Correlation between # POIs and ref. mult. determined centrally:
 TString sCorrelationNoPOIsVsRefMultName = "fCorrelationNoPOIsVsRefMult";
 sCorrelationNoPOIsVsRefMultName += fAnalysisLabel->Data();
 fCorrelationNoPOIsVsRefMult = new TH2D(sCorrelationNoPOIsVsRefMultName.Data(),"# POIs vs. Reference Multiplicity",fnBinsMult,fMinMult,fMaxMult,fnBinsMult,fMinMult,fMaxMult);
 fCorrelationNoPOIsVsRefMult->SetTickLength(-0.01,"Y");
 fCorrelationNoPOIsVsRefMult->SetLabelSize(0.04);
 fCorrelationNoPOIsVsRefMult->SetLabelOffset(0.02,"Y");
 fCorrelationNoPOIsVsRefMult->SetStats(kTRUE);
 fCorrelationNoPOIsVsRefMult->GetXaxis()->SetTitle("# POIs");
 fCorrelationNoPOIsVsRefMult->GetYaxis()->SetTitle("Reference Multiplicity");
 fControlHistogramsList->Add(fCorrelationNoPOIsVsRefMult);
 //  b3) Correlation between # RPs and # POIs:
 TString sCorrelationNoRPsVsNoPOIsName = "fCorrelationNoRPsVsNoPOIs";
 sCorrelationNoRPsVsNoPOIsName += fAnalysisLabel->Data();
 fCorrelationNoRPsVsNoPOIs = new TH2D(sCorrelationNoRPsVsNoPOIsName.Data(),"# RPs vs. # POIs",fnBinsMult,fMinMult,fMaxMult,fnBinsMult,fMinMult,fMaxMult);
 fCorrelationNoRPsVsNoPOIs->SetTickLength(-0.01,"Y");
 fCorrelationNoRPsVsNoPOIs->SetLabelSize(0.04);
 fCorrelationNoRPsVsNoPOIs->SetLabelOffset(0.02,"Y");
 fCorrelationNoRPsVsNoPOIs->SetStats(kTRUE);
 fCorrelationNoRPsVsNoPOIs->GetXaxis()->SetTitle("# RPs");
 fCorrelationNoRPsVsNoPOIs->GetYaxis()->SetTitle("# POIs");
 fControlHistogramsList->Add(fCorrelationNoRPsVsNoPOIs);
 // b4) <2>, <4>, <6> and <8> vs multiplicity (#RPs, #POIs or external):
 TString sCorrelation[4] = {"#LT2#GT","#LT4#GT","#LT6#GT","#LT8#GT"};
 TString sMultiplicity = "";
 if(fMultiplicityIs==AliFlowCommonConstants::kRP)
 {
  sMultiplicity = "# RPs"; 
 } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
   {
    sMultiplicity = "Reference multiplicity (from ESD)";
   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
     {
      sMultiplicity = "# POIs";
     }
 for(Int_t ci=0;ci<4;ci++)
 {
  fCorrelation2468VsMult[ci] = new TH2D(Form("%s vs M",sCorrelation[ci].Data()),Form("%s vs M",sCorrelation[ci].Data()),fnBinsMult,fMinMult,fMaxMult,fnBinsForCorrelations,fMinValueOfCorrelation[ci],fMaxValueOfCorrelation[ci]); 
  fCorrelation2468VsMult[ci]->SetTickLength(-0.01,"Y");
  fCorrelation2468VsMult[ci]->SetLabelSize(0.04);
  fCorrelation2468VsMult[ci]->SetLabelOffset(0.02,"Y");
  fCorrelation2468VsMult[ci]->SetStats(kTRUE);
  fCorrelation2468VsMult[ci]->GetXaxis()->SetTitle(sMultiplicity.Data());
  fCorrelation2468VsMult[ci]->GetYaxis()->SetTitle(sCorrelation[ci].Data());
  fControlHistogramsList->Add(fCorrelation2468VsMult[ci]); 
 } // end of for(Int_t ci=0;ci<4;ci++)
 // b5) <2><4>, <2><6>, <2><8>, <4><6> etc. vs multiplicity (#RPs, #POIs or external):
 TString sCorrelationProduct[1] = {"#LT2#GT#LT4#GT"}; // TBI: add the other ones when needed first time
 for(Int_t cpi=0;cpi<1;cpi++) // TBI: hardwired 1
 {
  fCorrelationProduct2468VsMult[cpi] = new TH2D(Form("%s vs M",sCorrelationProduct[cpi].Data()),Form("%s vs M",sCorrelationProduct[cpi].Data()),fnBinsMult,fMinMult,fMaxMult,fnBinsForCorrelations,fMinValueOfCorrelationProduct[cpi],fMaxValueOfCorrelationProduct[cpi]);
  fCorrelationProduct2468VsMult[cpi]->SetTickLength(-0.01,"Y");
  fCorrelationProduct2468VsMult[cpi]->SetLabelSize(0.04);
  fCorrelationProduct2468VsMult[cpi]->SetLabelOffset(0.02,"Y");
  fCorrelationProduct2468VsMult[cpi]->SetStats(kTRUE);
  fCorrelationProduct2468VsMult[cpi]->GetXaxis()->SetTitle(sMultiplicity.Data());
  fCorrelationProduct2468VsMult[cpi]->GetYaxis()->SetTitle(sCorrelationProduct[cpi].Data());
  fControlHistogramsList->Add(fCorrelationProduct2468VsMult[cpi]); 
 } // end of for(Int_t cpi=0;cpi<4;cpi++)
 // b6) |Qn|^2/M, |Q2n|^2/M, |Qn|^4/(M(2M-1)), Re[Q2nQn^*Qn^*]/M, ... vs multiplicity (#RPs, #POIs or external)
 if(fUseQvectorTerms)
 {
  TString sQvectorTerms[4] = {"#frac{|Q_{n}|^{2}}{M}","#frac{|Q_{2n}|^{2}}{M}","#frac{|Q_{n}|^{4}}{M(2M-1)}","#frac{Re[Q_{2n}Q_{n}^{*}Q_{n}^{*}]}{M^{3/2}}"}; // TBI: add the other ones when needed first time
  for(Int_t qvti=0;qvti<4;qvti++) // TBI: hardwired 4
  {
   fQvectorTermsVsMult[qvti] = new TH2D(Form("%s vs M",sQvectorTerms[qvti].Data()),Form("%s vs M",sQvectorTerms[qvti].Data()),fnBinsMult,fMinMult,fMaxMult,fnBinsForCorrelations,fMinValueOfQvectorTerms[qvti],fMaxValueOfQvectorTerms[qvti]); 
   fQvectorTermsVsMult[qvti]->SetTickLength(-0.01,"Y");
   fQvectorTermsVsMult[qvti]->SetLabelSize(0.04);
   fQvectorTermsVsMult[qvti]->SetLabelOffset(0.02,"Y");
   fQvectorTermsVsMult[qvti]->SetStats(kTRUE);
   fQvectorTermsVsMult[qvti]->GetXaxis()->SetTitle(sMultiplicity.Data());
   fQvectorTermsVsMult[qvti]->GetYaxis()->SetTitle(sQvectorTerms[qvti].Data());
   fControlHistogramsList->Add(fQvectorTermsVsMult[qvti]); 
  } // end of for(Int_t qvti=0;qvti<4;qvti++)
 } // end of if(fUseQvectorTerms)

} // end of void AliFlowAnalysisWithQCumulants::BookEverythingForControlHistograms()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::BookEverythingForBootstrap()
{
 // Book all objects needed for bootstrap. 

 // a) Book profile to hold all flags for bootstrap;
 // b) Book local random generator;
 // c) Book all bootstrap objects;
 // d) Book all bootstrap objects 'vs M'.

 // a) Book profile to hold all flags for bootstrap;
 TString bootstrapFlagsName = "fBootstrapFlags";
 bootstrapFlagsName += fAnalysisLabel->Data();
 fBootstrapFlags = new TProfile(bootstrapFlagsName.Data(),"Flags for bootstrap",3,0,3);
 fBootstrapFlags->SetTickLength(-0.01,"Y");
 fBootstrapFlags->SetMarkerStyle(25);
 fBootstrapFlags->SetLabelSize(0.04);
 fBootstrapFlags->SetLabelOffset(0.02,"Y");
 fBootstrapFlags->SetStats(kFALSE);
 fBootstrapFlags->GetXaxis()->SetBinLabel(1,"fUseBootstrap");
 fBootstrapFlags->GetXaxis()->SetBinLabel(2,"fUseBootstrapVsM");
 fBootstrapFlags->GetXaxis()->SetBinLabel(3,"fnSubsamples");
 fBootstrapList->Add(fBootstrapFlags);

 // b) Book local random generator:
 if(fUseBootstrap||fUseBootstrapVsM)
 { 
  fRandom = new TRandom3(0); // if uiSeed is 0, the seed is determined uniquely in space and time via TUUID
 }

 // c) Book all bootstrap objects:
 TString correlationFlag[4] = {"#LT#LT2#GT#GT","#LT#LT4#GT#GT","#LT#LT6#GT#GT","#LT#LT8#GT#GT"};
 TString cumulantFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
 if(fUseBootstrap)
 {
  // ....
  TString bootstrapCorrelationsName = "fBootstrapCorrelations";
  bootstrapCorrelationsName += fAnalysisLabel->Data();
  fBootstrapCorrelations = new TProfile2D(bootstrapCorrelationsName.Data(),"Bootstrap Correlations",4,0.,4.,fnSubsamples,0,fnSubsamples); // x-axis => <2>, <4>, <6>, <8>; y-axis => subsample # 
  fBootstrapCorrelations->SetStats(kFALSE);
  for(Int_t ci=0;ci<4;ci++) // correlation index
  {
   fBootstrapCorrelations->GetXaxis()->SetBinLabel(ci+1,correlationFlag[ci].Data());
  } // end of for(Int_t ci=0;ci<4;ci++) // correlation index
  for(Int_t ss=0;ss<fnSubsamples;ss++)
  {
   fBootstrapCorrelations->GetYaxis()->SetBinLabel(ss+1,Form("#%d",ss));
  } // end of for(Int_t ss=0;ss<fnSubsamples;ss++)
  fBootstrapProfilesList->Add(fBootstrapCorrelations);
  // ....
  TString bootstrapCumulantsName = "fBootstrapCumulants";
  bootstrapCumulantsName += fAnalysisLabel->Data();
  fBootstrapCumulants = new TH2D(bootstrapCumulantsName.Data(),"Bootstrap Cumulants",4,0.,4.,fnSubsamples,0,fnSubsamples); // x-axis => QC{2}, QC{4}, QC{6}, QC{8}; y-axis => subsample # 
  fBootstrapCumulants->SetStats(kFALSE);
  for(Int_t co=0;co<4;co++) // cumulant order
  {
   fBootstrapCumulants->GetXaxis()->SetBinLabel(co+1,cumulantFlag[co].Data());
  } // end of for(Int_t co=0;co<4;co++) // cumulant order
  for(Int_t ss=0;ss<fnSubsamples;ss++)
  {
   fBootstrapCumulants->GetYaxis()->SetBinLabel(ss+1,Form("#%d",ss));
  } // end of for(Int_t ss=0;ss<fnSubsamples;ss++)
  fBootstrapResultsList->Add(fBootstrapCumulants);
 } // end of if(fUseBootstrap)

 // d) Book all bootstrap objects 'vs M':
 TString sMultiplicity = "";
 if(fMultiplicityIs==AliFlowCommonConstants::kRP)
 {
  sMultiplicity = "# RPs"; 
 } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
   {
    sMultiplicity = "Reference multiplicity (from ESD)";
   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
     {
      sMultiplicity = "# POIs";
     }
 if(fUseBootstrapVsM)
 {
  // ....
  TString bootstrapCorrelationsVsMName = "fBootstrapCorrelationsVsM";
  bootstrapCorrelationsVsMName += fAnalysisLabel->Data();
  for(Int_t ci=0;ci<4;ci++) // correlation index
  {
   fBootstrapCorrelationsVsM[ci] = new TProfile2D(Form("%s, %s",bootstrapCorrelationsVsMName.Data(),correlationFlag[ci].Data()),
                                       Form("Bootstrap Correlations Vs. M, %s",correlationFlag[ci].Data()),
                                       fnBinsMult,fMinMult,fMaxMult,fnSubsamples,0,fnSubsamples); // index => <2>, <4>, <6>, <8>; x-axis => multiplicity; y-axis => subsample # 
   fBootstrapCorrelationsVsM[ci]->SetStats(kFALSE);
   fBootstrapCorrelationsVsM[ci]->GetXaxis()->SetTitle(sMultiplicity.Data());
   for(Int_t ss=0;ss<fnSubsamples;ss++)
   {
    fBootstrapCorrelationsVsM[ci]->GetYaxis()->SetBinLabel(ss+1,Form("#%d",ss));
   } // end of for(Int_t ss=0;ss<fnSubsamples;ss++)
   fBootstrapProfilesList->Add(fBootstrapCorrelationsVsM[ci]);
  } // end of for(Int_t ci=0;ci<4;ci++) // correlation index 
  // ....
  TString bootstrapCumulantsVsMName = "fBootstrapCumulantsVsM";
  bootstrapCumulantsVsMName += fAnalysisLabel->Data();
  for(Int_t co=0;co<4;co++) // cumulant order
  {
   fBootstrapCumulantsVsM[co] = new TH2D(Form("%s, %s",bootstrapCumulantsVsMName.Data(),cumulantFlag[co].Data()),
                                       Form("Bootstrap Cumulants Vs. M, %s",cumulantFlag[co].Data()),
                                       fnBinsMult,fMinMult,fMaxMult,fnSubsamples,0,fnSubsamples); // index => <2>, <4>, <6>, <8>; x-axis => multiplicity; y-axis => subsample # 
   fBootstrapCumulantsVsM[co]->SetStats(kFALSE);
   fBootstrapCumulantsVsM[co]->GetXaxis()->SetTitle(sMultiplicity.Data());
   for(Int_t ss=0;ss<fnSubsamples;ss++)
   {
    fBootstrapCumulantsVsM[co]->GetYaxis()->SetBinLabel(ss+1,Form("#%d",ss));
   } // end of for(Int_t ss=0;ss<fnSubsamples;ss++)
   fBootstrapResultsList->Add(fBootstrapCumulantsVsM[co]);
  } // end of for(Int_t co=0;co<4;co++) // correlation index 
 } // end of if(fUseBootstrapVsM)

} // end of void AliFlowAnalysisWithQCumulants::BookEverythingForBootstrap()

//=======================================================================================================================

void AliFlowAnalysisWithQCumulants::BookEverythingForMixedHarmonics()
{
 // Book all objects for mixed harmonics.

 // a) Book profile to hold all flags for mixed harmonics;
 // b) Book all objects in TList fMixedHarmonicsProfiles;
 // c) Book all objects in TList fMixedHarmonicsResults;
 // d) Book all objects in TList fMixedHarmonicsErrorPropagation.

 // a) Book profile to hold all flags for mixed harmonics:
 TString mixedHarmonicsFlagsName = "fMixedHarmonicsFlags";
 mixedHarmonicsFlagsName += fAnalysisLabel->Data();
 fMixedHarmonicsFlags = new TProfile(mixedHarmonicsFlagsName.Data(),"Flags for Mixed Harmonics",4,0,4);
 fMixedHarmonicsFlags->SetTickLength(-0.01,"Y");
 fMixedHarmonicsFlags->SetMarkerStyle(25);
 fMixedHarmonicsFlags->SetLabelSize(0.04);
 fMixedHarmonicsFlags->SetLabelOffset(0.02,"Y");
 fMixedHarmonicsFlags->SetStats(kFALSE);
 fMixedHarmonicsFlags->GetXaxis()->SetBinLabel(1,"Calculate Mixed Harmonics");
 fMixedHarmonicsFlags->GetXaxis()->SetBinLabel(2,"Generic Harmonic");
 fMixedHarmonicsFlags->GetXaxis()->SetBinLabel(3,"Calculate vs Multiplicity");
 fMixedHarmonicsFlags->GetXaxis()->SetBinLabel(4,"Multiplicity Weight");
 fMixedHarmonicsList->Add(fMixedHarmonicsFlags);

 if(!fCalculateMixedHarmonics){return;}

 // b) Book all objects in TList fMixedHarmonicsProfiles:
 //  b1) 2-p correlations:
 TString s2pCorrelationsName = "f2pCorrelations";
 s2pCorrelationsName += fAnalysisLabel->Data();
 f2pCorrelations = new TProfile(s2pCorrelationsName.Data(),Form("2-particle correlations (n = %d)",fHarmonic),6,0,6,"s");
 f2pCorrelations->SetTickLength(-0.01,"Y");
 f2pCorrelations->SetMarkerStyle(25);
 f2pCorrelations->SetLabelSize(0.04);
 f2pCorrelations->SetLabelOffset(0.02,"Y");
 f2pCorrelations->SetStats(kFALSE);
 f2pCorrelations->Sumw2();
 f2pCorrelations->GetXaxis()->SetBinLabel(1,Form("#LT#LT2#GT#GT_{%dn|%dn}",1*fHarmonic,1*fHarmonic));
 f2pCorrelations->GetXaxis()->SetBinLabel(2,Form("#LT#LT2#GT#GT_{%dn|%dn}",2*fHarmonic,2*fHarmonic));
 f2pCorrelations->GetXaxis()->SetBinLabel(3,Form("#LT#LT2#GT#GT_{%dn|%dn}",3*fHarmonic,3*fHarmonic));
 f2pCorrelations->GetXaxis()->SetBinLabel(4,Form("#LT#LT2#GT#GT_{%dn|%dn}",4*fHarmonic,4*fHarmonic));
 f2pCorrelations->GetXaxis()->SetBinLabel(5,Form("#LT#LT2#GT#GT_{%dn|%dn}",5*fHarmonic,5*fHarmonic));
 f2pCorrelations->GetXaxis()->SetBinLabel(6,Form("#LT#LT2#GT#GT_{%dn|%dn}",6*fHarmonic,6*fHarmonic));
 fMixedHarmonicsProfiles->Add(f2pCorrelations);
 //  b2) 3-p correlations (3+6):
 TString s3pCorrelationsName = "f3pCorrelations";
 s3pCorrelationsName += fAnalysisLabel->Data();
 f3pCorrelations = new TProfile(s3pCorrelationsName.Data(),Form("3-particle correlations (n = %d)",fHarmonic),10,0,10,"s"); 
 f3pCorrelations->SetTickLength(-0.01,"Y");
 f3pCorrelations->SetMarkerStyle(25);
 f3pCorrelations->SetLabelSize(0.04);
 f3pCorrelations->SetLabelOffset(0.02,"Y");
 f3pCorrelations->SetStats(kFALSE);
 f3pCorrelations->Sumw2();
 // 3-p correlations sensitive to two distinct harmonics (3):
 f3pCorrelations->GetXaxis()->SetBinLabel(1,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",2*fHarmonic,1*fHarmonic,1*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(2,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,2*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(3,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,3*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(4,""); // empty
 // 3-p correlations sensitive to three distinct harmonics (6):
 f3pCorrelations->GetXaxis()->SetBinLabel(5,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",3*fHarmonic,2*fHarmonic,1*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(6,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",4*fHarmonic,3*fHarmonic,1*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(7,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,2*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(8,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",5*fHarmonic,4*fHarmonic,1*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(9,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,2*fHarmonic));
 f3pCorrelations->GetXaxis()->SetBinLabel(10,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicsProfiles->Add(f3pCorrelations);
 //  b3) 4-p correlations (6+15+2+10+8):
 TString s4pCorrelationsName = "f4pCorrelations";
 s4pCorrelationsName += fAnalysisLabel->Data();
 f4pCorrelations = new TProfile(s4pCorrelationsName.Data(),Form("4-particle correlations (n = %d)",fHarmonic),45,0,45,"s");
 f4pCorrelations->SetTickLength(-0.01,"Y");
 f4pCorrelations->SetMarkerStyle(25);
 f4pCorrelations->SetLabelSize(0.03);
 f4pCorrelations->SetLabelOffset(0.02,"Y");
 f4pCorrelations->SetStats(kFALSE);
 f4pCorrelations->Sumw2();
 // "same harmonic" (6):
 f4pCorrelations->GetXaxis()->SetBinLabel(1,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(2,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",2*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(3,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(4,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(5,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,5*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(6,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,6*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(7,""); // empty
 // "standard candles" (15):
 f4pCorrelations->GetXaxis()->SetBinLabel(8,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",2*fHarmonic,1*fHarmonic,2*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(9,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,1*fHarmonic,3*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(10,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(11,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,4*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(12,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(13,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(14,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,5*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(15,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,5*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(16,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,5*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(17,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,4*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(18,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,6*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(19,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,6*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(20,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,6*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(21,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,4*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(22,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,5*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(23,""); // empty
 // 4-p correlations sensitive to two distinct harmonics (2):
 f4pCorrelations->GetXaxis()->SetBinLabel(24,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",3*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(25,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(26,""); // empty
 // 4-p correlations sensitive to three distinct harmonics (10):
 f4pCorrelations->GetXaxis()->SetBinLabel(27,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(28,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(29,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(30,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(31,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(32,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(33,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,4*fHarmonic,4*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(34,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,1*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(35,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,4*fHarmonic,4*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(36,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,5*fHarmonic,5*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(37,""); // empty
 // 4-p correlations sensitive to four distinct harmonics (8):
 f4pCorrelations->GetXaxis()->SetBinLabel(38,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,3*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(39,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,4*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(40,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,4*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(41,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,4*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(42,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,5*fHarmonic,2*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(43,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(44,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,5*fHarmonic,3*fHarmonic));
 f4pCorrelations->GetXaxis()->SetBinLabel(45,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,5*fHarmonic,4*fHarmonic));
 fMixedHarmonicsProfiles->Add(f4pCorrelations);
 //  b3) 5-p correlations (30+9+30+11+3):
 TString s5pCorrelationsName = "f5pCorrelations";
 s5pCorrelationsName += fAnalysisLabel->Data();
 f5pCorrelations = new TProfile(s5pCorrelationsName.Data(),Form("5-particle correlations (n = %d)",fHarmonic),87,0,87,"s");
 f5pCorrelations->SetTickLength(-0.01,"Y");
 f5pCorrelations->SetMarkerStyle(25);
 f5pCorrelations->SetLabelSize(0.02);
 f5pCorrelations->SetLabelOffset(0.02,"Y");
 f5pCorrelations->SetStats(kFALSE);
 f5pCorrelations->Sumw2();
 // "standard candles" (30):
 f5pCorrelations->GetXaxis()->SetBinLabel(1,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,2*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(2,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(3,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(4,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(5,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,4*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(6,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,3*fHarmonic,4*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(7,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,1*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(8,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,5*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(9,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(10,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(11,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(12,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,5*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(13,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(14,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(15,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(16,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(17,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,4*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(18,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(19,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(20,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,5*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(21,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,5*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(22,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,5*fHarmonic,4*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(23,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,5*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(24,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,6*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(25,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,6*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(26,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(27,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(28,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,5*fHarmonic,4*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(29,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(30,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,4*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(31,""); // empty
 // 5-p correlations sensitive to two distinct harmonics (9):
 f5pCorrelations->GetXaxis()->SetBinLabel(32,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",2*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(33,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",2*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(34,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(35,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",4*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(36,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(37,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(38,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(39,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(40,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(41,""); // empty
 // 5-p correlations sensitive to three distinct harmonics (30):
 f5pCorrelations->GetXaxis()->SetBinLabel(42,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,1*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(43,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(44,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(45,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,1*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(46,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(47,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(48,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(49,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(50,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(51,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(52,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(53,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(54,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,1*fHarmonic,4*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(55,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(56,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(57,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(58,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,4*fHarmonic,4*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(59,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,3*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(60,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,4*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(61,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(62,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(63,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,1*fHarmonic,4*fHarmonic,4*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(64,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,5*fHarmonic,1*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(65,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(66,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(67,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,5*fHarmonic,5*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(68,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,5*fHarmonic,5*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(69,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,5*fHarmonic,5*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(70,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,4*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(71,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,5*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(72,""); // empty
 // 5-p correlations sensitive to four distinct harmonics (11):
 f5pCorrelations->GetXaxis()->SetBinLabel(73,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(74,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,1*fHarmonic,4*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(75,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,4*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(76,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,1*fHarmonic,4*fHarmonic,4*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(77,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(78,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,4*fHarmonic,4*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(79,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,1*fHarmonic,5*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(80,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,5*fHarmonic,2*fHarmonic,2*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(81,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(82,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,1*fHarmonic,5*fHarmonic,5*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(83,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,5*fHarmonic,4*fHarmonic,3*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(84,""); // empty
 // 5-p correlations sensitive to five distinct harmonics (3):
 f5pCorrelations->GetXaxis()->SetBinLabel(85,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,4*fHarmonic,3*fHarmonic,1*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(86,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,1*fHarmonic,5*fHarmonic,4*fHarmonic));
 f5pCorrelations->GetXaxis()->SetBinLabel(87,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,5*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicsProfiles->Add(f5pCorrelations);
 //  b4) 6-p correlations (??+??+??+??+??):
 TString s6pCorrelationsName = "f6pCorrelations";
 s6pCorrelationsName += fAnalysisLabel->Data();
 f6pCorrelations = new TProfile(s6pCorrelationsName.Data(),Form("6-particle correlations (n = %d)",fHarmonic),1,0.,1.);
 f6pCorrelations->SetTickLength(-0.01,"Y");
 f6pCorrelations->SetMarkerStyle(25);
 f6pCorrelations->SetLabelSize(0.02);
 f6pCorrelations->SetLabelOffset(0.02,"Y");
 f6pCorrelations->SetStats(kFALSE);
 f6pCorrelations->Sumw2(); 
 //fMixedHarmonicsProfiles->Add(f6pCorrelations); // TBI
 //  b5) 7-p correlations (??+??+??+??+??):
 TString s7pCorrelationsName = "f7pCorrelations";
 s7pCorrelationsName += fAnalysisLabel->Data();
 f7pCorrelations = new TProfile(s7pCorrelationsName.Data(),Form("7-particle correlations (n = %d)",fHarmonic),1,0.,1.);
 f7pCorrelations->SetTickLength(-0.01,"Y");
 f7pCorrelations->SetMarkerStyle(25);
 f7pCorrelations->SetLabelSize(0.02);
 f7pCorrelations->SetLabelOffset(0.02,"Y");
 f7pCorrelations->SetStats(kFALSE);
 f7pCorrelations->Sumw2(); 
 //fMixedHarmonicsProfiles->Add(f7pCorrelations); // TBI
 //  b6) 8-p correlations (??+??+??+??+??):
 TString s8pCorrelationsName = "f8pCorrelations";
 s8pCorrelationsName += fAnalysisLabel->Data();
 f8pCorrelations = new TProfile(s8pCorrelationsName.Data(),Form("8-particle correlations (n = %d)",fHarmonic),1,0.,1.);
 f8pCorrelations->SetTickLength(-0.01,"Y");
 f8pCorrelations->SetMarkerStyle(25);
 f8pCorrelations->SetLabelSize(0.02);
 f8pCorrelations->SetLabelOffset(0.02,"Y");
 f8pCorrelations->SetStats(kFALSE);
 f8pCorrelations->Sumw2(); 
 //fMixedHarmonicsProfiles->Add(f8pCorrelations); // TBI

 // c) Book all objects in TList fMixedHarmonicsResults:
 // QC{2}:
 f2pCumulants = f2pCorrelations->ProjectionX("f2pCumulants");
 f2pCumulants->SetTitle(Form("2-particle cumulants (n = %d)",fHarmonic));
 f2pCumulants->SetStats(kFALSE);
 f2pCumulants->SetMarkerStyle(kFullSquare);
 f2pCumulants->SetMarkerColor(kBlack);
 f2pCumulants->SetLineColor(kBlack);
 fMixedHarmonicsResults->Add(f2pCumulants);
 // QC{3}:
 f3pCumulants = f3pCorrelations->ProjectionX("f3pCumulants");
 f3pCumulants->SetTitle(Form("3-particle cumulants (n = %d)",fHarmonic));
 f3pCumulants->SetStats(kFALSE);
 f3pCumulants->SetMarkerStyle(kFullSquare);
 f3pCumulants->SetMarkerColor(kGreen+2);
 f3pCumulants->SetLineColor(kGreen+2);
 fMixedHarmonicsResults->Add(f3pCumulants);
 // QC{4}:
 f4pCumulants = f4pCorrelations->ProjectionX("f4pCumulants");
 f4pCumulants->SetTitle(Form("4-particle cumulants (n = %d)",fHarmonic));
 f4pCumulants->SetStats(kFALSE);
 f4pCumulants->SetMarkerStyle(kFullSquare);
 f4pCumulants->SetMarkerColor(kRed);
 f4pCumulants->SetLineColor(kRed);
 fMixedHarmonicsResults->Add(f4pCumulants);
 // QC{5}:
 f5pCumulants = f5pCorrelations->ProjectionX("f5pCumulants");
 f5pCumulants->SetTitle(Form("5-particle cumulants (n = %d)",fHarmonic));
 f5pCumulants->SetStats(kFALSE);
 f5pCumulants->SetMarkerStyle(kFullSquare);
 f5pCumulants->SetMarkerColor(kBlue);
 f5pCumulants->SetLineColor(kBlue);
 fMixedHarmonicsResults->Add(f5pCumulants);

 // d) Book all objects in TList fMixedHarmonicsErrorPropagation: 
 // Sum of linear and quadratic event weights for mixed harmonics => [0=linear 1,1=quadratic]: 
 TString mixedHarmonicEventWeightsName = "fMixedHarmonicEventWeights";
 mixedHarmonicEventWeightsName += fAnalysisLabel->Data();
 TString powerFlag[2] = {"linear","quadratic"}; 
 for(Int_t power=0;power<2;power++)
 {
  fMixedHarmonicEventWeights[power] = new TH1D(Form("%s: %s",mixedHarmonicEventWeightsName.Data(),powerFlag[power].Data()),Form("Sum of %s event weights for correlations",powerFlag[power].Data()),8,0.,8.);
  fMixedHarmonicEventWeights[power]->SetLabelSize(0.04);
  fMixedHarmonicEventWeights[power]->SetMarkerStyle(25);
  fMixedHarmonicEventWeights[power]->SetStats(kFALSE);
  if(power == 0)
  {
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum w_{#LT1#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum w_{#LT2#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum w_{#LT3#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum w_{#LT4#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(5,"#sum w_{#LT5#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(6,"#sum w_{#LT6#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(7,"#sum w_{#LT7#GT}");
   (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(8,"#sum w_{#LT8#GT}");
  } else if (power == 1) 
    {
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum w_{#LT1#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum w_{#LT2#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum w_{#LT3#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum w_{#LT4#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(5,"#sum w_{#LT5#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(6,"#sum w_{#LT6#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(7,"#sum w_{#LT7#GT}^{2}");
     (fMixedHarmonicEventWeights[power]->GetXaxis())->SetBinLabel(8,"#sum w_{#LT8#GT}^{2}");
    }
  fMixedHarmonicsErrorPropagation->Add(fMixedHarmonicEventWeights[power]);
 } // end of for(Int_t power=0;power<2;power++)
 
 // Sums of products of event weights for mixed harmonics:
 TString mixedHarmonicProductOfEventWeightsName = "fMixedHarmonicProductOfEventWeights";
 mixedHarmonicProductOfEventWeightsName += fAnalysisLabel->Data();
 fMixedHarmonicProductOfEventWeights = new TH2D(mixedHarmonicProductOfEventWeightsName.Data(),"Sums of products of event weights",8,0.,8.,8,0.,8.);
 fMixedHarmonicProductOfEventWeights->SetStats(kFALSE);
 fMixedHarmonicProductOfEventWeights->GetXaxis()->SetLabelSize(0.05);
 fMixedHarmonicProductOfEventWeights->GetYaxis()->SetLabelSize(0.05);
 for(Int_t b=1;b<=8;b++)
 {
  fMixedHarmonicProductOfEventWeights->GetXaxis()->SetBinLabel(b,Form("w_{#LT%i#GT}",b)); 
  fMixedHarmonicProductOfEventWeights->GetYaxis()->SetBinLabel(b,Form("w_{#LT%i#GT}",b)); 
 } 
 fMixedHarmonicsErrorPropagation->Add(fMixedHarmonicProductOfEventWeights);

 // Averages of products of mixed harmonics correlations:
 TString mixedHarmonicProductOfCorrelationsName = "fMixedHarmonicProductOfCorrelations";
 mixedHarmonicProductOfCorrelationsName += fAnalysisLabel->Data();
 fMixedHarmonicProductOfCorrelations = new TProfile2D(mixedHarmonicProductOfCorrelationsName.Data(),"Averages of products of mixed correlators",139,0.,139.,139,0.,139.);
 fMixedHarmonicProductOfCorrelations->Sumw2();
 fMixedHarmonicProductOfCorrelations->SetStats(kFALSE);
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetLabelSize(0.015);
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetLabelSize(0.015);
 // x-axis:
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(1,Form("#LT#LT2#GT#GT_{%dn|%dn}",1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(2,Form("#LT#LT2#GT#GT_{%dn|%dn}",2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(3,Form("#LT#LT2#GT#GT_{%dn|%dn}",3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(4,Form("#LT#LT2#GT#GT_{%dn|%dn}",4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(5,Form("#LT#LT2#GT#GT_{%dn|%dn}",5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(6,Form("#LT#LT2#GT#GT_{%dn|%dn}",6*fHarmonic,6*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(7,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(8,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(9,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(10,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",3*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(11,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",4*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(12,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(13,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",5*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(14,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(15,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(16,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(17,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",2*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(18,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(19,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(20,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(21,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,6*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(22,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",2*fHarmonic,1*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(23,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,1*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(24,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(25,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(26,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(27,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(28,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(29,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,5*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(30,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,5*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(31,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(32,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,6*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(33,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,6*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(34,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,6*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(35,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(36,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(37,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",3*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(38,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(39,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(40,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(41,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(42,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(43,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(44,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(45,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(46,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(47,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(48,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(49,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(50,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(51,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(52,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(53,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,5*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(54,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(55,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,5*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(56,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,5*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(57,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,2*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(58,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(59,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(60,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(61,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,4*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(62,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,3*fHarmonic,4*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(63,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,1*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(64,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,5*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(65,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(66,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(67,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(68,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,5*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(69,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(70,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(71,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(72,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(73,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,4*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(74,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(75,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(76,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,5*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(77,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,5*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(78,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,5*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(79,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,5*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(80,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,6*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(81,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,6*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(82,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(83,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(84,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,5*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(85,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(86,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(87,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",2*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(88,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",2*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(89,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(90,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",4*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(91,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(92,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(93,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(94,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(95,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(96,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,1*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(97,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(98,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",3*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(99,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,1*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(100,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(101,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(102,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(103,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(104,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(105,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(106,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(107,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(108,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,1*fHarmonic,4*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(109,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(110,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(111,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,4*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(112,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,4*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(113,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(114,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(115,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(116,Form("#LT#LT5#GT#GT_{%dn|%dn,%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(117,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,1*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(118,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,5*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(119,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(120,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(121,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(122,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,5*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(123,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,5*fHarmonic,5*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(124,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(125,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(126,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(127,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,1*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(128,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,4*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(129,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,1*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(130,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,1*fHarmonic,3*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(131,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,4*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(132,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,1*fHarmonic,5*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(133,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,3*fHarmonic,5*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(134,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,5*fHarmonic,4*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(135,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,1*fHarmonic,5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(136,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,6*fHarmonic,5*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(137,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,4*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(138,Form("#LT#LT5#GT#GT_{%dn,%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,1*fHarmonic,5*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetXaxis()->SetBinLabel(139,Form("#LT#LT5#GT#GT_{%dn,%dn|%dn,%dn,%dn}",6*fHarmonic,4*fHarmonic,5*fHarmonic,3*fHarmonic,2*fHarmonic));
 // y-axis:
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(1,Form("#LT#LT2#GT#GT_{%dn|%dn}",1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(2,Form("#LT#LT2#GT#GT_{%dn|%dn}",2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(3,Form("#LT#LT2#GT#GT_{%dn|%dn}",3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(4,Form("#LT#LT2#GT#GT_{%dn|%dn}",4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(5,Form("#LT#LT2#GT#GT_{%dn|%dn}",5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(6,Form("#LT#LT2#GT#GT_{%dn|%dn}",6*fHarmonic,6*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(7,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(8,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(9,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(10,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",3*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(11,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",4*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(12,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(13,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",5*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(14,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(15,Form("#LT#LT3#GT#GT_{%dn|%dn,%dn}",6*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(16,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",1*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(17,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",2*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(18,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,3*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(19,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,4*fHarmonic,4*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(20,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,5*fHarmonic,5*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(21,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,6*fHarmonic,6*fHarmonic,6*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(22,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",2*fHarmonic,1*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(23,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,1*fHarmonic,3*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(24,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,2*fHarmonic,3*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(25,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,1*fHarmonic,4*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(26,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,4*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(27,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,3*fHarmonic,4*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(28,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,5*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(29,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,2*fHarmonic,5*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(30,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,5*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(31,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,4*fHarmonic,5*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(32,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,1*fHarmonic,6*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(33,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,2*fHarmonic,6*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(34,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,3*fHarmonic,6*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(35,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,4*fHarmonic,6*fHarmonic,4*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(36,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",6*fHarmonic,5*fHarmonic,6*fHarmonic,5*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(37,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",3*fHarmonic,1*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(38,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",6*fHarmonic,2*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(39,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",3*fHarmonic,1*fHarmonic,2*fHarmonic,2*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(40,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",4*fHarmonic,2*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(41,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",4*fHarmonic,2*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(42,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",5*fHarmonic,2*fHarmonic,2*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(43,Form("#LT#LT4#GT#GT_{%dn|%dn,%dn,%dn}",5*fHarmonic,3*fHarmonic,1*fHarmonic,1*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(44,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,1*fHarmonic,3*fHarmonic,3*fHarmonic));
 fMixedHarmonicProductOfCorrelations->GetYaxis()->SetBinLabel(45,Form("#LT#LT4#GT#GT_{%dn,%dn|%dn,%dn}",5*fHarmonic,3*fHarmonic,4*