]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx
generatio of various plaes in on the fly
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / AliFlowAnalysisWithQCumulants.cxx
index fdc046ce3d97f70d447d4e85456d633d93867abb..31e12c1dd59ac1e1abcf92b3b5f66367afec8d5e 100644 (file)
@@ -111,6 +111,7 @@ AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():
  fEtaWeights(NULL),
  // 2b.) event weights:
  fMultiplicityWeight(NULL),
+ fMultiplicityIs(AliFlowCommonConstants::kRP),
  // 3.) integrated flow:
  fIntFlowList(NULL), 
  fIntFlowProfiles(NULL),
@@ -125,16 +126,21 @@ AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():
  fPropagateErrorAlsoFromNIT(kFALSE), 
  fCalculateCumulantsVsM(kFALSE),
  fCalculateAllCorrelationsVsM(kFALSE), 
- fMultiplicityIsRefMultiplicity(kFALSE),
  fMinimumBiasReferenceFlow(kTRUE), 
  fForgetAboutCovariances(kFALSE), 
  fStorePhiDistributionForOneEvent(kFALSE),
+ fExactNoRPs(0),
+ fUse2DHistograms(kFALSE),
+ fFillProfilesVsMUsingWeights(kTRUE),
+ fUseQvectorTerms(kFALSE),
  fReQ(NULL),
  fImQ(NULL),
  fSpk(NULL),
  fIntFlowCorrelationsEBE(NULL),
  fIntFlowEventWeightsForCorrelationsEBE(NULL),
  fIntFlowCorrelationsAllEBE(NULL),
+ fNumberOfRPsEBE(0.),
+ fNumberOfPOIsEBE(0.),
  fReferenceMultiplicityEBE(0.),  
  fAvMultiplicity(NULL),
  fIntFlowCorrelationsPro(NULL),
@@ -170,6 +176,7 @@ AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():
  fDistributionsList(NULL),
  fDistributionsFlags(NULL),
  fStoreDistributions(kFALSE),
+ fnBinsForCorrelations(10000),
  // 7.) various:
  fVariousList(NULL),
  fPhiDistributionForOneEvent(NULL),
@@ -208,7 +215,14 @@ AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():
  f7pCumulants(NULL),
  f8pCumulants(NULL),
  fMixedHarmonicProductOfEventWeights(NULL),
- fMixedHarmonicProductOfCorrelations(NULL)
+ fMixedHarmonicProductOfCorrelations(NULL),
+ // 10.) Control histograms:
+ fControlHistogramsList(NULL), 
+ fControlHistogramsFlags(NULL),
+ fStoreControlHistograms(kFALSE),
+ fCorrelationNoRPsVsRefMult(NULL), 
+ fCorrelationNoPOIsVsRefMult(NULL),
+ fCorrelationNoRPsVsNoPOIs(NULL)
  {
   // constructor  
   
@@ -222,7 +236,7 @@ AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():
   
   // multiplicity weight:
   fMultiplicityWeight = new TString("combinations");
-    
+
   // analysis label;
   fAnalysisLabel = new TString();
       
@@ -233,6 +247,7 @@ AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():
   this->InitializeArraysForVarious();
   this->InitializeArraysForNestedLoops();
   this->InitializeArraysForMixedHarmonics();
+  this->InitializeArraysForControlHistograms();
   
  } // end of constructor
  
@@ -256,7 +271,8 @@ void AliFlowAnalysisWithQCumulants::Init()
  // 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.
+ // g) Store flags for mixed harmonics;
+ // h) Store flags for control histograms.
   
  //save old value and prevent histograms from being added to directory
  //to avoid name clashes in case multiple analaysis objects are used
@@ -279,6 +295,7 @@ void AliFlowAnalysisWithQCumulants::Init()
  this->BookEverythingForVarious();
  this->BookEverythingForNestedLoops();
  this->BookEverythingForMixedHarmonics();
+ this->BookEverythingForControlHistograms();
 
  // d) Store flags for integrated and differential flow:
  this->StoreIntFlowFlags();
@@ -289,7 +306,9 @@ void AliFlowAnalysisWithQCumulants::Init()
  this->StoreHarmonic();
  // g) Store flags for mixed harmonics:
  this->StoreMixedHarmonicsFlags();
+ // h) Store flags for control histograms:
+ this->StoreControlHistogramsFlags();
+
  TH1::AddDirectory(oldHistAddStatus);
 } // end of void AliFlowAnalysisWithQCumulants::Init()
 
@@ -325,28 +344,32 @@ void AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
  Double_t wPt  = 1.; // pt weight
  Double_t wEta = 1.; // eta weight
  Double_t wTrack = 1.; // track weight
- Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of RPs (i.e. number of reference particles)
+ Int_t nCounterNoRPs = 0; // needed only for shuffling
+ fNumberOfRPsEBE = anEvent->GetNumberOfRPs(); // number of RPs (i.e. number of reference particles)
+ if(fExactNoRPs > 0 && fNumberOfRPsEBE<fExactNoRPs){return;}
+ fNumberOfPOIsEBE = anEvent->GetNumberOfPOIs(); // number of POIs (i.e. number of particles of interest)
  fReferenceMultiplicityEBE = anEvent->GetReferenceMultiplicity(); // reference multiplicity for current event
  Double_t ptEta[2] = {0.,0.}; // 0 = dPt, 1 = dEta
   
  // c) Fill the common control histograms and call the method to fill fAvMultiplicity:
  this->FillCommonControlHistograms(anEvent);                                                               
- this->FillAverageMultiplicities(nRP);                                                                  
+ this->FillAverageMultiplicities((Int_t)(fNumberOfRPsEBE)); 
+ if(fStoreControlHistograms){this->FillControlHistograms(anEvent);}                                                              
                                                                                                                                                                                                                                                                                         
  // d) Loop over data and calculate e-b-e quantities Q_{n,k}, S_{p,k} and s_{p,k}:
- Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI where:
-                                           //  nRP   = # of reference particles;
-                                           //  nPOI  = # of particles of interest.
+ Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks
  AliFlowTrackSimple *aftsTrack = NULL;
  Int_t n = fHarmonic; // shortcut for the harmonic 
  for(Int_t i=0;i<nPrim;i++) 
  { 
+  if(fExactNoRPs > 0 && nCounterNoRPs>fExactNoRPs){continue;}
   aftsTrack=anEvent->GetTrack(i);
   if(aftsTrack)
   {
    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())){continue;} // safety measure: consider only tracks which are RPs or POIs
    if(aftsTrack->InRPSelection()) // RP condition:
    {    
+    nCounterNoRPs++;
     dPhi = aftsTrack->Phi();
     dPt  = aftsTrack->Pt();
     dEta = aftsTrack->Eta();
@@ -521,29 +544,29 @@ void AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
  {
   if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
   {
-   if(nRP>1){this->CalculateIntFlowCorrelations();} // without using particle weights
+   if(fNumberOfRPsEBE>1){this->CalculateIntFlowCorrelations();} // without using particle weights
   } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
     {
-     if(nRP>1){this->CalculateIntFlowCorrelationsUsingParticleWeights();} // with using particle weights   
+     if(fNumberOfRPsEBE>1){this->CalculateIntFlowCorrelationsUsingParticleWeights();} // with using particle weights   
     }        
   // Whether or not using particle weights the following is calculated in the same way:  
-  if(nRP>3){this->CalculateIntFlowProductOfCorrelations();}
-  if(nRP>1){this->CalculateIntFlowSumOfEventWeights();}
-  if(nRP>1){this->CalculateIntFlowSumOfProductOfEventWeights();}  
+  if(fNumberOfRPsEBE>3){this->CalculateIntFlowProductOfCorrelations();}
+  if(fNumberOfRPsEBE>1){this->CalculateIntFlowSumOfEventWeights();}
+  if(fNumberOfRPsEBE>1){this->CalculateIntFlowSumOfProductOfEventWeights();}  
   // Non-isotropic terms:
   if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
   {
-   if(nRP>0){this->CalculateIntFlowCorrectionsForNUASinTerms();}
-   if(nRP>0){this->CalculateIntFlowCorrectionsForNUACosTerms();}
+   if(fNumberOfRPsEBE>0){this->CalculateIntFlowCorrectionsForNUASinTerms();}
+   if(fNumberOfRPsEBE>0){this->CalculateIntFlowCorrectionsForNUACosTerms();}
   } else // to if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights))
     {
-     if(nRP>0){this->CalculateIntFlowCorrectionsForNUASinTermsUsingParticleWeights();}
-     if(nRP>0){this->CalculateIntFlowCorrectionsForNUACosTermsUsingParticleWeights();}     
+     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(nRP>0){this->CalculateIntFlowProductOfCorrectionTermsForNUA();}     
-  if(nRP>0){this->CalculateIntFlowSumOfEventWeightsNUA();}     
-  if(nRP>0){this->CalculateIntFlowSumOfProductOfEventWeightsNUA();}     
+  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)
@@ -706,7 +729,9 @@ void AliFlowAnalysisWithQCumulants::Finish()
  fStorePhiDistributionForOneEvent = (Bool_t)fIntFlowFlags->GetBinContent(13);
  fFillMultipleControlHistograms = (Bool_t)fIntFlowFlags->GetBinContent(14); 
  fCalculateAllCorrelationsVsM = (Bool_t)fIntFlowFlags->GetBinContent(15);
- fMultiplicityIsRefMultiplicity = (Bool_t)fIntFlowFlags->GetBinContent(16);
+ 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);
@@ -842,7 +867,7 @@ void AliFlowAnalysisWithQCumulants::EvaluateIntFlowNestedLoops(AliFlowEventSimpl
 {
  // Evaluate all correlators for reference flow with nested loops.
  
- Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = nRP + nPOI 
+ Int_t nPrim = anEvent->NumberOfTracks(); // number of primaries
  if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10 
  {
   // Without using particle weights:
@@ -893,7 +918,7 @@ void AliFlowAnalysisWithQCumulants::EvaluateDiffFlowNestedLoops(AliFlowEventSimp
 
  if(!fCalculateDiffFlow){return;}
 
- Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = nRP + nPOI 
+ Int_t nPrim = anEvent->NumberOfTracks(); // number of primaries
  if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10
  {
   // Without using particle weights:
@@ -986,14 +1011,17 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUACosTerms()
 
  // Multiplicity bin of an event (relevant for all histos vs M): 
  Double_t dMultiplicityBin = 0.;
- if(!fMultiplicityIsRefMultiplicity)
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
  {
-  dMultiplicityBin = dMult+0.5;
- } else
+  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): ****
  //                                  *************************************************************
@@ -1092,13 +1120,16 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUASinTerms()
         
  // Multiplicity bin of an event (relevant for all histos vs M): 
  Double_t dMultiplicityBin = 0.;
- if(!fMultiplicityIsRefMultiplicity)
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
  {
-  dMultiplicityBin = dMult+0.5;
- } else
+  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): ****
@@ -1186,8 +1217,9 @@ void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)
  // 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 nested loops' histograms;
- // h) Get pointers for mixed harmonics histograms.
+ // g) Get pointers for mixed harmonics histograms;
+ // h) Get pointers for nested loops' histograms;
+ // i) Get pointers for control histograms.
  
  if(outputListHistos)
  {     
@@ -1205,6 +1237,7 @@ void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)
   this->GetPointersForOtherDiffCorrelators();  
   this->GetPointersForMixedHarmonicsHistograms(); 
   this->GetPointersForNestedLoopsHistograms(); 
+  this->GetPointersForControlHistograms();
  } else 
    {
     printf("\n WARNING (QC): outputListHistos is NULL in AFAWQC::GOH() !!!!\n\n");
@@ -1629,7 +1662,7 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
  // 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",16,0,16);
+ fIntFlowFlags = new TProfile(intFlowFlagsName.Data(),"Flags for Integrated Flow",20,0.,20.);
  fIntFlowFlags->SetTickLength(-0.01,"Y");
  fIntFlowFlags->SetMarkerStyle(25);
  fIntFlowFlags->SetLabelSize(0.04);
@@ -1650,7 +1683,11 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
  fIntFlowFlags->GetXaxis()->SetBinLabel(13,"fStorePhiDistributionForOneEvent");
  fIntFlowFlags->GetXaxis()->SetBinLabel(14,"fFillMultipleControlHistograms");
  fIntFlowFlags->GetXaxis()->SetBinLabel(15,"Calculate all correlations vs M");
- fIntFlowFlags->GetXaxis()->SetBinLabel(16,"fMultiplicityIsRefMultiplicity");
+ 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:
@@ -1747,8 +1784,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
                                                  fnBinsMult,fMinMult,fMaxMult,"s");   
    fIntFlowCorrelationsVsMPro[ci]->Sumw2();                                                                                       
    fIntFlowCorrelationsVsMPro[ci]->GetYaxis()->SetTitle(correlationFlag[ci].Data());
-   fIntFlowCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");} 
+   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";
@@ -1758,8 +1803,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
                                                         fnBinsMult,fMinMult,fMaxMult,"s");   
    fIntFlowSquaredCorrelationsVsMPro[ci]->Sumw2();                                                                                              
    fIntFlowSquaredCorrelationsVsMPro[ci]->GetYaxis()->SetTitle(squaredCorrelationFlag[ci].Data());
-   fIntFlowSquaredCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowSquaredCorrelationsVsMPro[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");} 
+   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)
@@ -1911,8 +1964,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
    if(fIntFlowCorrelationsAllVsMPro[n])
    {
     fIntFlowCorrelationsAllVsMPro[n]->Sumw2();
-    fIntFlowCorrelationsAllVsMPro[n]->GetXaxis()->SetTitle("M");  
-    if(fMultiplicityIsRefMultiplicity){fIntFlowCorrelationsAllVsMPro[n]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}  
+    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++)
@@ -1958,8 +2019,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
    fIntFlowProductOfCorrelationsVsMPro[pi] = new TProfile(Form("%s, %s",intFlowProductOfCorrelationsVsMProName.Data(),productFlag[pi].Data()),
                                                           Form("%s versus multiplicity",productFlag[pi].Data()),
                                                           fnBinsMult,fMinMult,fMaxMult);             
-   fIntFlowProductOfCorrelationsVsMPro[pi]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowProductOfCorrelationsVsMPro[pi]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}
+   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) 
@@ -2053,8 +2122,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
                                               Form("%s vs multiplicity",correlationFlag[ci].Data()),
                                               fnBinsMult,fMinMult,fMaxMult);                                            
    fIntFlowCorrelationsVsMHist[ci]->GetYaxis()->SetTitle(correlationFlag[ci].Data());
-   fIntFlowCorrelationsVsMHist[ci]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowCorrelationsVsMHist[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}
+   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) 
@@ -2177,8 +2254,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
                                          Form("%s vs multiplicity",covarianceFlag[ci].Data()),
                                          fnBinsMult,fMinMult,fMaxMult);
    fIntFlowCovariancesVsM[ci]->GetYaxis()->SetTitle(covarianceFlag[ci].Data());
-   fIntFlowCovariancesVsM[ci]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowCovariancesVsM[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}
+   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) 
@@ -2198,8 +2283,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
                                                        Form("%s vs multiplicity",sumFlag[power][si].Data()),
                                                        fnBinsMult,fMinMult,fMaxMult);    
     fIntFlowSumOfEventWeightsVsM[si][power]->GetYaxis()->SetTitle(sumFlag[power][si].Data());  
-    fIntFlowSumOfEventWeightsVsM[si][power]->GetXaxis()->SetTitle("M"); 
-    if(fMultiplicityIsRefMultiplicity){fIntFlowSumOfEventWeightsVsM[si][power]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");} 
+    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++)   
@@ -2218,8 +2311,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
    fIntFlowSumOfProductOfEventWeightsVsM[pi] = new TH1D(Form("%s, %s",intFlowSumOfProductOfEventWeightsVsMName.Data(),sopowFlag[pi].Data()),
                                                         Form("%s versus multiplicity",sopowFlag[pi].Data()),
                                                         fnBinsMult,fMinMult,fMaxMult); 
-   fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}
+   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++) 
@@ -2370,8 +2471,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
    fIntFlowQcumulantsVsM[co] = new TH1D(Form("%s, %s",intFlowQcumulantsVsMName.Data(),cumulantFlag[co].Data()),
                                         Form("%s vs multiplicity",cumulantFlag[co].Data()),
                                         fnBinsMult,fMinMult,fMaxMult);
-   fIntFlowQcumulantsVsM[co]->GetXaxis()->SetTitle("M"); 
-   if(fMultiplicityIsRefMultiplicity){fIntFlowQcumulantsVsM[co]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}                                    
+   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
@@ -2413,8 +2522,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
    fIntFlowVsM[co] = new TH1D(Form("%s, %s",intFlowVsMName.Data(),flowFlag[co].Data()),
                               Form("%s vs multiplicity",flowFlag[co].Data()),
                               fnBinsMult,fMinMult,fMaxMult);
-   fIntFlowVsM[co]->GetXaxis()->SetTitle("M");
-   if(fMultiplicityIsRefMultiplicity){fIntFlowVsM[co]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}                                      
+   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
@@ -2440,8 +2557,16 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
    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);
-   fIntFlowDetectorBiasVsM[ci]->GetXaxis()->SetTitle("M");    
-   if(fMultiplicityIsRefMultiplicity){fIntFlowDetectorBiasVsM[ci]->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");}                                 
+   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
@@ -2451,6 +2576,120 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
 
 //=======================================================================================================================
 
+void AliFlowAnalysisWithQCumulants::BookEverythingForControlHistograms()
+{
+ // Book all objects for control histograms.
+
+ // a) Book profile to hold all flags for control histograms;
+ // b) Book all control histograms.
+
+ // a) Book profile to hold all flags for control histograms:
+ TString controlHistogramsFlagsName = "fControlHistogramsFlags";
+ controlHistogramsFlagsName += fAnalysisLabel->Data();
+ fControlHistogramsFlags = new TProfile(controlHistogramsFlagsName.Data(),"Flags for Control Histograms",2,0,2);
+ fControlHistogramsFlags->SetTickLength(-0.01,"Y");
+ fControlHistogramsFlags->SetMarkerStyle(25);
+ fControlHistogramsFlags->SetLabelSize(0.04);
+ fControlHistogramsFlags->SetLabelOffset(0.02,"Y");
+ fControlHistogramsFlags->SetStats(kFALSE);
+ fControlHistogramsFlags->GetXaxis()->SetBinLabel(1,"fStoreControlHistograms");
+ fControlHistogramsFlags->GetXaxis()->SetBinLabel(2,"fUseQvectorTerms");
+ fControlHistogramsList->Add(fControlHistogramsFlags);
+
+ if(!fStoreControlHistograms){return;}
+
+ // b) Book all control histograms:
+ //  b1) Correlation between # RPs and ref. mult. determined centrally:
+ TString sCorrelationNoRPsVsRefMultName = "fCorrelationNoRPsVsRefMult";
+ sCorrelationNoRPsVsRefMultName += fAnalysisLabel->Data();
+ fCorrelationNoRPsVsRefMult = new TH2D(sCorrelationNoRPsVsRefMultName.Data(),"# RPs vs. Reference Multiplicity",fnBinsMult,fMinMult,fMaxMult,fnBinsMult,fMinMult,fMaxMult);
+ fCorrelationNoRPsVsRefMult->SetTickLength(-0.01,"Y");
+ fCorrelationNoRPsVsRefMult->SetLabelSize(0.04);
+ fCorrelationNoRPsVsRefMult->SetLabelOffset(0.02,"Y");
+ fCorrelationNoRPsVsRefMult->SetStats(kTRUE);
+ fCorrelationNoRPsVsRefMult->GetXaxis()->SetTitle("# RPs");
+ fCorrelationNoRPsVsRefMult->GetYaxis()->SetTitle("Reference Multiplicity");
+ fControlHistogramsList->Add(fCorrelationNoRPsVsRefMult);
+ //  b2) Correlation between # POIs and ref. mult. determined centrally:
+ TString sCorrelationNoPOIsVsRefMultName = "fCorrelationNoPOIsVsRefMult";
+ sCorrelationNoPOIsVsRefMultName += fAnalysisLabel->Data();
+ fCorrelationNoPOIsVsRefMult = new TH2D(sCorrelationNoPOIsVsRefMultName.Data(),"# POIs vs. Reference Multiplicity",fnBinsMult,fMinMult,fMaxMult,fnBinsMult,fMinMult,fMaxMult);
+ fCorrelationNoPOIsVsRefMult->SetTickLength(-0.01,"Y");
+ fCorrelationNoPOIsVsRefMult->SetLabelSize(0.04);
+ fCorrelationNoPOIsVsRefMult->SetLabelOffset(0.02,"Y");
+ fCorrelationNoPOIsVsRefMult->SetStats(kTRUE);
+ fCorrelationNoPOIsVsRefMult->GetXaxis()->SetTitle("# POIs");
+ fCorrelationNoPOIsVsRefMult->GetYaxis()->SetTitle("Reference Multiplicity");
+ fControlHistogramsList->Add(fCorrelationNoPOIsVsRefMult);
+ //  b3) Correlation between # RPs and # POIs:
+ TString sCorrelationNoRPsVsNoPOIsName = "fCorrelationNoRPsVsNoPOIs";
+ sCorrelationNoRPsVsNoPOIsName += fAnalysisLabel->Data();
+ fCorrelationNoRPsVsNoPOIs = new TH2D(sCorrelationNoRPsVsNoPOIsName.Data(),"# RPs vs. # POIs",fnBinsMult,fMinMult,fMaxMult,fnBinsMult,fMinMult,fMaxMult);
+ fCorrelationNoRPsVsNoPOIs->SetTickLength(-0.01,"Y");
+ fCorrelationNoRPsVsNoPOIs->SetLabelSize(0.04);
+ fCorrelationNoRPsVsNoPOIs->SetLabelOffset(0.02,"Y");
+ fCorrelationNoRPsVsNoPOIs->SetStats(kTRUE);
+ fCorrelationNoRPsVsNoPOIs->GetXaxis()->SetTitle("# RPs");
+ fCorrelationNoRPsVsNoPOIs->GetYaxis()->SetTitle("# POIs");
+ fControlHistogramsList->Add(fCorrelationNoRPsVsNoPOIs);
+ // b4) <2>, <4>, <6> and <8> vs multiplicity (#RPs, #POIs or external):
+ TString sCorrelation[4] = {"#LT2#GT","#LT4#GT","#LT6#GT","#LT8#GT"};
+ TString sMultiplicity = "";
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
+ {
+  sMultiplicity = "# RPs"; 
+ } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
+   {
+    sMultiplicity = "Reference multiplicity (from ESD)";
+   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
+     {
+      sMultiplicity = "# POIs";
+     }
+ for(Int_t ci=0;ci<4;ci++)
+ {
+  fCorrelation2468VsMult[ci] = new TH2D(Form("%s vs M",sCorrelation[ci].Data()),Form("%s vs M",sCorrelation[ci].Data()),fnBinsMult,fMinMult,fMaxMult,fnBinsForCorrelations,fMinValueOfCorrelation[ci],fMaxValueOfCorrelation[ci]); 
+  fCorrelation2468VsMult[ci]->SetTickLength(-0.01,"Y");
+  fCorrelation2468VsMult[ci]->SetLabelSize(0.04);
+  fCorrelation2468VsMult[ci]->SetLabelOffset(0.02,"Y");
+  fCorrelation2468VsMult[ci]->SetStats(kTRUE);
+  fCorrelation2468VsMult[ci]->GetXaxis()->SetTitle(sMultiplicity.Data());
+  fCorrelation2468VsMult[ci]->GetYaxis()->SetTitle(sCorrelation[ci].Data());
+  fControlHistogramsList->Add(fCorrelation2468VsMult[ci]); 
+ } // end of for(Int_t ci=0;ci<4;ci++)
+ // b5) <2><4>, <2><6>, <2><8>, <4><6> etc. vs multiplicity (#RPs, #POIs or external):
+ TString sCorrelationProduct[1] = {"#LT2#GT#LT4#GT"}; // TBI: add the other ones when needed first time
+ for(Int_t cpi=0;cpi<1;cpi++) // TBI: hardwired 1
+ {
+  fCorrelationProduct2468VsMult[cpi] = new TH2D(Form("%s vs M",sCorrelationProduct[cpi].Data()),Form("%s vs M",sCorrelationProduct[cpi].Data()),fnBinsMult,fMinMult,fMaxMult,fnBinsForCorrelations,fMinValueOfCorrelationProduct[cpi],fMaxValueOfCorrelationProduct[cpi]);
+  fCorrelationProduct2468VsMult[cpi]->SetTickLength(-0.01,"Y");
+  fCorrelationProduct2468VsMult[cpi]->SetLabelSize(0.04);
+  fCorrelationProduct2468VsMult[cpi]->SetLabelOffset(0.02,"Y");
+  fCorrelationProduct2468VsMult[cpi]->SetStats(kTRUE);
+  fCorrelationProduct2468VsMult[cpi]->GetXaxis()->SetTitle(sMultiplicity.Data());
+  fCorrelationProduct2468VsMult[cpi]->GetYaxis()->SetTitle(sCorrelationProduct[cpi].Data());
+  fControlHistogramsList->Add(fCorrelationProduct2468VsMult[cpi]); 
+ } // end of for(Int_t cpi=0;cpi<4;cpi++)
+ // b6) |Qn|^2/M, |Q2n|^2/M, |Qn|^4/(M(2M-1)), Re[Q2nQn^*Qn^*]/M, ... vs multiplicity (#RPs, #POIs or external)
+ if(fUseQvectorTerms)
+ {
+  TString sQvectorTerms[4] = {"#frac{|Q_{n}|^{2}}{M}","#frac{|Q_{2n}|^{2}}{M}","#frac{|Q_{n}|^{4}}{M(2M-1)}","#frac{Re[Q_{2n}Q_{n}^{*}Q_{n}^{*}]}{M^{3/2}}"}; // TBI: add the other ones when needed first time
+  for(Int_t qvti=0;qvti<4;qvti++) // TBI: hardwired 4
+  {
+   fQvectorTermsVsMult[qvti] = new TH2D(Form("%s vs M",sQvectorTerms[qvti].Data()),Form("%s vs M",sQvectorTerms[qvti].Data()),fnBinsMult,fMinMult,fMaxMult,fnBinsForCorrelations,-100.,100.); // TBI hardwired -100 and 100
+   fQvectorTermsVsMult[qvti]->SetTickLength(-0.01,"Y");
+   fQvectorTermsVsMult[qvti]->SetLabelSize(0.04);
+   fQvectorTermsVsMult[qvti]->SetLabelOffset(0.02,"Y");
+   fQvectorTermsVsMult[qvti]->SetStats(kTRUE);
+   fQvectorTermsVsMult[qvti]->GetXaxis()->SetTitle(sMultiplicity.Data());
+   fQvectorTermsVsMult[qvti]->GetYaxis()->SetTitle(sQvectorTerms[qvti].Data());
+   fControlHistogramsList->Add(fQvectorTermsVsMult[qvti]); 
+  } // end of for(Int_t qvti=0;qvti<4;qvti++)
+ } // end of if(fUseQvectorTerms)
+
+} // end of void AliFlowAnalysisWithQCumulants::BookEverythingForControlHistograms()
+
+//=======================================================================================================================
+
 void AliFlowAnalysisWithQCumulants::BookEverythingForMixedHarmonics()
 {
  // Book all objects for mixed harmonics.
@@ -3162,6 +3401,27 @@ void AliFlowAnalysisWithQCumulants::InitializeArraysForMixedHarmonics()
 
 //=======================================================================================================================
 
+void AliFlowAnalysisWithQCumulants::InitializeArraysForControlHistograms()
+{
+ // Initialize arrays of all objects relevant for control histograms.
+
+ for(Int_t ci=0;ci<4;ci++) // correlation index 
+ {
+  fCorrelation2468VsMult[ci] = NULL;    
+ }
+ for(Int_t cpi=0;cpi<1;cpi++) // correlation product index TBI: hardwired 1
+ {
+  fCorrelationProduct2468VsMult[cpi] = NULL;    
+ }
+ for(Int_t qwti=0;qwti<4;qwti++) // q-vector terms index TBI: hardwired 4
+ {
+  fQvectorTermsVsMult[qwti] = NULL;    
+ }
+
+} // end of void AliFlowAnalysisWithQCumulants::InitializeArraysForControlHistograms()
+
+//=======================================================================================================================
+
 void AliFlowAnalysisWithQCumulants::BookEverythingForNestedLoops()
 {
  // Book all objects relevant for calculations with nested loops.
@@ -3384,14 +3644,19 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
 
  // Multiplicity bin of an event (relevant for all histos vs M): 
  Double_t dMultiplicityBin = 0.;
- if(!fMultiplicityIsRefMultiplicity)
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
  {
-  dMultiplicityBin = dMult+0.5;
- } else
+   //Printf("RP multiplicity: %lf",fNumberOfRPsEBE);
+   dMultiplicityBin = fNumberOfRPsEBE+0.5;
+ } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
    {
+     //Printf("Reference multiplicity: %lf",fReferenceMultiplicityEBE);
     dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
-   }
-  
+   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
+     {
+      dMultiplicityBin = fNumberOfPOIsEBE+0.5;
+     }
  // Real parts of expressions involving various combinations of Q-vectors which appears
  // simultaneously in several equations for multiparticle correlations bellow: 
  // Re[Q_{2n}Q_{n}^*Q_{n}^*]
@@ -3580,13 +3845,13 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowCorrelationsEBE->SetBinContent(1,two1n1n); // <2>  
   // Testing other multiplicity weights:
   Double_t mWeight2p = 0.;
-  if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+  if(fMultiplicityWeight->Contains("combinations"))
   {
    mWeight2p = dMult*(dMult-1.);
-  } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+  } else if(fMultiplicityWeight->Contains("unit"))
     {
      mWeight2p = 1.;    
-    } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+    } else if(fMultiplicityWeight->Contains("multiplicity"))
       {
        mWeight2p = dMult;           
       }          
@@ -3595,16 +3860,27 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowSquaredCorrelationsPro->Fill(0.5,two1n1n*two1n1n,mWeight2p);
   if(fCalculateCumulantsVsM)
   {
-   fIntFlowCorrelationsVsMPro[0]->Fill(dMultiplicityBin,two1n1n,mWeight2p);
-   fIntFlowSquaredCorrelationsVsMPro[0]->Fill(dMultiplicityBin,two1n1n*two1n1n,mWeight2p);
-  } 
+   if(fFillProfilesVsMUsingWeights)  
+   {
+    fIntFlowCorrelationsVsMPro[0]->Fill(dMultiplicityBin,two1n1n,mWeight2p);
+    fIntFlowSquaredCorrelationsVsMPro[0]->Fill(dMultiplicityBin,two1n1n*two1n1n,mWeight2p);
+   } else
+     {
+      fIntFlowCorrelationsVsMPro[0]->Fill(dMultiplicityBin,two1n1n);
+      fIntFlowSquaredCorrelationsVsMPro[0]->Fill(dMultiplicityBin,two1n1n*two1n1n);
+     }
+  } // end of if(fCalculateCumulantsVsM)
   if(fCalculateAllCorrelationsVsM)
   {
    fIntFlowCorrelationsAllVsMPro[0]->Fill(dMultiplicityBin,two1n1n,mWeight2p);
    fIntFlowCorrelationsAllVsMPro[1]->Fill(dMultiplicityBin,two2n2n,mWeight2p);
    fIntFlowCorrelationsAllVsMPro[2]->Fill(dMultiplicityBin,two3n3n,mWeight2p);
    fIntFlowCorrelationsAllVsMPro[3]->Fill(dMultiplicityBin,two4n4n,mWeight2p);
-  }  
+  } 
+  if(fStoreControlHistograms)
+  {
+   fCorrelation2468VsMult[0]->Fill(dMultiplicityBin,two1n1n);
+  } 
  } // end of if(dMult>1)
  
  // 3-particle:
@@ -3723,13 +3999,13 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowCorrelationsEBE->SetBinContent(2,four1n1n1n1n); // <4>
   // Testing other multiplicity weights:
   Double_t mWeight4p = 0.;
-  if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+  if(fMultiplicityWeight->Contains("combinations"))
   {
    mWeight4p = dMult*(dMult-1.)*(dMult-2.)*(dMult-3.);
-  } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+  } else if(fMultiplicityWeight->Contains("unit"))
     {
      mWeight4p = 1.;    
-    } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+    } else if(fMultiplicityWeight->Contains("multiplicity"))
       {
        mWeight4p = dMult;           
       }      
@@ -3738,9 +4014,21 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowSquaredCorrelationsPro->Fill(1.5,four1n1n1n1n*four1n1n1n1n,mWeight4p);
   if(fCalculateCumulantsVsM)
   {
-   fIntFlowCorrelationsVsMPro[1]->Fill(dMultiplicityBin,four1n1n1n1n,mWeight4p);
-   fIntFlowSquaredCorrelationsVsMPro[1]->Fill(dMultiplicityBin,four1n1n1n1n*four1n1n1n1n,mWeight4p);
-  }   
+   if(fFillProfilesVsMUsingWeights)  
+   {
+    fIntFlowCorrelationsVsMPro[1]->Fill(dMultiplicityBin,four1n1n1n1n,mWeight4p);
+    fIntFlowSquaredCorrelationsVsMPro[1]->Fill(dMultiplicityBin,four1n1n1n1n*four1n1n1n1n,mWeight4p);
+   } else
+     {
+      fIntFlowCorrelationsVsMPro[1]->Fill(dMultiplicityBin,four1n1n1n1n);
+      fIntFlowSquaredCorrelationsVsMPro[1]->Fill(dMultiplicityBin,four1n1n1n1n*four1n1n1n1n);
+     }
+  } // end of if(fCalculateCumulantsVsM) 
+  if(fStoreControlHistograms)
+  {
+   fCorrelation2468VsMult[1]->Fill(dMultiplicityBin,four1n1n1n1n);
+   fCorrelationProduct2468VsMult[0]->Fill(dMultiplicityBin,two1n1n*four1n1n1n1n);
+  } 
  } // end of if(dMult>3)
 
  // 5-particle:
@@ -3870,13 +4158,13 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowCorrelationsEBE->SetBinContent(3,six1n1n1n1n1n1n); // <6>
   // Testing other multiplicity weights:
   Double_t mWeight6p = 0.;
-  if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+  if(fMultiplicityWeight->Contains("combinations"))
   {
    mWeight6p = dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.);
-  } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+  } else if(fMultiplicityWeight->Contains("unit"))
     {
      mWeight6p = 1.;    
-    } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+    } else if(fMultiplicityWeight->Contains("multiplicity"))
       {
        mWeight6p = dMult;           
       }
@@ -3885,9 +4173,20 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowSquaredCorrelationsPro->Fill(2.5,six1n1n1n1n1n1n*six1n1n1n1n1n1n,mWeight6p);
   if(fCalculateCumulantsVsM)
   {
-   fIntFlowCorrelationsVsMPro[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n,mWeight6p);
-   fIntFlowSquaredCorrelationsVsMPro[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n*six1n1n1n1n1n1n,mWeight6p);
-  }    
+   if(fFillProfilesVsMUsingWeights)  
+   {
+    fIntFlowCorrelationsVsMPro[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n,mWeight6p);
+    fIntFlowSquaredCorrelationsVsMPro[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n*six1n1n1n1n1n1n,mWeight6p);
+   } else
+     {
+      fIntFlowCorrelationsVsMPro[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n);
+      fIntFlowSquaredCorrelationsVsMPro[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n*six1n1n1n1n1n1n);
+     }
+  } // end of if(fCalculateCumulantsVsM)
+  if(fStoreControlHistograms)
+  {
+   fCorrelation2468VsMult[2]->Fill(dMultiplicityBin,six1n1n1n1n1n1n);
+  }     
  } // end of if(dMult>5)
  
  // 7-particle:
@@ -3965,13 +4264,13 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowCorrelationsEBE->SetBinContent(4,eight1n1n1n1n1n1n1n1n); // <8>
   // Testing other multiplicity weights:
   Double_t mWeight8p = 0.;
-  if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+  if(fMultiplicityWeight->Contains("combinations"))
   {
    mWeight8p = dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.);
-  } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+  } else if(fMultiplicityWeight->Contains("unit"))
     {
      mWeight8p = 1.;    
-    } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+    } else if(fMultiplicityWeight->Contains("multiplicity"))
       {
        mWeight8p = dMult;           
       }        
@@ -3980,9 +4279,20 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   fIntFlowSquaredCorrelationsPro->Fill(3.5,eight1n1n1n1n1n1n1n1n*eight1n1n1n1n1n1n1n1n,mWeight8p);  
   if(fCalculateCumulantsVsM)
   {
-   fIntFlowCorrelationsVsMPro[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n,mWeight8p);
-   fIntFlowSquaredCorrelationsVsMPro[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n*eight1n1n1n1n1n1n1n1n,mWeight8p);
-  }    
+   if(fFillProfilesVsMUsingWeights)  
+   {
+    fIntFlowCorrelationsVsMPro[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n,mWeight8p);
+    fIntFlowSquaredCorrelationsVsMPro[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n*eight1n1n1n1n1n1n1n1n,mWeight8p);
+   } else
+     {
+      fIntFlowCorrelationsVsMPro[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n);
+      fIntFlowSquaredCorrelationsVsMPro[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n*eight1n1n1n1n1n1n1n1n);
+     } 
+  } // end of if(fCalculateCumulantsVsM)
+  if(fStoreControlHistograms)
+  {
+   fCorrelation2468VsMult[3]->Fill(dMultiplicityBin,eight1n1n1n1n1n1n1n1n);
+  } // end of if(fStoreControlHistograms)  
  } // end of if(dMult>7) 
  
  // EXTRA correlations for v3{5} study:
@@ -4380,6 +4690,19 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
   }          
  } // end of if(dMult>5.)
  
+ // |Qn|^2/M, |Q2n|^2/M, |Qn|^4/(M(2M-1)), Re[Q2nQn^*Qn^*]/M, ... vs multiplicity (#RPs, #POIs or external):
+ if(fUseQvectorTerms)
+ {
+  Double_t dM = dMultiplicityBin-0.5;
+  if(dM>1.) // TBI re-think this if statement 
+  {
+   fQvectorTermsVsMult[0]->Fill(dMultiplicityBin,(pow(dReQ1n,2.)+pow(dImQ1n,2.))/dM);
+   fQvectorTermsVsMult[1]->Fill(dMultiplicityBin,(pow(dReQ2n,2.)+pow(dImQ2n,2.))/dM);
+   fQvectorTermsVsMult[2]->Fill(dMultiplicityBin,(pow(pow(dReQ1n,2.)+pow(dImQ1n,2.),2.))/(dM*(2.*dM-1.)));
+   fQvectorTermsVsMult[3]->Fill(dMultiplicityBin,reQ2nQ1nstarQ1nstar/pow(dM,1.5)); // TBI a bit of heuristic inserted here, re-think the rescaling factor 
+  } // end of if(dM>1.) // TBI re-think this if statement 
+ } // end of if(fUseQvectorTerms)
+
 } // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
 
 //=====================================================================================================
@@ -4611,7 +4934,7 @@ void AliFlowAnalysisWithQCumulants::CalculateMixedHarmonics()
  Double_t d6pMultiplicityWeight = 0.; // weight for <6>_{...} to get <<6>>_{...}
  Double_t d7pMultiplicityWeight = 0.; // weight for <7>_{...} to get <<7>>_{...}
  Double_t d8pMultiplicityWeight = 0.; // weight for <8>_{...} to get <<8>>_{...}
- if(!strcmp(fMultiplicityWeight->Data(),"combinations")) // default multiplicity weight
+ if(fMultiplicityWeight->Contains("combinations")) // default multiplicity weight
  {
   d2pMultiplicityWeight = dMult*(dMult-1.);
   d3pMultiplicityWeight = dMult*(dMult-1.)*(dMult-2.);
@@ -4620,7 +4943,7 @@ void AliFlowAnalysisWithQCumulants::CalculateMixedHarmonics()
   d6pMultiplicityWeight = dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.);
   d7pMultiplicityWeight = dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.);
   d8pMultiplicityWeight = dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.);
- } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+ } else if(fMultiplicityWeight->Contains("unit"))
    {
     d2pMultiplicityWeight = 1.;
     d3pMultiplicityWeight = 1.;
@@ -4629,7 +4952,7 @@ void AliFlowAnalysisWithQCumulants::CalculateMixedHarmonics()
     d6pMultiplicityWeight = 1.;
     d7pMultiplicityWeight = 1.;
     d8pMultiplicityWeight = 1.;
-   } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+   } else if(fMultiplicityWeight->Contains("multiplicity"))
      {
       d2pMultiplicityWeight = dMult;
       d3pMultiplicityWeight = dMult;
@@ -12649,19 +12972,19 @@ void AliFlowAnalysisWithQCumulants::StorePhiDistributionForOneEvent(AliFlowEvent
 void AliFlowAnalysisWithQCumulants::CalculateIntFlowProductOfCorrelations()
 {
  // Calculate averages of products of correlations for integrated flow.
- // multiplicity:
- Double_t dMult = (*fSpk)(0,0);
+
  // Multiplicity bin of an event (relevant for all histos vs M): 
  Double_t dMultiplicityBin = 0.;
- if(!fMultiplicityIsRefMultiplicity)
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
  {
-  dMultiplicityBin = dMult+0.5;
- } else
+  dMultiplicityBin = fNumberOfRPsEBE+0.5;
+ } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
    {
     dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
-   }
+   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
+     {
+      dMultiplicityBin = fNumberOfPOIsEBE+0.5;
+     }
 
  Int_t counter = 0;
  
@@ -12677,11 +13000,19 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowProductOfCorrelations()
    // products versus multiplicity:  // [0=<<2><4>>,1=<<2><6>>,2=<<2><8>>,3=<<4><6>>,4=<<4><8>>,5=<<6><8>>]
    if(fCalculateCumulantsVsM)
    {
-    fIntFlowProductOfCorrelationsVsMPro[counter]->Fill(dMultiplicityBin, // to be improved: dMult => sum of weights ?
-                                                       fIntFlowCorrelationsEBE->GetBinContent(ci1)*
-                                                       fIntFlowCorrelationsEBE->GetBinContent(ci2),
-                                                       fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci1)*
-                                                       fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci2));
+    if(fFillProfilesVsMUsingWeights) 
+    {
+     fIntFlowProductOfCorrelationsVsMPro[counter]->Fill(dMultiplicityBin, // to be improved: dMult => sum of weights ?
+                                                        fIntFlowCorrelationsEBE->GetBinContent(ci1)*
+                                                        fIntFlowCorrelationsEBE->GetBinContent(ci2),
+                                                        fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci1)*
+                                                        fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci2));
+    } else
+      {
+       fIntFlowProductOfCorrelationsVsMPro[counter]->Fill(dMultiplicityBin, // to be improved: dMult => sum of weights ?
+                                                          fIntFlowCorrelationsEBE->GetBinContent(ci1)*
+                                                          fIntFlowCorrelationsEBE->GetBinContent(ci2));
+      }
    } // end of if(fCalculateCumulantsVsM)
    counter++;                                                                                                                        
   }
@@ -12961,7 +13292,7 @@ void AliFlowAnalysisWithQCumulants::CalculateCovariancesIntFlow()
   {
    for(Int_t c2=c1+1;c2<4;c2++)
    {
-    productOfCorrelationsVsM[c1][c2] = fIntFlowProductOfCorrelationsVsMPro[productOfCorrelationsLabelVsM-1]->GetBinContent(b);
+    productOfCorrelationsVsM[c1][c2] = fIntFlowProductOfCorrelationsVsMPro[productOfCorrelationsLabelVsM-1]->GetBinContent(b); 
     if(TMath::Abs(fIntFlowSumOfEventWeightsVsM[c1][0]->GetBinContent(b)) > 1.e-44 && TMath::Abs(fIntFlowSumOfEventWeightsVsM[c2][0]->GetBinContent(b)) > 1.e-44)
     {
      denominatorVsM[c1][c2] = 1.-(fIntFlowSumOfProductOfEventWeightsVsM[sumOfProductOfEventWeightsLabel1VsM-1]->GetBinContent(b))
@@ -13964,33 +14295,123 @@ void AliFlowAnalysisWithQCumulants::CalculateCumulantsIntFlow()
  
  // Versus multiplicity: 
  if(!fCalculateCumulantsVsM){return;}
- Int_t nBins = fIntFlowCorrelationsVsMPro[0]->GetNbinsX(); // to be improved (hardwired 0) 
+ Int_t nBins = fIntFlowCorrelationsVsMPro[0]->GetNbinsX(); // tbi (hardwired 0) 
  Double_t value[4] = {0.}; // QCs vs M
  Double_t error[4] = {0.}; // error of QCs vs M
  Double_t dSum1[4] = {0.}; // sum value_i/(error_i)^2
  Double_t dSum2[4] = {0.}; // sum 1/(error_i)^2
+ two = 0.; // TBI can be done safer, lines below included 
+ twoError = 0.;
+ four = 0.;
+ fourError = 0.;
+ six = 0.;
+ sixError = 0.;
+ eight = 0.;
+ eightError = 0.;
+ wCov24 = 0.;
+ wCov26 = 0.; 
+ wCov28 = 0.; 
+ wCov46 = 0.; 
+ wCov48 = 0.; 
+ wCov68 = 0.; 
  for(Int_t b=1;b<=nBins;b++)
  {
   // Correlations:
-  two = fIntFlowCorrelationsVsMHist[0]->GetBinContent(b); // <<2>> 
-  four = fIntFlowCorrelationsVsMHist[1]->GetBinContent(b); // <<4>>  
-  six = fIntFlowCorrelationsVsMHist[2]->GetBinContent(b); // <<6>> 
-  eight = fIntFlowCorrelationsVsMHist[3]->GetBinContent(b); // <<8>>  
+  if(!fUse2DHistograms)
+  {
+   if(!fUseQvectorTerms)
+   {
+
+     //cout<<"TProfile"<<endl;
+
+    two = fIntFlowCorrelationsVsMHist[0]->GetBinContent(b); // <<2>>  
+    four = fIntFlowCorrelationsVsMHist[1]->GetBinContent(b); // <<4>>  
+    six = fIntFlowCorrelationsVsMHist[2]->GetBinContent(b); // <<6>> 
+    eight = fIntFlowCorrelationsVsMHist[3]->GetBinContent(b); // <<8>>  
+   } // end of if(!fUseQvectorTerms) 
+  } else
+    {
+     // TH2D:
+
+     cout<<"TH2D"<<endl;
+
+     two = fCorrelation2468VsMult[0]->ProjectionY("2",b,b)->GetMean(); // <<2>>  
+     four = fCorrelation2468VsMult[1]->ProjectionY("4",b,b)->GetMean(); // <<4>>
+     six = fCorrelation2468VsMult[2]->ProjectionY("6",b,b)->GetMean(); // <<6>>
+     eight = fCorrelation2468VsMult[3]->ProjectionY("8",b,b)->GetMean(); // <<8>>
+    } 
+  if(fUseQvectorTerms)
+  {
+   Double_t dM = fIntFlowCorrelationsVsMPro[0]->GetBinLowEdge(b);
+   if(dM>3.) // TBI re-think this if statement
+   {
+
+    cout<<"Q-vector terms"<<endl;
+
+    two = (fQvectorTermsVsMult[0]->ProjectionY("qvt0a",b,b)->GetMean()-1.)/(dM-1.);    
+    Double_t dTerm1 = (dM*(2.*dM-1.))*fQvectorTermsVsMult[2]->ProjectionY("qvt2",b,b)->GetMean();
+    Double_t dTerm2 = dM*fQvectorTermsVsMult[1]->ProjectionY("qvt1",b,b)->GetMean();
+    Double_t dTerm3 = -2.*pow(dM,1.5)*fQvectorTermsVsMult[3]->ProjectionY("qvt3",b,b)->GetMean();
+    Double_t dTerm4 = -4.*(dM-2.)*dM*fQvectorTermsVsMult[0]->ProjectionY("qvt0b",b,b)->GetMean();
+    Double_t dTerm5 = 2.*dM*(dM-3.); 
+    four = (dTerm1+dTerm2+dTerm3+dTerm4+dTerm5)/(dM*(dM-1.)*(dM-2.)*(dM-3.));
+   } // end of if(dM>3.) // TBI rethink this if statement
+  } // end of else if(fUseQvectorTerms) 
   // Statistical errors of average 2-, 4-, 6- and 8-particle azimuthal correlations:
-  twoError = fIntFlowCorrelationsVsMHist[0]->GetBinError(b); // statistical error of <2>  
-  fourError = fIntFlowCorrelationsVsMHist[1]->GetBinError(b); // statistical error of <4>   
-  sixError = fIntFlowCorrelationsVsMHist[2]->GetBinError(b); // statistical error of <6> 
-  eightError = fIntFlowCorrelationsVsMHist[3]->GetBinError(b); // statistical error of <8> 
+  if(!fUse2DHistograms)
+  {
+   if(!fUseQvectorTerms)
+   {
+    twoError = fIntFlowCorrelationsVsMHist[0]->GetBinError(b); // <2> error  
+    fourError = fIntFlowCorrelationsVsMHist[1]->GetBinError(b); // <4> error    
+    sixError = fIntFlowCorrelationsVsMHist[2]->GetBinError(b); // <6> error  
+    eightError = fIntFlowCorrelationsVsMHist[3]->GetBinError(b); // <8> error  
+   }
+  } else
+    {
+     // TH2D:
+     twoError = fCorrelation2468VsMult[0]->ProjectionY("2",b,b)->GetMeanError(); // <2> error
+     fourError = fCorrelation2468VsMult[1]->ProjectionY("4",b,b)->GetMeanError(); // <4> error 
+     sixError = fCorrelation2468VsMult[2]->ProjectionY("6",b,b)->GetMeanError(); // <6> error 
+     eightError = fCorrelation2468VsMult[3]->ProjectionY("8",b,b)->GetMeanError(); // <8> error 
+    } // end of else
+  if(fUseQvectorTerms)
+  {
+   Double_t dM = fIntFlowCorrelationsVsMPro[0]->GetBinLowEdge(b);
+   if(dM>3.) // TBI re-think this if statement
+   {
+    twoError = (fQvectorTermsVsMult[0]->ProjectionY("qvt0a",b,b)->GetMeanError())/(dM-1.); 
+   }
+  } // end of if(fUseQvectorTerms)
+
   // Covariances (multiplied by prefactor depending on weights - see comments in CalculateCovariancesIntFlow()):
   if(!fForgetAboutCovariances)
   {
-   wCov24 = fIntFlowCovariancesVsM[0]->GetBinContent(b); // Cov(<2>,<4>) * prefactor(w_<2>,w_<4>)
-   wCov26 = fIntFlowCovariancesVsM[1]->GetBinContent(b); // Cov(<2>,<6>) * prefactor(w_<2>,w_<6>)
-   wCov28 = fIntFlowCovariancesVsM[2]->GetBinContent(b); // Cov(<2>,<8>) * prefactor(w_<2>,w_<8>)
-   wCov46 = fIntFlowCovariancesVsM[3]->GetBinContent(b); // Cov(<4>,<6>) * prefactor(w_<4>,w_<6>)
-   wCov48 = fIntFlowCovariancesVsM[4]->GetBinContent(b); // Cov(<4>,<8>) * prefactor(w_<4>,w_<8>)
-   wCov68 = fIntFlowCovariancesVsM[5]->GetBinContent(b); // Cov(<6>,<8>) * prefactor(w_<6>,w_<8>) 
-  }
+   if(!fUse2DHistograms)
+   {
+    if(!fUseQvectorTerms)
+    {
+     wCov24 = fIntFlowCovariancesVsM[0]->GetBinContent(b); // Cov(<2>,<4>) * prefactor(w_<2>,w_<4>)
+     wCov26 = fIntFlowCovariancesVsM[1]->GetBinContent(b); // Cov(<2>,<6>) * prefactor(w_<2>,w_<6>)
+     wCov28 = fIntFlowCovariancesVsM[2]->GetBinContent(b); // Cov(<2>,<8>) * prefactor(w_<2>,w_<8>)
+     wCov46 = fIntFlowCovariancesVsM[3]->GetBinContent(b); // Cov(<4>,<6>) * prefactor(w_<4>,w_<6>)
+     wCov48 = fIntFlowCovariancesVsM[4]->GetBinContent(b); // Cov(<4>,<8>) * prefactor(w_<4>,w_<8>)
+     wCov68 = fIntFlowCovariancesVsM[5]->GetBinContent(b); // Cov(<6>,<8>) * prefactor(w_<6>,w_<8>) 
+    }
+   } else
+     {
+      // TH2D:
+       Int_t nEntries = (Int_t)(fCorrelationProduct2468VsMult[0]->ProjectionY("Cov(2,4)",b,b)->GetEntries());
+      if(nEntries>0.)
+      {      
+       wCov24 = (fCorrelationProduct2468VsMult[0]->ProjectionY("Cov(2,4)",b,b)->GetMean()
+              - (fCorrelation2468VsMult[0]->ProjectionY("2cov",b,b)->GetMean()) 
+              * (fCorrelation2468VsMult[1]->ProjectionY("4cov",b,b)->GetMean()))
+              / nEntries; // w*Cov(<2>,<4>)
+      } // end of if(nEntries>0. && dn2pCombinations>0. && dn4pCombinations>0.)
+      // ... TBI add other covariances as well 
+     }
+  } // end of if(!fForgetAboutCovariances)
   // Q-cumulants: 
   qc2 = 0.; // QC{2}
   qc4 = 0.; // QC{4}
@@ -14014,7 +14435,8 @@ void AliFlowAnalysisWithQCumulants::CalculateCumulantsIntFlow()
   qc2Error = twoError;                                             
   // Statistical error of QC{4}:              
   qc4ErrorSquared = 16.*pow(two,2.)*pow(twoError,2)+pow(fourError,2.)
-                  - 8.*two*wCov24;                     
+                  - 8.*two*wCov24;  
+  if(fUseQvectorTerms){qc4ErrorSquared = 0.01*pow(qc4,2.);} // TBI for the time being use bootstrap and ignore this error completely, this method needs an urgent clean-up, as the code is getting fragile here :'(                   
   if(qc4ErrorSquared>0.)
   {
    qc4Error = pow(qc4ErrorSquared,0.5);
@@ -15053,7 +15475,12 @@ void AliFlowAnalysisWithQCumulants::InitializeArraysForDistributions()
  fMaxValueOfCorrelation[2] = 0.0000006; // <6>_max  
  fMinValueOfCorrelation[3] = -0.000000006; // <8>_min 
  fMaxValueOfCorrelation[3] = 0.000000003; // <8>_max 
+
+ // c) Initialize default min and max values of correlation products:
+ //    (Remark: The default values bellow were chosen for v2=5% and M=500)
+ fMinValueOfCorrelationProduct[0] = -0.01; // <2><4>_min 
+ fMaxValueOfCorrelationProduct[0] = 0.04; // <2><4>_max 
+
 } // end of void AliFlowAnalysisWithQCumulants::InitializeArraysForDistributions()
 
 //=======================================================================================================================
@@ -15076,26 +15503,27 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForDistributions()
  // a) Book profile to hold all flags for distributions of correlations;
  // b) Book all histograms to hold distributions of correlations.
  
- TString correlationIndex[4] = {"<2>","<4>","<6>","<8>"}; // to be improved (should I promote this to data members?)
+ TString correlationIndex[4] = {"#LT2#GT","#LT4#GT","#LT6#GT","#LT8#GT"}; // TBI (should I promote this to data members?)
   
  // a) Book profile to hold all flags for distributions of correlations:
  TString distributionsFlagsName = "fDistributionsFlags";
  distributionsFlagsName += fAnalysisLabel->Data();
- fDistributionsFlags = new TProfile(distributionsFlagsName.Data(),"Flags for Distributions of Correlations",9,0,9);
+ fDistributionsFlags = new TProfile(distributionsFlagsName.Data(),"Flags for Distributions of Correlations",10,0,10);
  fDistributionsFlags->SetTickLength(-0.01,"Y");
  fDistributionsFlags->SetMarkerStyle(25);
  fDistributionsFlags->SetLabelSize(0.05);
  fDistributionsFlags->SetLabelOffset(0.02,"Y");
  fDistributionsFlags->SetStats(kFALSE);
  fDistributionsFlags->GetXaxis()->SetBinLabel(1,"Store or not?");
- fDistributionsFlags->GetXaxis()->SetBinLabel(2,"<2>_{min}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(3,"<2>_{max}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(4,"<4>_{min}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(5,"<4>_{max}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(6,"<6>_{min}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(7,"<6>_{max}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(8,"<8>_{min}");
- fDistributionsFlags->GetXaxis()->SetBinLabel(9,"<8>_{max}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(2,"#LT2#GT_{min}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(3,"#LT2#GT_{max}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(4,"#LT4#GT_{min}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(5,"#LT4#GT_{max}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(6,"#LT6#GT_{min}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(7,"#LT6#GT_{max}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(8,"#LT8#GT_{min}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(9,"#LT8#GT_{max}");
+ fDistributionsFlags->GetXaxis()->SetBinLabel(10,"fnBinsForCorrelations");
  fDistributionsList->Add(fDistributionsFlags);
  
  // b) Book all histograms to hold distributions of correlations.
@@ -15105,7 +15533,7 @@ void AliFlowAnalysisWithQCumulants::BookEverythingForDistributions()
   distributionsName += fAnalysisLabel->Data();
   for(Int_t di=0;di<4;di++) // distribution index
   {
-   fDistributions[di] = new TH1D(Form("Distribution of %s",correlationIndex[di].Data()),Form("Distribution of %s",correlationIndex[di].Data()),10000,fMinValueOfCorrelation[di],fMaxValueOfCorrelation[di]); 
+   fDistributions[di] = new TH1D(Form("Distribution of %s",correlationIndex[di].Data()),Form("Distribution of %s",correlationIndex[di].Data()),fnBinsForCorrelations,fMinValueOfCorrelation[di],fMaxValueOfCorrelation[di]); 
    fDistributions[di]->SetXTitle(correlationIndex[di].Data());
    fDistributionsList->Add(fDistributions[di]);
   } // end of for(Int_t di=0;di<4;di++) // distribution index
@@ -15149,6 +15577,7 @@ void AliFlowAnalysisWithQCumulants::StoreFlagsForDistributions()
   fDistributionsFlags->Fill(1.5+2.*(Double_t)di,fMinValueOfCorrelation[di]);
   fDistributionsFlags->Fill(2.5+2.*(Double_t)di,fMaxValueOfCorrelation[di]);
  }
+ fDistributionsFlags->Fill(9.5,fnBinsForCorrelations);
      
 } // end of void AliFlowAnalysisWithQCumulants::StoreFlagsForDistributions()
 
@@ -15192,7 +15621,8 @@ void AliFlowAnalysisWithQCumulants::BookAndNestAllLists()
  //  e) Book and nest list for various unclassified objects; 
  //  f) Book and nest list for other differential correlators;
  //  g) Book and nest list for nested loops;
- //  h) Book and nest lists for mixed harmonics.
+ //  h) Book and nest lists for mixed harmonics;
+ //  i) Book and nest lists for control histograms. 
  
  // a) Book and nest all lists for integrated flow:
  //  Base list for integrated flow:
@@ -15276,6 +15706,13 @@ void AliFlowAnalysisWithQCumulants::BookAndNestAllLists()
  fMixedHarmonicsErrorPropagation->SetOwner(kTRUE);
  if(fCalculateMixedHarmonics){fMixedHarmonicsList->Add(fMixedHarmonicsErrorPropagation);}
 
+ // i) Book and nest lists for control histograms:
+ //  Base list for mixed harmonics:
+ fControlHistogramsList = new TList();
+ fControlHistogramsList->SetName("Control Histograms");
+ fControlHistogramsList->SetOwner(kTRUE);
+ fHistList->Add(fControlHistogramsList);
+
 } // end of void AliFlowAnalysisWithQCumulants::BookAndNestAllLists()
 
 //=======================================================================================================================
@@ -15554,18 +15991,19 @@ void AliFlowAnalysisWithQCumulants::CommonConstants(TString method)
 
 void AliFlowAnalysisWithQCumulants::CrossCheckSettings()
 {
- // a) Cross check if the choice for multiplicity weights make sense.
- // a) Cross check if the choice for multiplicity weights make sense:
- if(strcmp(fMultiplicityWeight->Data(),"combinations") && 
-    strcmp(fMultiplicityWeight->Data(),"unit") &&
-    strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+ // a) Cross-check if the choice for multiplicity weights make sense;
+ // b) Cross-check if the choice for multiplicity itself make sense.
+
+ // a) Cross-check if the choice for multiplicity weights make sense:
+ if((!fMultiplicityWeight->Contains("combinations")) && 
+    (!fMultiplicityWeight->Contains("unit")) &&
+    (!fMultiplicityWeight->Contains("multiplicity")) )
  {
   cout<<"WARNING (QC): Multiplicity weight can be either \"combinations\", \"unit\""<<endl;
   cout<<"              or \"multiplicity\". Certainly not \""<<fMultiplicityWeight->Data()<<"\"."<<endl;
   exit(0);
  }   
+
 } // end of void AliFlowAnalysisWithQCumulants::CrossCheckSettings()
 
 //=======================================================================================================================
@@ -15574,18 +16012,20 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfEventWeights()
 {
  // Calculate sum of linear and quadratic event weights for correlations.
  
- // multiplicity:
- Double_t dMult = (*fSpk)(0,0);
+ // TBI re-think what is the right multiplicity when particle weights are used! 
 
  // Multiplicity bin of an event (relevant for all histos vs M): 
  Double_t dMultiplicityBin = 0.;
- if(!fMultiplicityIsRefMultiplicity)
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
  {
-  dMultiplicityBin = dMult+0.5;
- } else
+  dMultiplicityBin = fNumberOfRPsEBE+0.5;
+ } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
    {
-    dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
-   }
+     dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
+   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
+     {
+      dMultiplicityBin = fNumberOfPOIsEBE+0.5;
+     }
                         
  for(Int_t p=0;p<2;p++) // power-1
  {
@@ -15594,10 +16034,16 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfEventWeights()
    fIntFlowSumOfEventWeights[p]->Fill(ci+0.5,pow(fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci+1),p+1)); 
    if(fCalculateCumulantsVsM)
    {
-    fIntFlowSumOfEventWeightsVsM[ci][p]->Fill(dMultiplicityBin,pow(fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci+1),p+1)); // to be improved: dMult => sum of weights?
-   }
-  }
- }
+    if(fFillProfilesVsMUsingWeights)
+    {
+     fIntFlowSumOfEventWeightsVsM[ci][p]->Fill(dMultiplicityBin,pow(fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci+1),p+1)); // to be improved: dMult => sum of weights?
+    } else
+      {
+       fIntFlowSumOfEventWeightsVsM[ci][p]->Fill(dMultiplicityBin); // to be improved: dMult => sum of weights?
+      }
+   } // end of if(fCalculateCumulantsVsM)
+  } // end of for(Int_t ci=0;ci<4;ci++) // correlation index
+ } // end of for(Int_t p=0;p<2;p++) // power-1
   
 } // end of void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfEventWeights()
 
@@ -15626,18 +16072,20 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfProductOfEventWeights()
 {
  // Calculate sum of product of event weights for correlations.
   
- // multiplicity:
- Double_t dMult = (*fSpk)(0,0);
+ // TBI re-think what is the right multiplicity when particle weights are used!
 
  // Multiplicity bin of an event (relevant for all histos vs M): 
  Double_t dMultiplicityBin = 0.;
- if(!fMultiplicityIsRefMultiplicity)
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP)
  {
-  dMultiplicityBin = dMult+0.5;
- } else
+  dMultiplicityBin = fNumberOfRPsEBE+0.5;
+ } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
    {
     dMultiplicityBin = fReferenceMultiplicityEBE+0.5;
-   }
+   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
+     {
+      dMultiplicityBin = fNumberOfPOIsEBE+0.5;
+     }
   
  Int_t counter = 0;
  
@@ -15650,13 +16098,19 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfProductOfEventWeights()
                                             fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci2));
    if(fCalculateCumulantsVsM)
    {                                                                                    
-    fIntFlowSumOfProductOfEventWeightsVsM[counter]->Fill(dMultiplicityBin, // to be improved: dMult => sum of weights?
-                                                         fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci1)*
-                                                         fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci2));
+    if(fFillProfilesVsMUsingWeights)
+    {
+     fIntFlowSumOfProductOfEventWeightsVsM[counter]->Fill(dMultiplicityBin, // to be improved: dMult => sum of weights?
+                                                          fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci1)*
+                                                          fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci2));
+    } else
+      {
+       fIntFlowSumOfProductOfEventWeightsVsM[counter]->Fill(dMultiplicityBin);
+      }
    } // end of if(fCalculateCumulantsVsM)
    counter++;                                         
-  }
- }
+  } // end of for(Int_t ci2=ci1+1;ci2<=4;ci2++)
+ } // end of for(Int_t ci1=1;ci1<4;ci1++)
 
 } // end of void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfProductOfEventWeights()
 
@@ -15879,10 +16333,10 @@ void AliFlowAnalysisWithQCumulants::CalculateDiffFlowCorrelations(TString type,
     two1n1nPtEta = (p1n0kRe*dReQ1n+p1n0kIm*dImQ1n-mq)
                  / (mp*dMult-mq);
     // determine multiplicity weight:
-    if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+    if(fMultiplicityWeight->Contains("combinations"))
     {
      mWeight2pPrime = mp*dMult-mq;
-    } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+    } else if(fMultiplicityWeight->Contains("unit"))
       {
        mWeight2pPrime = 1.;    
       } 
@@ -15929,10 +16383,10 @@ void AliFlowAnalysisWithQCumulants::CalculateDiffFlowCorrelations(TString type,
                       / ((mp-mq)*dMult*(dMult-1.)*(dMult-2.)
                           + mq*(dMult-1.)*(dMult-2.)*(dMult-3.)); 
     // determine multiplicity weight:
-    if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+    if(fMultiplicityWeight->Contains("combinations"))
     {
      mWeight4pPrime = (mp-mq)*dMult*(dMult-1.)*(dMult-2.) + mq*(dMult-1.)*(dMult-2.)*(dMult-3.);
-    } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+    } else if(fMultiplicityWeight->Contains("unit"))
       {
        mWeight4pPrime = 1.;    
       }     
@@ -16106,10 +16560,10 @@ void AliFlowAnalysisWithQCumulants::CalculateOtherDiffCorrelators(TString type,
                + 2.*mq)
                / ((mp*dMult-2.*mq)*(dMult-1.));
     // determine multiplicity weight:
-    if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+    if(fMultiplicityWeight->Contains("combinations"))
     {
      mWeightTaeneyYan = (mp*dMult-2.*mq)*(dMult-1.);
-    } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+    } else if(fMultiplicityWeight->Contains("unit"))
       {
        mWeightTaeneyYan = 1.;    
       } 
@@ -16231,10 +16685,10 @@ void AliFlowAnalysisWithQCumulants::Calculate2DDiffFlowCorrelations(TString type
     two1n1nPtEta = (p1n0kRe*dReQ1n+p1n0kIm*dImQ1n-mq)
                  / (mp*dMult-mq);
     // Determine multiplicity weight:
-    if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+    if(fMultiplicityWeight->Contains("combinations"))
     {
      mWeight2pPrime = mp*dMult-mq;
-    } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+    } else if(fMultiplicityWeight->Contains("unit"))
       {
        mWeight2pPrime = 1.;    
       } 
@@ -16263,10 +16717,10 @@ void AliFlowAnalysisWithQCumulants::Calculate2DDiffFlowCorrelations(TString type
                       / ((mp-mq)*dMult*(dMult-1.)*(dMult-2.)
                           + mq*(dMult-1.)*(dMult-2.)*(dMult-3.)); 
     // Determine multiplicity weight:
-    if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+    if(fMultiplicityWeight->Contains("combinations"))
     {
      mWeight4pPrime = (mp-mq)*dMult*(dMult-1.)*(dMult-2.) + mq*(dMult-1.)*(dMult-2.)*(dMult-3.);
-    } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+    } else if(fMultiplicityWeight->Contains("unit"))
       {
        mWeight4pPrime = 1.;    
       }     
@@ -17170,13 +17624,13 @@ void AliFlowAnalysisWithQCumulants::StoreIntFlowFlags()
  // particle weights used or not:
  fIntFlowFlags->Fill(0.5,(Int_t)fUsePhiWeights||fUsePtWeights||fUseEtaWeights||fUseTrackWeights);
  // which event weights were used:
- if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+ if(fMultiplicityWeight->Contains("combinations"))
  {
   fIntFlowFlags->Fill(1.5,0); // 0 = "combinations" (default)
- } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+ } else if(fMultiplicityWeight->Contains("unit"))
    {
     fIntFlowFlags->Fill(1.5,1); // 1 = "unit"   
-   } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+   } else if(fMultiplicityWeight->Contains("multiplicity"))
      {
       fIntFlowFlags->Fill(1.5,2); // 2 = "multiplicity"        
      } 
@@ -17193,7 +17647,22 @@ void AliFlowAnalysisWithQCumulants::StoreIntFlowFlags()
  fIntFlowFlags->Fill(12.5,(Int_t)fStorePhiDistributionForOneEvent); 
  fIntFlowFlags->Fill(13.5,(Int_t)fFillMultipleControlHistograms);  
  fIntFlowFlags->Fill(14.5,(Int_t)fCalculateAllCorrelationsVsM);  
- fIntFlowFlags->Fill(15.5,(Int_t)fMultiplicityIsRefMultiplicity);  
+ // which multiplicity was used:
+ if(fMultiplicityIs==AliFlowCommonConstants::kRP) // # of Reference Particles
+ {
+  fIntFlowFlags->Fill(15.5,0); // 0 = # of Reference Particles
+ } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal)
+   {
+    fIntFlowFlags->Fill(15.5,1); // 1 = ref. mult. from ESD
+   } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI)
+     {
+      fIntFlowFlags->Fill(15.5,2); // 2 = # of Particles of Interest
+     } 
+ fIntFlowFlags->Fill(16.5,(Int_t)fExactNoRPs); 
+ fIntFlowFlags->Fill(17.5,(Int_t)fUse2DHistograms); 
+ fIntFlowFlags->Fill(18.5,(Int_t)fFillProfilesVsMUsingWeights); 
+ fIntFlowFlags->Fill(19.5,(Int_t)fUseQvectorTerms); 
+
 } // end of void AliFlowAnalysisWithQCumulants::StoreIntFlowFlags()
 
 //=======================================================================================================================
@@ -17234,13 +17703,13 @@ void AliFlowAnalysisWithQCumulants::StoreMixedHarmonicsFlags()
  fMixedHarmonicsFlags->Fill(2.5,(Int_t)fCalculateMixedHarmonicsVsM);
  // Which multiplicity weight was used?:
 
- if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
+ if(fMultiplicityWeight->Contains("combinations"))
  {
   fMixedHarmonicsFlags->Fill(3.5,0); // 0 = "combinations" (default)
- } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
+ } else if(fMultiplicityWeight->Contains("unit"))
    {
     fMixedHarmonicsFlags->Fill(3.5,1); // 1 = "unit"   
-   } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
+   } else if(fMultiplicityWeight->Contains("multiplicity"))
      {
       fMixedHarmonicsFlags->Fill(3.5,2); // 2 = "multiplicity"        
      } 
@@ -17249,6 +17718,23 @@ void AliFlowAnalysisWithQCumulants::StoreMixedHarmonicsFlags()
 
 //=======================================================================================================================
 
+void AliFlowAnalysisWithQCumulants::StoreControlHistogramsFlags()
+{
+ // Store all flags for control histograms in profile fControlHistogramsFlags.
+
+ if(!fControlHistogramsFlags)
+ {
+  cout<<"WARNING: fControlHistogramsFlags is NULL in AFAWQC::SCHF() !!!!"<<endl;
+  exit(0);
+ } 
+
+ fControlHistogramsFlags->Fill(0.5,(Int_t)fStoreControlHistograms);
+ fControlHistogramsFlags->Fill(1.5,(Int_t)fUseQvectorTerms);
+
+} // end of void AliFlowAnalysisWithQCumulants::StoreControlHistogramsFlags()
+
+//=======================================================================================================================
+
 void AliFlowAnalysisWithQCumulants::GetPointersForCommonHistograms() 
 {
  // Access all pointers to common control and common result histograms and profiles.
@@ -19669,6 +20155,99 @@ void AliFlowAnalysisWithQCumulants::GetPointersForMixedHarmonicsHistograms()
 
 //=======================================================================================================================
 
+void AliFlowAnalysisWithQCumulants::GetPointersForControlHistograms()
+{
+ // Get pointers to all control histograms.
+   
+ // a) Get pointer to base list for control histograms;
+ // b) Get pointer to TProfile fControlHistogramsFlags holding all flags for control histograms;
+ // c) Get pointers to TH2D *fCorrelation2468VsMult[4], TH2D *fCorrelationProduct2468VsMult[1] and TH2D *fQvectorTermsVsMult[4].
+
+ // a) Get pointer to base list for control histograms:
+ TList *controlHistogramsList = dynamic_cast<TList*>(fHistList->FindObject("Control Histograms"));
+ if(controlHistogramsList) 
+ {
+  this->SetControlHistogramsList(controlHistogramsList);
+ } else
+   {
+    cout<<"WARNING: controlHistogramsList is NULL in AFAWQC::GPFMHH() !!!!"<<endl;
+    exit(0);
+   }
+
+ // b) Get pointer to TProfile fControlHistogramsFlags holding all flags for control histograms:
+ TString controlHistogramsFlagsName = "fControlHistogramsFlags";
+ controlHistogramsFlagsName += fAnalysisLabel->Data();
+ TProfile *controlHistogramsFlags = dynamic_cast<TProfile*>
+                                 (controlHistogramsList->FindObject(controlHistogramsFlagsName.Data()));
+ if(controlHistogramsFlags)
+ {
+  this->SetControlHistogramsFlags(controlHistogramsFlags);  
+  fStoreControlHistograms = (Bool_t)controlHistogramsFlags->GetBinContent(1); 
+  fUseQvectorTerms = (Bool_t)controlHistogramsFlags->GetBinContent(2); 
+ } else 
+   {
+    cout<<"WARNING: controlHistogramsFlags is NULL in AFAWQC::GPFMHH() !!!!"<<endl;
+    exit(0);
+   }
+  
+ if(!fStoreControlHistograms){return;}
+
+ // c) Get pointers to TH2D *fCorrelation2468VsMult[4], TH2D *fCorrelationProduct2468VsMult[1] and TH2D *fQvectorTermsVsMult[4]:
+ TString sCorrelation[4] = {"#LT2#GT","#LT4#GT","#LT6#GT","#LT8#GT"};
+ TString sCorrelation2468VsMultName = "fCorrelation2468VsMult";
+ sCorrelation2468VsMultName += fAnalysisLabel->Data();
+ for(Int_t ci=0;ci<4;ci++)
+ {
+  TH2D *hCorrelation2468VsMult = dynamic_cast<TH2D*>(controlHistogramsList->FindObject(Form("%s vs M",sCorrelation[ci].Data())));
+  if(hCorrelation2468VsMult) 
+  {
+   this->SetCorrelation2468VsMult(hCorrelation2468VsMult,ci);
+  } else 
+    {
+     cout<<"WARNING: hCorrelation2468VsMult is NULL in AFAWQC::GPFCH() !!!!"<<endl; 
+     cout<<"ci = "<<ci<<endl;
+     exit(0);
+    }                                   
+ } // end of for(Int_t ci=0;ci<4;ci++)
+ TString sCorrelationProduct[1] = {"#LT2#GT#LT4#GT"}; // TBI: add the other ones when needed first time
+ TString sCorrelationProduct2468VsMultName = "fCorrelationProduct2468VsMult";
+ sCorrelationProduct2468VsMultName += fAnalysisLabel->Data();
+ for(Int_t cpi=0;cpi<1;cpi++) // TBI: hardwired 1
+ {
+  TH2D *hCorrelationProduct2468VsMult = dynamic_cast<TH2D*>(controlHistogramsList->FindObject(Form("%s vs M",sCorrelationProduct[cpi].Data())));
+  if(hCorrelationProduct2468VsMult) 
+  {
+   this->SetCorrelationProduct2468VsMult(hCorrelationProduct2468VsMult,cpi);
+  } else 
+    {
+     cout<<"WARNING: hCorrelationProduct2468VsMult is NULL in AFAWQC::GPFCH() !!!!"<<endl; 
+     cout<<"cpi = "<<cpi<<endl;
+     exit(0);
+    }  
+ } // end of for(Int_t cpi=0;cpi<1;cpi++) // TBI: hardwired 1
+
+ if(!fUseQvectorTerms){return;}
+ 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}}"};
+ TString sQvectorTermsVsMultName = "fQvectorTermsVsMult";
+ sQvectorTermsVsMultName += fAnalysisLabel->Data();
+ for(Int_t qwti=0;qwti<4;qwti++) // TBI: hardwired 4
+ {
+  TH2D *hQvectorTermsVsMult = dynamic_cast<TH2D*>(controlHistogramsList->FindObject(Form("%s vs M",sQvectorTerms[qwti].Data())));
+  if(hQvectorTermsVsMult) 
+  {
+   this->SetQvectorTermsVsMult(hQvectorTermsVsMult,qwti);
+  } else 
+    {
+     cout<<"WARNING: hQvectorTermsVsMult is NULL in AFAWQC::GPFCH() !!!!"<<endl; 
+     cout<<"qwti = "<<qwti<<endl;
+     exit(0);
+    }  
+ } // end of for(Int_t qwti=0;qwti<1;qwti++) // TBI: hardwired 4
+
+} // end of void AliFlowAnalysisWithQCumulants::GetPointersForControlHistograms()
+
+//=======================================================================================================================
+
 void AliFlowAnalysisWithQCumulants::StoreHarmonic()
 {
  // Store flow harmonic in common control histograms.
@@ -19873,7 +20452,7 @@ void AliFlowAnalysisWithQCumulants::FillCommonControlHistograms(AliFlowEventSimp
 {
  // Fill common control histograms.
  
- Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of RPs (i.e. number of particles used to determine the reaction plane)
+ Int_t nRP = anEvent->GetNumberOfRPs(); // number of Reference Particles 
  fCommonHists->FillControlHistograms(anEvent); 
  if(fFillMultipleControlHistograms)
  {
@@ -19899,6 +20478,22 @@ void AliFlowAnalysisWithQCumulants::FillCommonControlHistograms(AliFlowEventSimp
 
 //=======================================================================================================================
 
+void AliFlowAnalysisWithQCumulants::FillControlHistograms(AliFlowEventSimple *anEvent)
+{
+ // Fill common control histograms.
+ Int_t nRPs = anEvent->GetNumberOfRPs(); // number of Reference Particles
+ Int_t nPOIs = anEvent->GetNumberOfPOIs(); // number of Particles Of Interest
+ Int_t nRefMult = anEvent->GetReferenceMultiplicity(); // reference multiplicity for current event (TBI: This call is not really needed here, use fReferenceMultiplicityEBE instead)
+
+ fCorrelationNoRPsVsRefMult->Fill(nRPs,nRefMult);
+ fCorrelationNoPOIsVsRefMult->Fill(nPOIs,nRefMult);
+ fCorrelationNoRPsVsNoPOIs->Fill(nRPs,nPOIs);
+} // end of void AliFlowAnalysisWithQCumulants::FillControlHistograms(AliFlowEventSimple *anEvent)
+
+//=======================================================================================================================
+
 void AliFlowAnalysisWithQCumulants::ResetEventByEventQuantities()
 {
  // Reset all event by event quantities.