PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx

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