AliFlowAnalysisCRC.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.                  *
 **************************************************************************/

/*************************************************
 * Charge-rapidity correlations with Q-cumulants *
 *                                               *
 * author: Jacopo Margutti                       *
 *         (margutti@nikhef.nl)                  *
 *************************************************/

#define AliFlowAnalysisCRC_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 "TProfile3D.h"
#include "TMath.h"
#include "TArrow.h"
#include "TPaveLabel.h"
#include "TCanvas.h"
#include "AliFlowEventSimple.h"
#include "AliFlowVector.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowAnalysisCRC.h"
#include "AliLog.h"
#include "TRandom.h"
#include "TF1.h"
#include "TNtuple.h"
#include "THnSparse.h"
#include <complex>
#include <cmath>

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(AliFlowAnalysisCRC)

AliFlowAnalysisCRC::AliFlowAnalysisCRC(const char* name,
                                           Int_t nCen,
                                           Double_t CenWidth):
TNamed(name,name),
// 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(1),
fAnalysisLabel(NULL),
// 2a.) particle weights:
fUsePhiWeights(kFALSE),
fUsePtWeights(kFALSE),
fUseEtaWeights(kFALSE),
fUseTrackWeights(kFALSE),
fUsePhiEtaWeights(kFALSE),
fUsePhiEtaWeightsChDep(kFALSE),
fUsePhiEtaWeightsVtxDep(kFALSE),
fUsePhiEtaWeightsChPtDep(kFALSE),
fUseZDCESEMulWeights(kFALSE),
fUseZDCESESpecWeights(kFALSE),
fUseParticleWeights(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),
fStoreVarious(kFALSE),
fExactNoRPs(0),
fUse2DHistograms(kFALSE),
fFillProfilesVsMUsingWeights(kTRUE),
fUseQvectorTerms(kFALSE),
fReQ(NULL),
fImQ(NULL),
fSpk(NULL),
fReQGF(NULL),
fImQGF(NULL),
fIntFlowCorrelationsEBE(NULL),
fIntFlowEventWeightsForCorrelationsEBE(NULL),
fIntFlowCorrelationsAllEBE(NULL),
fNumberOfRPsEBE(0.),
fNumberOfPOIsEBE(0.),
fReferenceMultiplicityEBE(0.),
fReferenceMultiplicityRecEBE(0.),
fCentralityEBE(0.),
fNewCentralityEBE(0.),
fZDCESEclEbE(0),
fNewMetricLEBE(0.),
fNewMetricDEBE(0.),
fNewMetricL2EBE(0.),
fNewMetricD2EBE(0.),
fCentralityCL1EBE(0.),
fNITSCL1EBE(0.),
fCentralityTRKEBE(0.),
fZNCQ0(0.),
fZNAQ0(0.),
fZNCen(0.),
fZNAen(0.),
fEnNucl(1.),
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),
// 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),
// 12.) Charge-Eta Asymmetry:
fCRCList(NULL),
fTempList(NULL),
fCRCFlags(NULL),
fCalculateCRC(kTRUE),
fCalculateCRCPt(kFALSE),
fCalculateCME(kFALSE),
fCalculateCRCInt(kFALSE),
fCalculateCRC2(kFALSE),
fCalculateCRCVZ(kFALSE),
fCalculateCRCZDC(kFALSE),
fCalculateFlowQC(kFALSE),
fCalculateFlowZDC(kFALSE),
fCalculateFlowVZ(kFALSE),
fCalculateEbEFlow(kFALSE),
fStoreZDCQVecVtxPos(kFALSE),
fUsePhiEtaCuts(kFALSE),
fUseVZERO(kFALSE),
fUseZDC(kFALSE),
fRecenterZDC(kFALSE),
fNUAforCRC(kFALSE),
fUseCRCRecenter(kFALSE),
fDivSigma(kTRUE),
fInvertZDC(kFALSE),
fCRCTestSin(kFALSE),
fVtxRbR(kFALSE),
fCRCEtaMin(0.),
fCRCEtaMax(0.),
fRunNum(0),
fCachedRunNum(0),
fRunBin(0),
fCenBin(0),
fCorrWeightTPC(kUnit),
fCorrWeightVZ(kUnit),
fCorrWeightZDC(kUnit),
fCRCIntList(NULL),
fCRCnCen(nCen),
fCRCCenBinWidth(CenWidth),
fCRCIntRbRList(NULL),
fCRCnRun(211),
fDataSet(kAny),
fInteractionRate(kAll),
fSelectCharge(kAllCh),
fPOIExtraWeights(kNone),
fCRCQVecList(NULL),
fCRCQVecListTPC(NULL),
fCRCQVecWeightsList(NULL),
fCRCZDCCalibList(NULL),
fCRCZDC2DCutList(NULL),
fCRCVZEROCalibList(NULL),
fCRCZDCResList(NULL),
fZDCESEList(NULL),
fCRCVZList(NULL),
fCRCZDCList(NULL),
fCRCZDCRbRList(NULL),
fCRCPtList(NULL),
fCMEList(NULL),
fCMETPCList(NULL),
fCMEZDCList(NULL),
fCRC2List(NULL),
fCRC2nEtaBins(6),
fFlowSPZDCList(NULL),
fFlowQCList(NULL),
//fFlowQCOrdMagList(NULL),
fFlowQCHOList(NULL),
fFlowGFList(NULL),
fFlowQCCenBin(100),
fFlowQCDeltaEta(0.4),
fFlowSPVZList(NULL),
fVariousList(NULL),
fEbEFlowList(NULL),
fCenWeightEbE(0.),
fQAZDCCuts(kFALSE),
fQAZDCCutsFlag(kTRUE),
fUseTracklets(kFALSE),
fMinMulZN(0),
fMaxDevZN(5.),
fZDCGainAlpha(0.395),
fbFlagIsPosMagField(kFALSE)
{
  // constructor
  
  // base list to hold all output objects:
  fHistList = new TList();
  fHistList->SetName("cobjQC");
  fHistList->SetOwner(kTRUE);
  
  // base list to hold all temp objects:
  fTempList = new TList();
  fTempList->SetName("temp");
  fTempList->SetOwner(kTRUE);
  
  // 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();
  
  fRunList = TArrayI();
  fAvVtxPosX = TArrayD();
  fAvVtxPosY = TArrayD();
  fAvVtxPosZ = TArrayD();
  fnEvRbR = TArrayI();
  
  // CRC
  this->InitializeCostantsForCRC();
  this->InitializeArraysForParticleWeights();
  this->InitializeArraysForCRC();
  this->InitializeArraysForCRCVZ();
  this->InitializeArraysForCRCZDC();
  this->InitializeArraysForCRC2();
  this->InitializeArraysForQVec();
  this->InitializeArraysForCRCPt();
  this->InitializeArraysForCME();
  this->InitializeArraysForFlowEbE();
  this->InitializeArraysForFlowQC();
  this->InitializeArraysForFlowQCHighOrders();
  this->InitializeArraysForFlowGF();
  this->InitializeArraysForFlowSPZDC();
  this->InitializeArraysForFlowSPVZ();
  this->InitializeArraysForEbEFlow();
  
  // printf("Arrays initialized \n");
  
} // end of constructor

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

AliFlowAnalysisCRC::~AliFlowAnalysisCRC()
{
  // destructor
  delete fHistList;
  delete fTempList;
  if(fCRCQVecWeightsList) delete fCRCQVecWeightsList;
  if(fCRCZDCCalibList)    delete fCRCZDCCalibList;
  if(fCRCZDC2DCutList)    delete fCRCZDC2DCutList;
  if(fCRCVZEROCalibList)  delete fCRCVZEROCalibList;
  if(fCRCZDCResList)      delete fCRCZDCResList;
  if(fZDCESEList)         delete fZDCESEList;
  delete[] fCRCPtvarPtBins;
  delete[] fCRCPtBins;
  delete[] fZDCEPweightEbE;
  delete[] fCorrMap;
  delete[] fchisqVA;
  delete[] fchisqVC;
  if(fPhiExclZoneHist) delete fPhiExclZoneHist;
} // end of AliFlowAnalysisCRC::~AliFlowAnalysisCRC()

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

void AliFlowAnalysisCRC::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);
  fRandom = new TRandom3(0); // if uiSeed is 0, the seed is determined uniquely in space and time via TUUID
  
  // 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->BookEverythingForNestedLoops();
  this->BookEverythingForMixedHarmonics();
  this->BookEverythingForControlHistograms();
  this->BookEverythingForBootstrap();
  this->SetRunList();
  if(fCalculateCRC) {
    this->BookEverythingForCRC();
    this->BookEverythingForCRCVZ();
    this->BookEverythingForCRCZDC();
    this->BookEverythingForCRCPt();
    this->BookEverythingForCRC2();
    this->BookEverythingForQVec();
    this->BookEverythingForCME();
    this->BookEverythingForFlowEbE();
    this->BookEverythingForFlowQC();
    this->BookEverythingForFlowQCHighOrders();
    this->BookEverythingForFlowGF();
    this->BookEverythingForFlowSPZDC();
    this->BookEverythingForFlowSPVZ();
    this->BookEverythingForEbEFlow();
  }
  this->BookEverythingForVarious();
  
  this->SetCentralityWeights();
  
  // 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();
  if(fCalculateCRC) {
    // j) Store CRC flags:
    this->StoreCRCFlags();
  }
  
  TH1::AddDirectory(oldHistAddStatus);
  
  // printf("Stuff booked \n");
  
} // end of void AliFlowAnalysisCRC::Init()

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

void AliFlowAnalysisCRC::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};
  // d.1) initialize particle weights
  // 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 various;
  // 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
  Double_t wPhiEta = 1.;
  Double_t wt = 1.;
  Double_t wProbCut = 1.;
  Int_t nCounterNoRPs = 0; // needed only for shuffling
  fNumberOfRPsEBE = anEvent->GetNumberOfRPs(); // number of RPs (i.e. number of reference particles)
  fNumberOfPOIsEBE = anEvent->GetNumberOfPOIs(); // number of POIs (i.e. number of particles of interest)
  fReferenceMultiplicityEBE = anEvent->GetReferenceMultiplicity(); // reference multiplicity for current event
  fCentralityEBE = anEvent->GetCentrality(); // centrality percentile for current event
  fCentralityCL1EBE = anEvent->GetCentralityCL1(); // centrality percentile for current event (alternative estimation)
  fCentralityTRKEBE = anEvent->GetCentralityTRK(); // centrality percentile for current event (alternative estimation)
  fNITSCL1EBE = anEvent->GetNITSCL1();
  
//  printf("begin AliFlowAnalysisCRC::Make \n");
  
  if(fExactNoRPs > 0 && fNumberOfRPsEBE<fExactNoRPs){return;}
  if(!fCentralityEBE){return;}
  if(!fNumberOfRPsEBE || !fNumberOfPOIsEBE){return;}
  fhCenvsMul[0]->Fill(fCentralityEBE,fReferenceMultiplicityEBE);
  if((fDataSet==k2015 || fDataSet==k2015v6 || fDataSet==k2015pidfix) && !MultCut2015o()){return;}
  
  if(!fRefMultRbRPro) {
    fReferenceMultiplicityRecEBE = fReferenceMultiplicityEBE-fMultCutAv->GetBinContent(fMultCutAv->FindBin(fCentralityEBE));
  } else {
    Int_t runbin = fRefMultRbRPro->GetXaxis()->FindBin(Form("%d",fRunNum));
    Int_t cenbin = fRefMultRbRPro->GetYaxis()->FindBin(fCentralityEBE);
    fReferenceMultiplicityRecEBE = fReferenceMultiplicityEBE-fRefMultRbRPro->GetBinContent(runbin,cenbin);
  }
  
  // centrality flattening with weights
  fCenWeightEbE = 1.;
  if(fCenWeigCalHist) fCenWeightEbE = fCenWeigCalHist->GetBinContent(fCenWeigCalHist->FindBin(fCentralityEBE));
  
  // primary vertex position (x,y,z)
  anEvent->GetVertexPosition(fVtxPos);
  if(fAvVtxPosX[fRunBin]) fVtxPosCor[0] = fVtxPos[0]-fAvVtxPosX[fRunBin];
  if(fAvVtxPosY[fRunBin]) fVtxPosCor[1] = fVtxPos[1]-fAvVtxPosY[fRunBin];
  if(fAvVtxPosZ[fRunBin]) fVtxPosCor[2] = fVtxPos[2]-fAvVtxPosZ[fRunBin];
  
  Double_t ptEta[2] = {0.,0.}; // 0 = dPt, 1 = dEta
  Int_t dCharge = 0; // charge
  
  // 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
  
  // d.1) Initialize particle weights
  Int_t cw = 0;
  
  if (fDataSet==kAny) {fRunBin = 0;}
  else {fRunBin = GetCRCRunBin(fRunNum);}
  if(fRunBin<0 || fRunBin>=fCRCnRun) {return;}
  fCenBin = GetCRCCenBin(fCentralityEBE);
  if(fCenBin<0 || fCenBin>=fCRCnCen) {return;}
  
  fhCenvsMul[2]->Fill(fNumberOfRPsEBE,fNumberOfPOIsEBE);
  fhCenvsMul[3]->Fill(fReferenceMultiplicityEBE,fNumberOfRPsEBE);
  fhCenvsMul[4]->Fill(fReferenceMultiplicityEBE,fNumberOfPOIsEBE);
  if(fhAvAbsOrbit) {
    UInt_t TimeStamp = (UInt_t)anEvent->GetAbsOrbit();
    fhAvAbsOrbit->Fill(fRunBin+0.5,(Double_t)TimeStamp);
  }
  
  if(fRunNum!=fCachedRunNum) {
    fbFlagIsPosMagField = kFALSE;
    Int_t dRun15hPos[] = {246390, 246391, 246392, 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246495, 246493, 246488, 246487, 246434, 246431, 246428, 246424};
    for (Int_t i=0; i<40; i++) {
      if(fRunNum==dRun15hPos[i]) fbFlagIsPosMagField = kTRUE;
    }
  }
  
  // VZERO *********************************************************************************************************
  
  if(fUseVZERO) {
    
    for(Int_t h=0; h<fCRCnHar; h++) {
      // Get Q vectors for the subevents
      AliFlowVector vQarray[2];
      anEvent->GetV02Qsub(vQarray,h+1);
      fVZFlowVect[0][h] = vQarray[0];
      fVZFlowVect[1][h] = vQarray[1];
      // re-center VZERO Q-vectors
      if(fCRCVZEROCalibList) this->RecenterCRCQVecVZERO();
      
      // fill Q-vector RbR
//      if(fCRCVZEROQVec[fRunBin][h]) {
//        fCRCVZEROQVec[fRunBin][h]->Fill(0.5,fCentralityEBE,fVZFlowVect[0][h].X());
//        fCRCVZEROQVec[fRunBin][h]->Fill(1.5,fCentralityEBE,fVZFlowVect[0][h].Y());
//        fCRCVZEROQVec[fRunBin][h]->Fill(2.5,fCentralityEBE,fVZFlowVect[1][h].X());
//        fCRCVZEROQVec[fRunBin][h]->Fill(3.5,fCentralityEBE,fVZFlowVect[1][h].Y());
//        fCRCVZEROQVec[fRunBin][h]->Fill(4.5,fCentralityEBE,fVZFlowVect[0][h].X()*fVZFlowVect[1][h].X());
//        fCRCVZEROQVec[fRunBin][h]->Fill(5.5,fCentralityEBE,fVZFlowVect[0][h].Y()*fVZFlowVect[1][h].Y());
//        fCRCVZEROQVec[fRunBin][h]->Fill(6.5,fCentralityEBE,fVZFlowVect[0][h].X()*fVZFlowVect[1][h].Y());
//        fCRCVZEROQVec[fRunBin][h]->Fill(7.5,fCentralityEBE,fVZFlowVect[0][h].Y()*fVZFlowVect[1][h].X());
//      }
    } // end of for(Int_t h=0; h<fCRCnHar; h++)
  } // end of if(fUseVZERO)
  
  // ZDC *********************************************************************************************************
  
  if(fUseZDC) {
    // Get Q vectors for the subevents
    AliFlowVector vQarray[2];
    anEvent->GetZDC2Qsub(vQarray);
    fZDCFlowVect[0] = vQarray[0];
    fZDCFlowVect[1] = vQarray[1];
    fZNCQ0 = anEvent->GetZNCQ0()/fEnNucl;
    fZNAQ0 = anEvent->GetZNAQ0()/fEnNucl;
    fZNCen = anEvent->GetZNCEnergy()/fEnNucl;
    fZNAen = anEvent->GetZNAEnergy()/fEnNucl;
  } // end of if(fUseZDC)
  
  this->CalculateCRCQVec();
  
  if(fUseZDC) {
    this->PassQAZDCCuts();
    if(fRecenterZDC) {
      this->RecenterCRCQVecZDC();
    }
  }
  // ZDC-C (eta < -8.8)
  Double_t ZCRe = fZDCFlowVect[0].X();
  Double_t ZCIm = fZDCFlowVect[0].Y();
  Double_t ZCM  = fZDCFlowVect[0].GetMult();
  // ZDC-A (eta > 8.8)
  Double_t ZARe = fZDCFlowVect[1].X();
  Double_t ZAIm = fZDCFlowVect[1].Y();
  Double_t ZAM  = fZDCFlowVect[1].GetMult();
  if( fInvertZDC ) ZARe = -ZARe;
  
  // ZDC QA cuts
  Bool_t bPassZDCcuts = kTRUE;
  if( ZCM<=0. || ZAM<=0. || sqrt(ZCRe*ZCRe+ZCIm*ZCIm)<1.E-6 || sqrt(ZARe*ZARe+ZAIm*ZAIm)<1.E-6 ) bPassZDCcuts=kFALSE;
  if( !std::isfinite(fZDCFlowVect[0].Mod()) || !std::isfinite(fZDCFlowVect[1].Mod())) bPassZDCcuts=kFALSE;
  if(fQAZDCCuts && !fQAZDCCutsFlag) bPassZDCcuts=kFALSE;
  if(bPassZDCcuts) fEventCounter->Fill(1.5);
  
  // EbE flow *********************************************************************************************************
  
  if(fCalculateEbEFlow) {
    if(fNumberOfPOIsEBE>0 && fNumberOfPOIsEBE<=2500) {
      fEBEFlowMulBin = (Int_t)(fNumberOfPOIsEBE/100);
    } else {
      fEBEFlowMulBin = -1;
    }
  }
  
  // run-by-run corrections ********************************************************************************************
  
  if(fRunNum!=fCachedRunNum) {
    if(fPOIExtraWeights==AliFlowAnalysisCRC::kEtaPhiRbR) {
      if(fWeightsList->FindObject(Form("fCRCQVecPhiHistRbR[%d]",fRunNum))) {
        fPhiEtaRbRWeights = (TH3D*)(fWeightsList->FindObject(Form("fCRCQVecPhiRbRHist[%d]",fRunNum)));
      } else {
        AliWarning("WARNING: POIExtraWeights (kEtaPhiRbR) not found ! \n");
      }
    }
    if(fPOIExtraWeights==AliFlowAnalysisCRC::kEtaPhiChRbR) {
      for (Int_t i=0; i<2; i++) {
        if(fWeightsList->FindObject(Form("fCRCQVecPhiRbRHistCh[%d][%d]",fRunNum,i))) {
          fPhiEtaRbRWeightsCh[i] = (TH3D*)(fWeightsList->FindObject(Form("fCRCQVecPhiRbRHistCh[%d][%d]",fRunNum,i)));
        } else {
          AliWarning("WARNING: POIExtraWeights (kEtaPhiChRbR) not found ! \n");
        }
      }
    }
    if(fPOIExtraWeights==AliFlowAnalysisCRC::kEtaPhiVtxRbR) {
      for (Int_t cb=0; cb<fCRCnCen; cb++) {
        if(fWeightsList->FindObject(Form("CRCQVecPhiHistVtx[%d][%d]",cb,fRunNum))) {
          fPhiEtaWeightsVtx[cb] = (TH3D*)(fWeightsList->FindObject(Form("CRCQVecPhiHistVtx[%d][%d]",cb,fRunNum)));
        } else {
          AliWarning("WARNING: POIExtraWeights (kEtaPhiVtxRbR) not found ! \n");
        }
      }
    }
  }
  
  // loop over particles **********************************************************************************************
  
  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() || aftsTrack->InPOISelection(2))){continue;} // safety measure: consider only tracks which are RPs or POIs
      
      Bool_t IsSplitMergedTracks = kFALSE;
      if(fRemoveSplitMergedTracks) {
        IsSplitMergedTracks = EvaulateIfSplitMergedTracks(anEvent,aftsTrack,i);
      }
      if(IsSplitMergedTracks) continue;
      
      // RPs *********************************************************************************************************
      
      if(aftsTrack->InRPSelection()) {
        nCounterNoRPs++;
        dPhi = aftsTrack->Phi();
        dPt  = aftsTrack->Pt();
        dEta = aftsTrack->Eta();
        dCharge = aftsTrack->Charge();
        
        if(fSelectCharge==kPosCh && dCharge<0.) continue;
        if(fSelectCharge==kNegCh && dCharge>0.) continue;
        
        cw = (dCharge > 0. ? 0 : 1);
        wPhi = 1.;
        wPt  = 1.;
        wEta = 1.;
        wTrack = 1.;
        wPhiEta = 1.;
        wProbCut = 0.;
        // pT weights
        if(fUsePtWeights && fPtWeightsHist[fCenBin]) {
          if(dPt>fPtWeightsHist[fCenBin]->GetXaxis()->GetXmin() && dPt<fPtWeightsHist[fCenBin]->GetXaxis()->GetXmax()) wt = fPtWeightsHist[fCenBin]->Interpolate(dPt);
          else if(dPt<fPtWeightsHist[fCenBin]->GetXaxis()->GetXmin())  wt = fPtWeightsHist[fCenBin]->Interpolate(fPtWeightsHist[fCenBin]->GetXaxis()->GetXmin());
          else if(dPt>fPtWeightsHist[fCenBin]->GetXaxis()->GetXmax())  wt = fPtWeightsHist[fCenBin]->Interpolate(fPtWeightsHist[fCenBin]->GetXaxis()->GetXmax());
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        
        // extra weights: eta, phi, ch, vtx
        if(fPOIExtraWeights==kEtaPhi && fPhiEtaWeights) // determine phieta weight for POI:
        {
          wt = fPhiEtaWeights->GetBinContent(fPhiEtaWeights->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        if(fPOIExtraWeights==kEtaPhiCh && fPhiEtaWeightsCh[cw]) // determine phieta weight for POI, ch dep:
        {
          wt = fPhiEtaWeightsCh[cw]->GetBinContent(fPhiEtaWeightsCh[cw]->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        if((fPOIExtraWeights==kEtaPhiVtx || fPOIExtraWeights==kEtaPhiVtxRbR) && fPhiEtaWeightsVtx[fCenBin]) // determine phieta weight for POI:
        {
          wt = fPhiEtaWeightsVtx[fCenBin]->GetBinContent(fPhiEtaWeightsVtx[fCenBin]->FindBin(fVtxPosCor[2],dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        Int_t ptbebe = (dPt>1.? 2 : (dPt>0.5 ? 1 : 0)); // hardcoded
        if(fPOIExtraWeights==kEtaPhiChPt && fPhiEtaWeightsChPt[cw][ptbebe]) // determine phieta weight for POI, ch dep:
        {
          wt = fPhiEtaWeightsChPt[cw][ptbebe]->GetBinContent(fPhiEtaWeightsChPt[cw][ptbebe]->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        // run-by-run
        if(fPOIExtraWeights==kEtaPhiRbR && fPhiEtaRbRWeights) // determine phieta weight for POI:
        {
          wt = fPhiEtaRbRWeights->GetBinContent(fPhiEtaRbRWeights->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        if(fPOIExtraWeights==kEtaPhiChRbR && fPhiEtaRbRWeightsCh[cw]) // determine phieta weight for POI, ch dep:
        {
          wt = fPhiEtaRbRWeightsCh[cw]->GetBinContent(fPhiEtaRbRWeightsCh[cw]->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        
        if(fPhiExclZoneHist) {
          if(fPhiExclZoneHist->GetBinContent(fPhiExclZoneHist->FindBin(dEta,dPhi))<0.5) continue;
        }
        
        // 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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1)*n*dPhi);
            (*fImQ)(m,k)+=pow(wPhiEta*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(wPhiEta*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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
                  fImRPQ1dEBE[0][pe][m][k]->Fill(ptEta[pe],pow(wPhiEta*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(wPhiEta*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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
                fImRPQ2dEBE[0][m][k]->Fill(dPt,dEta,pow(wPhiEta*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(wPhiEta*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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
                    fImRPQ1dEBE[2][pe][m][k]->Fill(ptEta[pe],pow(wPhiEta*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(wPhiEta*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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
                  fImRPQ2dEBE[2][m][k]->Fill(dPt,dEta,pow(wPhiEta*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(wPhiEta*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())
      
      // POIs ********************************************************************************************************
      
      if(aftsTrack->InPOISelection() || aftsTrack->InPOISelection(2)) {
        
        if(!fUseTracklets && aftsTrack->InPOISelection(2)) continue;
        if(fUseTracklets && !aftsTrack->InPOISelection(2)) continue;
        
        dPhi = aftsTrack->Phi();
        dPt  = aftsTrack->Pt();
        dEta = aftsTrack->Eta();
        dCharge = aftsTrack->Charge();
        Int_t ITStype = aftsTrack->ITStype();
        
        if(fSelectCharge==kPosCh && dCharge<0.) continue;
        if(fSelectCharge==kNegCh && dCharge>0.) continue;
        
        if(fCRCTestSin) {
          if(dCharge > 0.) dPt += 1.E-2;
          else dPt -= 1.E-2;
        }
        
        cw = (dCharge > 0. ? 0 : 1);
        wPhi = 1.;
        wPt  = 1.;
        wEta = 1.;
        wTrack = 1.;
        wPhiEta = 1.;
        wProbCut = 0.;
        
        // pT weights
        if(fUsePtWeights && fPtWeightsHist[fCenBin]) {
          if(dPt>fPtWeightsHist[fCenBin]->GetXaxis()->GetXmin() && dPt<fPtWeightsHist[fCenBin]->GetXaxis()->GetXmax()) wt = fPtWeightsHist[fCenBin]->Interpolate(dPt);
          else if(dPt<fPtWeightsHist[fCenBin]->GetXaxis()->GetXmin())  wt = fPtWeightsHist[fCenBin]->Interpolate(fPtWeightsHist[fCenBin]->GetXaxis()->GetXmin());
          else if(dPt>fPtWeightsHist[fCenBin]->GetXaxis()->GetXmax())  wt = fPtWeightsHist[fCenBin]->Interpolate(fPtWeightsHist[fCenBin]->GetXaxis()->GetXmax());
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        
        // extra weights: eta, phi, ch, vtx
        if(fPOIExtraWeights==kEtaPhi && fPhiEtaWeights) // determine phieta weight for POI:
        {
          wt = fPhiEtaWeights->GetBinContent(fPhiEtaWeights->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        if(fPOIExtraWeights==kEtaPhiCh && fPhiEtaWeightsCh[cw]) // determine phieta weight for POI, ch dep:
        {
          wt = fPhiEtaWeightsCh[cw]->GetBinContent(fPhiEtaWeightsCh[cw]->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        if((fPOIExtraWeights==kEtaPhiVtx || fPOIExtraWeights==kEtaPhiVtxRbR) && fPhiEtaWeightsVtx[fCenBin]) // determine phieta weight for POI:
        {
          wt = fPhiEtaWeightsVtx[fCenBin]->GetBinContent(fPhiEtaWeightsVtx[fCenBin]->FindBin(fVtxPosCor[2],dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        Int_t ptbebe = (dPt>1.? 2 : (dPt>0.5 ? 1 : 0)); // hardcoded
        if(fPOIExtraWeights==kEtaPhiChPt && fPhiEtaWeightsChPt[cw][ptbebe]) // determine phieta weight for POI, ch dep:
        {
          wt = fPhiEtaWeightsChPt[cw][ptbebe]->GetBinContent(fPhiEtaWeightsChPt[cw][ptbebe]->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        // run-by-run
        if(fPOIExtraWeights==kEtaPhiRbR && fPhiEtaRbRWeights) // determine phieta weight for POI:
        {
          wt = fPhiEtaRbRWeights->GetBinContent(fPhiEtaRbRWeights->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        if(fPOIExtraWeights==kEtaPhiChRbR && fPhiEtaRbRWeightsCh[cw]) // determine phieta weight for POI, ch dep:
        {
          wt = fPhiEtaRbRWeightsCh[cw]->GetBinContent(fPhiEtaRbRWeightsCh[cw]->FindBin(fCentralityEBE,dPhi,dEta));
          if(std::isfinite(1./wt)) wPhiEta *= 1./wt;
        }
        
        if(fUsePhiEtaCuts)
        {
          // test: remove region with low SPD efficiency
          if(dPhi>2.136283 && dPhi<2.324779) continue;
        }
        
        // Generic Framework: 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<21;m++) // to be improved - hardwired 6
        {
          for(Int_t k=0;k<9;k++) // to be improved - hardwired 9
          {
            (*fReQGF)(m,k) += pow(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos(m*dPhi);
            (*fImQGF)(m,k) += pow(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Sin(m*dPhi);
          }
        }
        
        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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
                fImRPQ1dEBE[1][pe][m][k]->Fill(ptEta[pe],pow(wPhiEta*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(wPhiEta*wPhi*wPt*wEta*wTrack,k)*TMath::Cos((m+1.)*n*dPhi),1.);
              fImRPQ2dEBE[1][m][k]->Fill(dPt,dEta,pow(wPhiEta*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
        
        // Charge-Rapidity Correlations
        for (Int_t h=0;h<fCRCnHar;h++) {
          
          fCRCQRe[cw][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
          fCRCQIm[cw][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
          fCRCMult[cw][h]->Fill(dEta,wPhiEta);
          
          fCRC2QRe[cw][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
          fCRC2QIm[cw][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
          fCRC2Mul[cw][h]->Fill(dEta,pow(wPhiEta,h));
          
          fCRCZDCQRe[cw][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
          fCRCZDCQIm[cw][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
          fCRCZDCMult[cw][h]->Fill(dEta,wPhiEta);
          
          if(fRandom->Integer(2)>0.5) {
            fCRC2QRe[2][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
            fCRC2QIm[2][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
            fCRC2Mul[2][h]->Fill(dEta,pow(wPhiEta,h));
          }
          
          if(fRandom->Integer(2)>0.5) {
            fCRCZDCQRe[2][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
            fCRCZDCQIm[2][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
            fCRCZDCMult[2][h]->Fill(dEta,wPhiEta);
          } else {
            fCRCZDCQRe[3][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
            fCRCZDCQIm[3][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
            fCRCZDCMult[3][h]->Fill(dEta,wPhiEta);
          }
          
          if(fCalculateCME) {
            Double_t SpecWeig = 1.;
            if(fUseZDCESESpecWeights && fZDCESESpecWeightsHist[fZDCESEclEbE] && dPt>0.2 && dPt<20.2) {
              Double_t weraw = fZDCESESpecWeightsHist[fZDCESEclEbE]->GetBinContent(fZDCESESpecWeightsHist[fZDCESEclEbE]->FindBin(fCentralityEBE,dPt));
              if(weraw > 0.) SpecWeig = 1./weraw;
            }
            fCMEQRe[cw][h]->Fill(dEta,SpecWeig*wPhiEta*TMath::Cos((h+1.)*dPhi));
            fCMEQIm[cw][h]->Fill(dEta,SpecWeig*wPhiEta*TMath::Sin((h+1.)*dPhi));
            fCMEMult[cw][h]->Fill(dEta,SpecWeig*wPhiEta);
            fCMEQRe[2+cw][h]->Fill(dEta,pow(SpecWeig*wPhiEta,2.)*TMath::Cos((h+1.)*dPhi));
            fCMEQIm[2+cw][h]->Fill(dEta,pow(SpecWeig*wPhiEta,2.)*TMath::Sin((h+1.)*dPhi));
            fCMEMult[2+cw][h]->Fill(dEta,pow(SpecWeig*wPhiEta,2.));
            
            // spectra
            fhCenvsSpec[fZDCESEclEbE]->Fill(fCentralityEBE,dPt,SpecWeig*wPhiEta);
            fhCenvsSpec[fZDCESEnCl]->Fill(fCentralityEBE,dPt,SpecWeig*wPhiEta);
          }
          
        } // end of for (Int_t h=0;h<fCRCnHar;h++)
        
        // Flow SP ZDC
        Bool_t bFillDis=kTRUE;
        for (Int_t k=0; k<fQVecPower; k++) {
          for (Int_t h=0;h<fFlowNHarmMax;h++) {
            
            if(fFlowQCDeltaEta>0.) {
              
              fPOIPtDiffQRe[k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
              fPOIPtDiffQIm[k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
              fPOIPtDiffMul[k][h]->Fill(dPt,pow(wPhiEta,k));
              
              fPOIPtDiffQReCh[cw][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
              fPOIPtDiffQImCh[cw][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
              fPOIPtDiffMulCh[cw][k][h]->Fill(dPt,pow(wPhiEta,k));
              
              fPOIPhiDiffQRe[k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
              fPOIPhiDiffQIm[k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
              fPOIPhiDiffMul[k][h]->Fill(dPhi,pow(wPhiEta,k));
              
              fPOIPhiEtaDiffQRe[ITStype][k][h]->Fill(dPhi,dEta,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
              fPOIPhiEtaDiffQIm[ITStype][k][h]->Fill(dPhi,dEta,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
              fPOIPhiEtaDiffMul[ITStype][k][h]->Fill(dPhi,dEta,pow(wPhiEta,k));
              
              if(fabs(dEta)>fFlowQCDeltaEta/2.) {
                Int_t keta = (dEta<0.?0:1);
                fPOIPtDiffQReEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                fPOIPtDiffQImEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                fPOIPtDiffMulEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k));
                fPOIPhiDiffQReEG[keta][k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                fPOIPhiDiffQImEG[keta][k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                fPOIPhiDiffMulEG[keta][k][h]->Fill(dPhi,pow(wPhiEta,k));
              }
              
            } else if(fFlowQCDeltaEta<0. && fFlowQCDeltaEta>-1.) {
              
              if(dEta>0.) {
                fPOIPtDiffQRe[k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                fPOIPtDiffQIm[k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                fPOIPtDiffMul[k][h]->Fill(dPt,pow(wPhiEta,k));
                
                fPOIPhiDiffQRe[k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                fPOIPhiDiffQIm[k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                fPOIPhiDiffMul[k][h]->Fill(dPhi,pow(wPhiEta,k));
                
                Double_t boundetagap = fabs(fFlowQCDeltaEta);
                
                if((dEta>0. && dEta<0.4-boundetagap/2.) || (dEta>0.4+boundetagap/2. && dEta<0.8)) {
                  Int_t keta;
                  if(dEta>0. && dEta<0.4-boundetagap/2.) keta = 0;
                  else keta = 1;
                  fPOIPtDiffQReEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                  fPOIPtDiffQImEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                  fPOIPtDiffMulEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k));
                }
              } else {
                bFillDis = kFALSE;
                continue;
              }
              
            } else if(fFlowQCDeltaEta<-1. && fFlowQCDeltaEta>-2.) {
              
              if(dEta<0.) {
                fPOIPtDiffQRe[k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                fPOIPtDiffQIm[k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                fPOIPtDiffMul[k][h]->Fill(dPt,pow(wPhiEta,k));
                
                fPOIPhiDiffQRe[k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                fPOIPhiDiffQIm[k][h]->Fill(dPhi,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                fPOIPhiDiffMul[k][h]->Fill(dPhi,pow(wPhiEta,k));
                
                Double_t boundetagap = fabs(fFlowQCDeltaEta)-1.;
                
                if((dEta<0. && dEta>-0.4+boundetagap/2.) || (dEta<-0.4-boundetagap/2. && dEta>-0.8)) {
                  Int_t keta;
                  if(dEta<0. && dEta>-0.4+boundetagap/2.) keta = 0;
                  else keta = 1;
                  fPOIPtDiffQReEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Cos((h+1.)*dPhi));
                  fPOIPtDiffQImEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k)*TMath::Sin((h+1.)*dPhi));
                  fPOIPtDiffMulEG[keta][k][h]->Fill(dPt,pow(wPhiEta,k));
                }
              } else {
                bFillDis = kFALSE;
                continue;
              }
              
            }
          }
        }
        
        for (Int_t h=0;h<fFlowNHarmMax;h++) {
          fEtaDiffQRe[cw][h]->Fill(dEta,wPhiEta*TMath::Cos((h+1.)*dPhi));
          fEtaDiffQIm[cw][h]->Fill(dEta,wPhiEta*TMath::Sin((h+1.)*dPhi));
          fEtaDiffMul[cw][h]->Fill(dEta,pow(wPhiEta,h+1));
          fPOIEtaPtQRe[cw][h]->Fill(dEta,dPt,wPhiEta*TMath::Cos((h+1.)*dPhi));
          fPOIEtaPtQIm[cw][h]->Fill(dEta,dPt,wPhiEta*TMath::Sin((h+1.)*dPhi));
          fPOIEtaPtMul[cw][h]->Fill(dEta,dPt,wPhiEta);
        }
        
        if(fCalculateEbEFlow) {
          if(fEBEFlowMulBin>=0) fEbEFlowAzimDis[fEBEFlowMulBin]->Fill(dPhi,wPhiEta);
        }
        
        // WARNING: do not invert order of SPZDC and QC, used in SC
        if(bPassZDCcuts && fCalculateFlowZDC && fUseZDC) this->CalculateFlowSPZDC(ZCRe,ZCIm,ZARe,ZAIm,dPhi,dEta,dPt,wPhiEta,dCharge,i);
        
        if(bFillDis && bPassZDCcuts && fCalculateFlowZDC && fUseZDC) {
          
          fFlowSPZDCv1etaPro[fCenBin][0][7]->Fill(dEta,TMath::Cos(dPhi)*ZARe+TMath::Sin(dPhi)*ZAIm,wPhiEta);
          fFlowSPZDCv1etaPro[fCenBin][0][8]->Fill(dEta,TMath::Cos(dPhi)*ZCRe+TMath::Sin(dPhi)*ZCIm,wPhiEta);
          if(cw==0) {
            fFlowSPZDCv1etaPro[fCenBin][0][9]->Fill(dEta,TMath::Cos(dPhi)*ZARe+TMath::Sin(dPhi)*ZAIm,wPhiEta);
            fFlowSPZDCv1etaPro[fCenBin][0][10]->Fill(dEta,TMath::Cos(dPhi)*ZCRe+TMath::Sin(dPhi)*ZCIm,wPhiEta);
          } else {
            fFlowSPZDCv1etaPro[fCenBin][0][11]->Fill(dEta,TMath::Cos(dPhi)*ZARe+TMath::Sin(dPhi)*ZAIm,wPhiEta);
            fFlowSPZDCv1etaPro[fCenBin][0][12]->Fill(dEta,TMath::Cos(dPhi)*ZCRe+TMath::Sin(dPhi)*ZCIm,wPhiEta);
          }
          
        }
        
        // all runs
        fCRCQVecPhiHist->Fill(fCentralityEBE,dPhi,dEta,wPhiEta);
        fCRCQVecPhiHistCh[cw]->Fill(fCentralityEBE,dPhi,dEta,wPhiEta);
        fCRCQVecPhiHistChPt[cw][ptbebe]->Fill(fCentralityEBE,dPhi,dEta,wPhiEta);
        for (Int_t h=0;h<6;h++) {
          fCRCQVecHarCosProCh[cw]->Fill(fCentralityEBE,(Double_t)h+0.5,dEta,TMath::Cos((h+1.)*dPhi),wPhiEta);
          fCRCQVecHarSinProCh[cw]->Fill(fCentralityEBE,(Double_t)h+0.5,dEta,TMath::Sin((h+1.)*dPhi),wPhiEta);
        }
        Double_t FillCw = (fbFlagIsPosMagField==kTRUE?(cw==0?0.5:1.5):(cw==0?2.5:3.5));
        if(fCentralityEBE>5. && fCentralityEBE<40.) {
          fCRCQVecPtHistMagField->Fill(dPt,FillCw,wPhiEta);
        }
        if(fVtxRbR) fCRCQVecPhiHistVtx[fCenBin][fRunBin]->Fill(fVtxPosCor[2],dPhi,dEta,wPhiEta);
        
        fFlowQCSpectraPubBin->Fill(fCentralityEBE,dPt,wPhiEta*fCenWeightEbE);
        fFlowQCSpectraCharge[cw]->Fill(fCentralityEBE,dPt,wPhiEta*fCenWeightEbE);
        
      } // 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);
    }
  }
  
  // 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)
  
  // i.2) Calculate CRC quantities:
  if(fCalculateCRC) {
    if(fCalculateCRCInt) this->CalculateCRCCorr();
    if(fCalculateCRC2) this->CalculateCRC2Cor();
    if(fCalculateCRCVZ && fUseVZERO) this->CalculateCRCVZERO();
    if(fCalculateCRCZDC && fUseZDC) this->CalculateCRCZDC();
    if(fCalculateCRCPt) this->CalculateCRCPtCorr();
    //    if(fUseVZERO && fUseZDC) this->CalculateVZvsZDC();
    if(fCalculateCME && fUseZDC) {
      this->CalculateCMETPC();
      this->CalculateCMEZDC();
    }
  }
  if(fCalculateFlowQC) {
    this->CalculateFlowSCZDC();
    this->CalculateFlowQC();
    this->CalculateFlowQCHighOrders();
    this->CalculateFlowGF();
  }
  if(fCalculateFlowVZ && fUseVZERO) this->CalculateFlowSPVZ();
  if(fCalculateEbEFlow) this->FitEbEFlow();
  
  // j) Distributions of correlations:
  if(fStoreDistributions){this->StoreDistributionsOfCorrelations();}
  
  // 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) Store multiplicity in various:
  if(fStoreVarious) this->FillVarious();
  
  // o) Reset all event-by-event quantities (very important !!!!):
  this->ResetEventByEventQuantities();
  
  // p) cache run number
  fCachedRunNum = fRunNum;
  
  fQAZDCCutsFlag = kTRUE;
  
//  printf("end AliFlowAnalysisCRC::Make \n");
  
} // end of AliFlowAnalysisCRC::Make(AliFlowEventSimple* anEvent)

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

void AliFlowAnalysisCRC::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;
  // n) Calculate cumulants for mixed harmonics;
  // o) Calculate charge-rapidity correlations;
  // p) Calculate cumulants for bootstrap;
  // q) Finalize various;
  
  // a) Check all pointers used in this method:
  this->CheckPointersUsedInFinish();
  
  // b) Access the constants:
  this->CommonConstants("Finish");
  
  if(fCommonHists && fCommonHists->GetHarmonic()) {
    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);
  fUsePhiEtaWeights = (Bool_t)fUseParticleWeights->GetBinContent(5);
  fUsePhiEtaWeightsChDep = (Bool_t)fUseParticleWeights->GetBinContent(6);
  fUsePhiEtaWeightsVtxDep = (Bool_t)fUseParticleWeights->GetBinContent(7);
  fUsePhiEtaWeightsChPtDep = (Bool_t)fUseParticleWeights->GetBinContent(8);
  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);
  fStoreVarious = (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);
  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)
  
  // n) Calculate cumulants for mixed harmonics:
  if(fCalculateMixedHarmonics){this->CalculateCumulantsMixedHarmonics();}
  
  // o) Calculate charge-rapidity correlations:
  if(fCalculateCRC) {
    if(fCalculateCRCInt) this->FinalizeCRCCorr();
    if(fCalculateCRC2) this->FinalizeCRC2Cor();
    if(fCalculateCRCVZ && fUseVZERO) this->FinalizeCRCVZERO();
    if(fCalculateCRCZDC && fUseZDC) this->FinalizeCRCZDC();
    if(fCalculateCRCPt) this->FinalizeCRCPtCorr();
    if(fCalculateCME && fUseZDC) {
      this->FinalizeCMETPC();
      this->FinalizeCMEZDC();
    }
  }
  // WARNING: do not invert order of SPZDC and QC, used in SC
  if(fCalculateFlowZDC && fUseZDC) this->FinalizeFlowSPZDC();
  if(fCalculateFlowQC) {
    this->FinalizeFlowQC();
    this->FinalizeFlowQCHighOrders();
    this->FinalizeFlowGF();
  }
  if(fCalculateFlowVZ && fUseVZERO) this->FinalizeFlowSPVZ();
  
  // p) Calculate cumulants for bootstrap:
  if(fUseBootstrap||fUseBootstrapVsM){this->CalculateCumulantsForBootstrap();}
  
  // q) Finalize various:
  if(fStoreVarious) this->FinalizeVarious();
  
} // end of AliFlowAnalysisCRC::Finish()

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

Bool_t AliFlowAnalysisCRC::EvaulateIfSplitMergedTracks(AliFlowEventSimple* anEvent, AliFlowTrackSimple* aftsTrack, Int_t it1)
{
  const Float_t kLimit1 = 0.02 * 3;
  Float_t bSign = (fbFlagIsPosMagField? -1 : 1);
  Int_t nTracks = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks
  
  Bool_t isNoSplit = kFALSE;
  
  //your cuts
  if (it1 < nTracks - 1) {
    for (Int_t itll2 = it1 + 1; itll2 < nTracks; itll2++) {
      
      AliFlowTrackSimple* aftsTrack2 = (AliFlowTrackSimple*)anEvent->GetTrack(itll2);
      if (!aftsTrack2) {
        delete aftsTrack2;
        continue;
      }
      
      if(!(aftsTrack2->InRPSelection() || aftsTrack2->InPOISelection() || aftsTrack2->InPOISelection(2))){continue;} // safety measure: consider only tracks which are RPs or POIs
      
      Double_t deta1 = aftsTrack->Eta() - aftsTrack2->Eta();
      
      if (TMath::Abs(deta1) < 0.02 * 2.5 * 3) {
        // phi in rad
        Float_t phi1rad1 = aftsTrack->Phi();
        Float_t phi2rad1 = aftsTrack2->Phi();
        
        // check first boundaries to see if is worth to loop and find the minimum
        Float_t dphistar11 = GetDPhiStar(phi1rad1, aftsTrack->Pt(), aftsTrack->Charge(), phi2rad1, aftsTrack2->Pt(), aftsTrack2->Charge(), 0.8, bSign);
        Float_t dphistar21 = GetDPhiStar(phi1rad1, aftsTrack->Pt(), aftsTrack->Charge(), phi2rad1, aftsTrack2->Pt(), aftsTrack2->Charge(), 2.5, bSign);
        Float_t dphistarminabs1 = 1e5;
        Float_t dphistarmin1 = 1e5;
        
        if (TMath::Abs(dphistar11) < kLimit1 || TMath::Abs(dphistar21) < kLimit1 || dphistar11 * dphistar21 < 0 ) {
          
          for (Double_t rad1 = 0.8; rad1 < 2.51; rad1 += 0.01) {
            
            Float_t dphistar1 = GetDPhiStar(phi1rad1, aftsTrack->Pt(), aftsTrack->Charge(), phi2rad1, aftsTrack2->Pt(), aftsTrack2->Charge(), rad1, bSign);
            Float_t dphistarabs1 = TMath::Abs(dphistar1);
            if (dphistarabs1 < dphistarminabs1) {
              dphistarmin1 = dphistar1;
              dphistarminabs1 = dphistarabs1;
            }
          }
          if (dphistarminabs1 < 0.02 && TMath::Abs(deta1) < 0.02) {
            //AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs,dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
            isNoSplit = kTRUE;
          }
        }
        
      }
      
    }
  }
  
  return isNoSplit;
}

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

void AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::EvaluateIntFlowNestedLoops(AliFlowEventSimple* anEvent)

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

void AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::EvaluateDiffFlowNestedLoops(AliFlowEventSimple* anEvent)

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

void AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::CalculateIntFlowCorrectionsForNUACosTerms()

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

void AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::CalculateIntFlowCorrectionsForNUASinTerms()

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

void AliFlowAnalysisCRC::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.
  // k) Get pointers for CRC histograms;
  
  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();
    this->GetPointersForCRC();
    this->GetPointersForCRCVZ();
    this->GetPointersForCRCZDC();
    this->GetPointersForCRCPt();
    this->GetPointersForCRC2();
    this->GetPointersForQVec();
    this->GetPointersForCME();
    this->GetPointersForFlowQC();
    this->GetPointersForFlowQCHighOrders();
    this->GetPointersForFlowGF();
    this->GetPointersForFlowSPZDC();
    this->GetPointersForFlowSPVZ();
    this->GetPointersForEbEFlow();
    this->GetPointersForVarious();
  } else
  {
    printf("\n WARNING (QC): outputListHistos is NULL in AFAWQC::GOH() !!!!\n\n");
    exit(0);
  }
  
} // end of void AliFlowAnalysisCRC::GetOutputHistograms(TList *outputListHistos)

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

TProfile* AliFlowAnalysisCRC::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* AliFlowAnalysisCRC::MakePtProjection(TProfile2D *profilePtEta)


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


TProfile* AliFlowAnalysisCRC::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* AliFlowAnalysisCRC::MakeEtaProjection(TProfile2D *profilePtEta)

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

void AliFlowAnalysisCRC::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++)
    {
      dVn[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinContent(1);
      dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinError(1);
      dVn[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinContent(1);
      dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinError(1);
      dVn[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinContent(1);
      dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinError(1);
      dVn[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinContent(1);
      dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinError(1);
    }
  } else if(type == "RP")
  {
    dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinContent(1);
    dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinError(1);
    dVn[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinContent(1);
    dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinError(1);
    dVn[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinContent(1);
    dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinError(1);
    dVn[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinContent(1);
    dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinError(1);
  } else if(type == "POI")
  {
    dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinContent(1);
    dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinError(1);
    dVn[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinContent(1);
    dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinError(1);
    dVn[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinContent(1);
    dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinError(1);
    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 AliFlowAnalysisCRC::PrintFinalResultsForIntegratedFlow(TString type="RF");

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

void AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::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:
  //  Common result histograms for QC{2}:
  TString commonHistResults2ndOrderName = "AliFlowCommonHistResults2ndOrderQC";
  commonHistResults2ndOrderName += fAnalysisLabel->Data();
  fCommonHistsResults2nd = new AliFlowCommonHistResults(commonHistResults2ndOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults2nd);
  //  Common result histograms for QC{4}:
  TString commonHistResults4thOrderName = "AliFlowCommonHistResults4thOrderQC";
  commonHistResults4thOrderName += fAnalysisLabel->Data();
  fCommonHistsResults4th = new AliFlowCommonHistResults(commonHistResults4thOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults4th);
  //  Common result histograms for QC{6}:
  TString commonHistResults6thOrderName = "AliFlowCommonHistResults6thOrderQC";
  commonHistResults6thOrderName += fAnalysisLabel->Data();
  fCommonHistsResults6th = new AliFlowCommonHistResults(commonHistResults6thOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults6th);
  //  Common result histograms for QC{8}:
  TString commonHistResults8thOrderName = "AliFlowCommonHistResults8thOrderQC";
  commonHistResults8thOrderName += fAnalysisLabel->Data();
  fCommonHistsResults8th = new AliFlowCommonHistResults(commonHistResults8thOrderName.Data(),"",fHarmonic);
  fHistList->Add(fCommonHistsResults8th);
  
} // end of void AliFlowAnalysisCRC::BookCommonHistograms()

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

void AliFlowAnalysisCRC::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 ",8,0,8);
  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->GetXaxis())->SetBinLabel(5,"w_{#phi,#eta}");
  (fUseParticleWeights->GetXaxis())->SetBinLabel(6,"w_{#phi,#eta} ch_dep");
  (fUseParticleWeights->GetXaxis())->SetBinLabel(7,"w_{#phi,#eta} vtx_dep");
  (fUseParticleWeights->GetXaxis())->SetBinLabel(8,"w_{#phi,#eta} ch_pT_dep");
  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);
  fUseParticleWeights->Fill(4.5,(Int_t)fUsePhiEtaWeights);
  fUseParticleWeights->Fill(5.5,(Int_t)fUsePhiEtaWeightsChDep);
  fUseParticleWeights->Fill(6.5,(Int_t)fUsePhiEtaWeightsVtxDep);
  fUseParticleWeights->Fill(7.5,(Int_t)fUsePhiEtaWeightsChPtDep);
  fWeightsList->Add(fUseParticleWeights);
  
  // // POIs
  // for(Int_t c=0; c<2; c++)
  // {
  //  fPhiWeightsPOIs[c] = new TH1F(Form("fPhiWeightsPOIs[%d][%d]",c,h),Form("fPhiWeightsPOIs[%d][%d]",c,h),fnBinsPhi,fPhiMin,fPhiMax);
  //  fPhiWeightsPOIs[c]->Sumw2();
  //  fEtaWeightsPOIs[c] = new TH1D(Form("fEtaWeightsPOIs[%d][%d]",c,h),Form("fEtaWeightsPOIs[%d][%d]",c,h),fnBinsEta,fEtaMin,fEtaMax);
  //  fEtaWeightsPOIs[c]->Sumw2();
  //  fPhiEtaWeightsPOIs[c] = new TH2D(Form("fPhiEtaWeightsPOIs[%d][%d]",c,h),Form("fPhiEtaWeightsPOIs[%d][%d]",c,h),fnBinsPhi,fPhiMin,fPhiMax,fnBinsEta,fEtaMin,fEtaMax);
  //  fPhiEtaWeightsPOIs[c]->Sumw2();
  //
  //  if(fUsePhiWeights)
  //  {
  //   if(fWeightsList->FindObject(Form("fPhiHist[%d][%d]",c,h)))
  //   {
  //    fPhiDistrRefPOIs[c] = dynamic_cast<TH1F*>(fWeightsList->FindObject(Form("fPhiHist[%d][%d]",c,h)));
  //    if(!fPhiDistrRefPOIs[c])
  //    {
  //     printf("\n WARNING (QC): fPhiDistrRefPOIs is NULL in AFAWQC::BAFWH() !!!!\n\n");
  //     exit(0);
  //    }
  //    if(TMath::Abs(fPhiDistrRefPOIs[c]->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(\"fPhiHist\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
  //    exit(0);
  //   }
  //  } // end of if(fUsePhiWeights)
  //
  //  if(fUsePtWeights)
  //  {
  //   if(fWeightsList->FindObject(Form("fPtHist[%d][%d]",c,h)))
  //   {
  //    fPtDistrRefPOIs[c] = dynamic_cast<TH1D*>(fWeightsList->FindObject(Form("fPtHist[%d][%d]",c,h)));
  //    if(!fPtDistrRefPOIs[c])
  //    {
  //     printf("\n WARNING (QC): fPtDistrRefPOIs is NULL in AFAWQC::BAFWH() !!!!\n\n");
  //     exit(0);
  //    }
  //    if(TMath::Abs(fPtDistrRefPOIs[c]->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(\"fPtHist\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
  //    exit(0);
  //   }
  //  } // end of if(fUsePtWeights)
  //
  //  if(fUseEtaWeights)
  //  {
  //   if(fWeightsList->FindObject(Form("fEtaHist[%d][%d]",c,h)))
  //   {
  //    fEtaDistrRefPOIs[c] = dynamic_cast<TH1D*>(fWeightsList->FindObject(Form("fEtaHist[%d][%d]",c,h)));
  //    if(!fEtaDistrRefPOIs[c])
  //    {
  //     printf("\n WARNING (QC): fEtaDistrRefPOIs is NULL in AFAWQC::BAFWH() !!!!\n\n");
  //     exit(0);
  //    }
  //    if(TMath::Abs(fEtaDistrRefPOIs[c]->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: fWeightsList->FindObject(\"fEtaHist\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
  //    exit(0);
  //   }
  //  } // end of if(fUseEtaWeights)
  //
  //
  // } // end of for(Int_t c=0; c<2; c++)
  
  //**********************************************************************************************************
  
  // RPs
  
  // fPhiWeightsRPs = new TH1F("fPhiWeightsRPs","fPhiWeightsRPs",fnBinsPhi,fPhiMin,fPhiMax);
  // fEtaWeightsRPs = new TH1D("fEtaWeightsRPs","fEtaWeightsRPs",fnBinsEta,fEtaMin,fEtaMax);
  //
  // if(fUsePhiWeights)
  // {
  //  if(fWeightsList->FindObject("fPhiDistrRPs"))
  //  {
  //   fPhiDistrRefRPs = dynamic_cast<TH1F*>(fWeightsList->FindObject("fPhiDistrRPs"));
  //   if(!fPhiDistrRefRPs)
  //   {
  //    printf("\n WARNING (QC): fPhiDistrRefRPs is NULL in AFAWQC::BAFWH() !!!!\n\n");
  //    exit(0);
  //   }
  //   if(TMath::Abs(fPhiDistrRefRPs->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(\"fPhiDistrRPs\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
  //   exit(0);
  //  }
  // } // end of if(fUsePhiWeights)
  //
  // if(fUsePtWeights)
  // {
  //  if(fWeightsList->FindObject("fPtDistrRPs"))
  //  {
  //   fPtDistrRefRPs = dynamic_cast<TH1D*>(fWeightsList->FindObject("fPtDistrRPs"));
  //   if(!fPtDistrRefRPs)
  //   {
  //    printf("\n WARNING (QC): fPtDistrRefRPs is NULL in AFAWQC::BAFWH() !!!!\n\n");
  //    exit(0);
  //   }
  //   if(TMath::Abs(fPtDistrRefRPs->GetBinWidth(1)-fPtBinWidth)>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(\"fPtDistrRPs\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
  //   exit(0);
  //  }
  // } // end of if(fUsePtWeights)
  //
  // if(fUseEtaWeights)
  // {
  //  if(fWeightsList->FindObject("fEtaDistrRPs"))
  //  {
  //   fEtaDistrRefRPs = dynamic_cast<TH1D*>(fWeightsList->FindObject("fEtaDistrRPs"));
  //   if(!fEtaDistrRefRPs)
  //   {
  //    printf("\n WARNING (QC): fEtaDistrRefRPs is NULL in AFAWQC::BAFWH() !!!!\n\n");
  //    exit(0);
  //   }
  //   if(TMath::Abs(fEtaDistrRefRPs->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: fWeightsList->FindObject(\"fEtaDistrRPs\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
  //   exit(0);
  //  }
  // } // end of if(fUseEtaWeights)
  //
  
} // end of AliFlowAnalysisCRC::BookAndFillWeightsHistograms()

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

void AliFlowAnalysisCRC::SetCentralityWeights()
{
  if(!fCenWeightsHist) return;
  fCenWeigCalHist = (TH1D*)(fCenWeightsHist->Clone("fCenWeigCalHist"));
  TF1 *CenFit = new TF1("CenFit","pol0", 0., 100.);
  if(fDataSet!=AliFlowAnalysisCRC::k2011) {
    fCenWeigCalHist->Fit("CenFit","QNR","",0.,50.);
    Double_t CenAv = CenFit->GetParameter(0);
    for(Int_t b=1; b<=fCenWeigCalHist->GetNbinsX(); b++) {
      Double_t newbin = fCenWeigCalHist->GetBinContent(b);
      if(newbin) {
        fCenWeigCalHist->SetBinContent(b,CenAv/newbin);
      } else {
        fCenWeigCalHist->SetBinContent(b,1.);
      }
    }
  }
  else {
    fCenWeigCalHist->Fit("CenFit","QNR","",0.,8.);
    Double_t CenAv = CenFit->GetParameter(0);
    fCenWeigCalHist->Fit("CenFit","QNR","",12.,50.);
    Double_t SemiCenAv = CenFit->GetParameter(0);
    for(Int_t b=1; b<=fCenWeigCalHist->GetNbinsX(); b++) {
      Double_t newbin = fCenWeigCalHist->GetBinContent(b);
      if(newbin) {
        if(b<=10) fCenWeigCalHist->SetBinContent(b,CenAv/newbin);
        if(b>10 && b<=50) fCenWeigCalHist->SetBinContent(b,SemiCenAv/newbin);
        if(b>50) fCenWeigCalHist->SetBinContent(b,1.);
      } else {
        fCenWeigCalHist->SetBinContent(b,1.);
      }
    }
  }
  fCenWeigCalHist->SetName("CenWeights");
  fVariousList->Add(fCenWeigCalHist);
} // end of AliFlowAnalysisCRC::SetCentralityWeights()

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

void AliFlowAnalysisCRC::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",20,0.,20.);
  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,"fStoreVarious");
  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");
  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);
  fReQGF = new TMatrixD(21,9);
  fImQGF = new TMatrixD(21,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 AliFlowAnalysisCRC::BookEverythingForIntegratedFlow()

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

void AliFlowAnalysisCRC::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 AliFlowAnalysisCRC::BookEverythingForControlHistograms()

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

void AliFlowAnalysisCRC::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);
  
  // 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 AliFlowAnalysisCRC::BookEverythingForBootstrap()

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

void AliFlowAnalysisCRC::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);