PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx

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