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 strong parity violation.                     *
  26  *                                                        *
  27  * Author: Ante Bilandzic (abilandzic@gmail.com)          *
  28  *********************************************************/
  29
  30 #define AliFlowAnalysisWithMixedHarmonics_cxx
  31
  32 #include "Riostream.h"
  33 #include "AliFlowCommonConstants.h"
  34 #include "AliFlowCommonHist.h"
  35 #include "AliFlowCommonHistResults.h"
  36
  37 #include "TMath.h"
  38 #include "TFile.h"
  39 #include "TList.h"
  40 #include "TProfile.h"
  41 #include "TProfile2D.h"
  42
  43 #include "AliFlowEventSimple.h"
  44 #include "AliFlowTrackSimple.h"
  45 #include "AliFlowAnalysisWithMixedHarmonics.h"
  46
  47 class TH1;
  48 class TList;
  49 ClassImp(AliFlowAnalysisWithMixedHarmonics)
  50
  51 //================================================================================================================
  52
  53 AliFlowAnalysisWithMixedHarmonics::AliFlowAnalysisWithMixedHarmonics():
  54 fHistList(NULL),
  55 fHistListName(NULL),
  56 fAnalysisLabel(NULL),
  57 fAnalysisSettings(NULL),
  58 fCorrelatorInteger(2),
  59 fNoOfMultipicityBins(10),
  60 fMultipicityBinWidth(2),
  61 fMinMultiplicity(1),
  62 fCorrectForDetectorEffects(kTRUE),
  63 fCommonHists(NULL),
  64 fnBinsPhi(0),
  65 fPhiMin(0),
  66 fPhiMax(0),
  67 fPhiBinWidth(0),
  68 fnBinsPt(0),
  69 fPtMin(0),
  70 fPtMax(0),
  71 fPtBinWidth(0),
  72 fnBinsEta(0),
  73 fEtaMin(0),
  74 fEtaMax(0),
  75 fEtaBinWidth(0),
  76 fWeightsList(NULL),
  77 fUsePhiWeights(kFALSE),
  78 fUsePtWeights(kFALSE),
  79 fUseEtaWeights(kFALSE),
  80 fUseParticleWeights(NULL),
  81 fPhiWeights(NULL),
  82 fPtWeights(NULL),
  83 fEtaWeights(NULL),
  84 fReQnk(NULL),
  85 fImQnk(NULL),
  86 fSpk(NULL),
  87 f3pCorrelatorEBE(NULL),
  88 fNonIsotropicTermsEBE(NULL),
  89 fProfileList(NULL),
  90 f3pCorrelatorPro(NULL),
  91 fNonIsotropicTermsPro(NULL),
  92 f3pCorrelatorVsMPro(NULL),
  93 fNonIsotropicTermsVsMPro(NULL),
  94 fResultsList(NULL),
  95 f3pCorrelatorHist(NULL),
  96 fDetectorBiasHist(NULL),
  97 fNonIsotropicTermsHist(NULL),
  98 f3pCorrelatorVsMHist(NULL),
  99 fDetectorBiasVsMHist(NULL),
 100 fNonIsotropicTermsVsMHist(NULL)
 101 {
 102  // Constructor.
 103   // Base list to hold all output objects:
 104  fHistList = new TList();
 105  fHistListName = new TString("cobjMH");
 106  fHistList->SetName(fHistListName->Data());
 107  fHistList->SetOwner(kTRUE);
 108
 109  // List to hold histograms with phi, pt and eta weights:
 110  fWeightsList = new TList();
 111
 112  // List to hold all all-event profiles:
 113  fProfileList = new TList();
 114
 115  // List to hold objects with final results:
 116  fResultsList = new TList();
 117
 118 } // AliFlowAnalysisWithMixedHarmonics::AliFlowAnalysisWithMixedHarmonics()
 119
 120 //================================================================================================================
 121
 122 AliFlowAnalysisWithMixedHarmonics::~AliFlowAnalysisWithMixedHarmonics()
 123 {
 124  // Destructor.
 125
 126  delete fHistList;
 127
 128 } // end of AliFlowAnalysisWithMixedHarmonics::~AliFlowAnalysisWithMixedHarmonics()
 129
 130 //================================================================================================================
 131
 132 void AliFlowAnalysisWithMixedHarmonics::Init()
 133 {
 134  // Initialize and book all objects.
 135
 136  // a) Cross check if the user settings make sense before starting;
 137  // b) Access all common constants;
 138  // c) Book and nest all lists in the base list fHistList;
 139  // d) Book common control histograms;
 140  // e) Book all event-by-event quantities;
 141  // f) Book all all-event quantities;
 142  // g) Book and fill histograms to hold phi, pt and eta weights.
 143
 144  //save old value and prevent histograms from being added to directory
 145  //to avoid name clashes in case multiple analaysis objects are used
 146  //in an analysis
 147  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 148  TH1::AddDirectory(kFALSE);
 149
 150  this->CrossCheckSettings();
 151  this->AccessConstants();
 152  this->BookAndNestAllLists();
 153  this->BookProfileHoldingSettings();
 154  this->BookCommonHistograms();
 155  this->BookAllEventByEventQuantities();
 156  this->BookAllAllEventQuantities();
 157  this->BookAndFillWeightsHistograms();
 158
 159  TH1::AddDirectory(oldHistAddStatus);
 160 } // end of void AliFlowAnalysisWithMixedHarmonics::Init()
 161
 162 //================================================================================================================
 163
 164 void AliFlowAnalysisWithMixedHarmonics::Make(AliFlowEventSimple* anEvent)
 165 {
 166  // Running over data only in this method.
 167
 168  // a) Define local variables;
 169  // b) Fill common control histograms;
 170  // c) Loop over data and calculate e-b-e quantities Q_{n,k} and S_{p,k};
 171  // d) Calculate 3-p azimuthal correlator and non-isotropic terms in terms of Q_{n,k} and S_{p,k};
 172  // e) Reset all event-by-event quantities.
 173
 174  this->CheckPointersUsedInMake();
 175
 176  // a) Define local variables:
 177  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 178  Double_t dPt  = 0.; // transverse momentum
 179  Double_t dEta = 0.; // pseudorapidity
 180  Double_t wPhi = 1.; // phi weight
 181  Double_t wPt  = 1.; // pt weight
 182  Double_t wEta = 1.; // eta weight
 183  AliFlowTrackSimple *aftsTrack = NULL; // simple track
 184
 185  // b) Fill common control histograms.
 186  fCommonHists->FillControlHistograms(anEvent);
 187
 188  // c) Loop over data and calculate e-b-e quantities:
 189  Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
 190                                            // nRP   = # of particles used to determine the reaction plane ("Reference Particles");
 191                                            // nPOI  = # of particles of interest for a detailed flow analysis ("Particles of Interest");
 192                                            // rest  = # of particles which are not niether RPs nor POIs.
 193  // Start loop over data:
 194  for(Int_t i=0;i<nPrim;i++)
 195  {
 196   aftsTrack=anEvent->GetTrack(i);
 197   if(aftsTrack)
 198   {
 199    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue; // consider only tracks which are either RPs or POIs
 200    Int_t n = fCorrelatorInteger; // integer n in the correlator cos[n(2phi1-phi2-phi3)]
 201    if(aftsTrack->InRPSelection()) // checking RP condition:
 202    {
 203     dPhi = aftsTrack->Phi();
 204     dPt  = aftsTrack->Pt();
 205     dEta = aftsTrack->Eta();
 206     if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi-weight for this particle:
 207     {
 208      wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 209     }
 210     if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt-weight for this particle:
 211     {
 212      wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 213     }
 214     if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta-weight for this particle:
 215     {
 216      wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 217     }
 218     // Calculate Re[Q_{m*n,k}] and Im[Q_{m*n,k}], (m = 1,2 and k = 0,1,2,3) for this event:
 219     for(Int_t m=0;m<2;m++)
 220     {
 221      for(Int_t k=0;k<4;k++) // to be improved (what is maximum k that I need?)
 222      {
 223       (*fReQnk)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1)*n*dPhi);
 224       (*fImQnk)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1)*n*dPhi);
 225      }
 226     }
 227     // Calculate partially S_{p,k} for this event (final calculation of S_{p,k} follows after the loop over data bellow):
 228     for(Int_t p=0;p<4;p++) // to be improved (what is maximum p that I need?)
 229     {
 230      for(Int_t k=0;k<4;k++) // to be improved (what is maximum k that I need?)
 231      {
 232       (*fSpk)(p,k)+=pow(wPhi*wPt*wEta,k);
 233      }
 234     }
 235    } // end of if(aftsTrack->InRPSelection())
 236   } else // to if(aftsTrack)
 237     {
 238      cout<<endl;
 239      cout<<" WARNING (MH): No particle! (i.e. aftsTrack is a NULL pointer in AFAWMH::Make().)"<<endl;
 240      cout<<endl;
 241     }
 242  } // end of for(Int_t i=0;i<nPrim;i++)
 243
 244  // Calculate the final expressions for S_{p,k}:
 245  for(Int_t p=0;p<4;p++) // to be improved (what is maximum p that I need?)
 246  {
 247   for(Int_t k=0;k<4;k++) // to be improved (what is maximum  that I need?)
 248   {
 249    (*fSpk)(p,k)=pow((*fSpk)(p,k),p+1);
 250   }
 251  }
 252
 253  // d) Calculate 3-p azimuthal correlator in terms of Q_{n,k} and S_{p,k}:
 254  if(anEvent->GetEventNSelTracksRP() >= 3)
 255  {
 256   this->Calculate3pCorrelator();
 257  }
 258  if(anEvent->GetEventNSelTracksRP() >= 0) // to be improved (is this correct if statement?)
 259  {
 260   this->CalculateNonIsotropicTerms();
 261  }
 262
 263  // e) Reset all event-by-event quantities:
 264  this->ResetEventByEventQuantities();
 265
 266 } // end of AliFlowAnalysisWithMixedHarmonics::Make(AliFlowEventSimple* anEvent)
 267 //================================================================================================================
 268
 269 void AliFlowAnalysisWithMixedHarmonics::Finish()
 270 {
 271  // Calculate the final results.
 272
 273  // a) Check all pointers used in this method;
 274  // b) Access settings for analysis with mixed harmonics;
 275  // c) Calculate errors for non-isotropic terms;
 276  // d) Correct for detector effects;
 277  // e) Quantify bias to 3-p correlator coming from detector effects.
 278
 279  this->CheckPointersUsedInFinish();
 280  this->AccessSettings();
 281  this->FinalizeNonIsotropicTerms();
 282  this->CorrectForDetectorEffects();
 283  this->QuantifyBiasFromDetectorEffects();
 284
 285  cout<<endl;
 286  cout<<"*************************************"<<endl;
 287  cout<<"*************************************"<<endl;
 288  cout<<" flow estimates from mixed harmonics "<<endl;
 289  cout<<endl;
 290  if(fCorrelatorInteger!=1)
 291  {
 292   cout<< "  cos["<<fCorrelatorInteger<<"(2phi1-phi2-phi3)] = "<<f3pCorrelatorPro->GetBinContent(1)<<endl;
 293  } else
 294    {
 295     cout<< "  cos(2phi1-phi2-phi3) = "<<f3pCorrelatorPro->GetBinContent(1)<<endl;
 296    }
 297
 298  cout<<endl;
 299  cout<<"*************************************"<<endl;
 300  cout<<"*************************************"<<endl;
 301
 302 } // end of AliFlowAnalysisWithMixedHarmonics::Finish()
 303
 304 //================================================================================================================
 305
 306 void AliFlowAnalysisWithMixedHarmonics::GetOutputHistograms(TList *outputListHistos)
 307 {
 308  // Get pointers to all objects saved in the output file.
 309
 310  // a) Get pointers for common control histograms.
 311  if(outputListHistos)
 312  {
 313   this->SetHistList(outputListHistos);
 314   if(!fHistList)
 315   {
 316    cout<<endl;
 317    cout<<" WARNING (MH): fHistList is NULL in AFAWMH::GOH() !!!!"<<endl;
 318    cout<<endl;
 319    exit(0);
 320   }
 321   this->GetPointersForBaseHistograms();
 322   this->GetPointersForCommonHistograms();
 323   this->GetPointersForAllEventProfiles();
 324   this->GetPointersForResultsHistograms();
 325  } else
 326    {
 327     cout<<endl;
 328     cout<<" WARNING (MH): outputListHistos is NULL in AFAWMH::GOH() !!!!"<<endl;
 329     cout<<endl;
 330     exit(0);
 331    }
 332
 333 } // end of void AliFlowAnalysisWithMixedHarmonics::GetOutputHistograms(TList *outputListHistos)
 334
 335 //================================================================================================================
 336
 337 void AliFlowAnalysisWithMixedHarmonics::GetPointersForBaseHistograms()
 338 {
 339  // Get pointers to base histograms.
 340
 341  TString analysisSettingsName = "fAnalysisSettings";
 342  TProfile *analysisSettings = dynamic_cast<TProfile*>(fHistList->FindObject(analysisSettingsName.Data()));
 343  if(analysisSettings)
 344  {
 345   this->SetAnalysisSettings(analysisSettings);
 346  } else
 347    {
 348     cout<<endl;
 349     cout<<" WARNING (MH): analysisSettings is NULL in AFAWMH::GPFBH() !!!!"<<endl;
 350     cout<<endl;
 351     exit(0);
 352    }
 353
 354 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForBaseHistograms()
 355
 356 //================================================================================================================
 357
 358 void AliFlowAnalysisWithMixedHarmonics::GetPointersForCommonHistograms()
 359 {
 360  // Get pointers to common control histograms.
 361
 362  TString commonHistsName = "AliFlowCommonHistMH";
 363  AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(fHistList->FindObject(commonHistsName.Data()));
 364  if(commonHist)
 365  {
 366   this->SetCommonHists(commonHist);
 367  } else
 368    {
 369     cout<<endl;
 370     cout<<" WARNING (MH): commonHist is NULL in AFAWMH::GPFCH() !!!!"<<endl;
 371     cout<<endl;
 372     exit(0);
 373    }
 374
 375 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForCommonHistograms()
 376
 377 //================================================================================================================
 378
 379 void AliFlowAnalysisWithMixedHarmonics::GetPointersForAllEventProfiles()
 380 {
 381  // Get pointers to histograms holding final results.
 382
 383  TList *profileList = NULL;
 384  profileList = dynamic_cast<TList*>(fHistList->FindObject("Profiles"));
 385  if(!profileList)
 386  {
 387   cout<<"WARNING: profileList is NULL in AFAWMH::GPFAEP() !!!!"<<endl;
 388   exit(0);
 389  }
 390
 391  TString s3pCorrelatorProName = "f3pCorrelatorPro";
 392  TProfile *h3pCorrelatorPro = dynamic_cast<TProfile*>(profileList->FindObject(s3pCorrelatorProName.Data()));
 393  if(h3pCorrelatorPro)
 394  {
 395   this->Set3pCorrelatorPro(h3pCorrelatorPro);
 396  }
 397  TString s3pCorrelatorVsMProName = "f3pCorrelatorVsMPro";
 398  TProfile *h3pCorrelatorVsMPro = dynamic_cast<TProfile*>(profileList->FindObject(s3pCorrelatorVsMProName.Data()));
 399  if(h3pCorrelatorVsMPro)
 400  {
 401   this->Set3pCorrelatorVsMPro(h3pCorrelatorVsMPro);
 402  }
 403  TString nonIsotropicTermsProName = "fNonIsotropicTermsPro";
 404  TProfile *nonIsotropicTermsPro = dynamic_cast<TProfile*>(profileList->FindObject(nonIsotropicTermsProName.Data()));
 405  if(nonIsotropicTermsPro)
 406  {
 407   this->SetNonIsotropicTermsPro(nonIsotropicTermsPro);
 408  }
 409  TString nonIsotropicTermsVsMProName = "fNonIsotropicTermsVsMPro";
 410  TProfile2D *nonIsotropicTermsVsMPro = dynamic_cast<TProfile2D*>(profileList->FindObject(nonIsotropicTermsVsMProName.Data()));
 411  if(nonIsotropicTermsVsMPro)
 412  {
 413   this->SetNonIsotropicTermsVsMPro(nonIsotropicTermsVsMPro);
 414  }
 415 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForAllEventProfiles()
 416
 417 //================================================================================================================
 418
 419 void AliFlowAnalysisWithMixedHarmonics::GetPointersForResultsHistograms()
 420 {
 421  // Get pointers to histograms holding final results.
 422
 423  TList *resultsList = NULL;
 424  resultsList = dynamic_cast<TList*>(fHistList->FindObject("Results"));
 425  if(!resultsList)
 426  {
 427   cout<<"WARNING: resultsList is NULL in AFAWMH::GPFRH() !!!!"<<endl;
 428   exit(0);
 429  }
 430
 431  TString s3pCorrelatorHistName = "f3pCorrelatorHist";
 432  TH1D *h3pCorrelatorHist = dynamic_cast<TH1D*>(resultsList->FindObject(s3pCorrelatorHistName.Data()));
 433  if(h3pCorrelatorHist)
 434  {
 435   this->Set3pCorrelatorHist(h3pCorrelatorHist);
 436  }
 437  TString s3pCorrelatorVsMHistName = "f3pCorrelatorVsMHist";
 438  TH1D *h3pCorrelatorVsMHist = dynamic_cast<TH1D*>(resultsList->FindObject(s3pCorrelatorVsMHistName.Data()));
 439  if(h3pCorrelatorVsMHist)
 440  {
 441   this->Set3pCorrelatorVsMHist(h3pCorrelatorVsMHist);
 442  }
 443  TString detectorBiasHistName = "fDetectorBiasHist";
 444  TH1D *detectorBiasHist = dynamic_cast<TH1D*>(resultsList->FindObject(detectorBiasHistName.Data()));
 445  if(detectorBiasHist)
 446  {
 447   this->SetDetectorBiasHist(detectorBiasHist);
 448  }
 449  TString detectorBiasVsMHistName = "fDetectorBiasVsMHist";
 450  TH1D *detectorBiasVsMHist = dynamic_cast<TH1D*>(resultsList->FindObject(detectorBiasVsMHistName.Data()));
 451  if(detectorBiasVsMHist)
 452  {
 453   this->SetDetectorBiasVsMHist(detectorBiasVsMHist);
 454  }
 455  TString nonIsotropicTermsHistName = "fNonIsotropicTermsHist";
 456  TH1D *nonIsotropicTermsHist = dynamic_cast<TH1D*>(resultsList->FindObject(nonIsotropicTermsHistName.Data()));
 457  if(nonIsotropicTermsHist)
 458  {
 459   this->SetNonIsotropicTermsHist(nonIsotropicTermsHist);
 460  }
 461  TString nonIsotropicTermsVsMHistName = "fNonIsotropicTermsVsMHist";
 462  TH2D *nonIsotropicTermsVsMHist = dynamic_cast<TH2D*>(resultsList->FindObject(nonIsotropicTermsVsMHistName.Data()));
 463  if(nonIsotropicTermsVsMHist)
 464  {
 465   this->SetNonIsotropicTermsVsMHist(nonIsotropicTermsVsMHist);
 466  }
 467
 468 } // end of void AliFlowAnalysisWithMixedHarmonics::GetPointersForResultsHistograms()
 469
 470 //================================================================================================================
 471
 472 void AliFlowAnalysisWithMixedHarmonics::WriteHistograms(TString outputFileName)
 473 {
 474  // Store the final results in output .root file.
 475  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 476  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
 477  delete output;
 478 }
 479
 480 //================================================================================================================
 481
 482 void AliFlowAnalysisWithMixedHarmonics::WriteHistograms(TDirectoryFile *outputFileName)
 483 {
 484  // Store the final results in output .root file.
 485  fHistList->SetName("cobjMH");
 486  fHistList->SetOwner(kTRUE);
 487  outputFileName->Add(fHistList);
 488  outputFileName->Write(outputFileName->GetName(),TObject::kSingleKey);
 489 }
 490
 491 //================================================================================================================
 492
 493 void AliFlowAnalysisWithMixedHarmonics::BookAndNestAllLists()
 494 {
 495  // Book and nest all list in base list fHistList.
 496
 497  // Weights:
 498  fWeightsList->SetName("Weights");
 499  fWeightsList->SetOwner(kTRUE);
 500  fHistList->Add(fWeightsList);
 501  // Profiles:
 502  fProfileList->SetName("Profiles");
 503  fProfileList->SetOwner(kTRUE);
 504  fHistList->Add(fProfileList);
 505  // Results:
 506  fResultsList->SetName("Results");
 507  fResultsList->SetOwner(kTRUE);
 508  fHistList->Add(fResultsList);
 509
 510 } // end of void AliFlowAnalysisWithMixedHarmonics::BookAndNestAllLists()
 511
 512 //================================================================================================================
 513
 514 void AliFlowAnalysisWithMixedHarmonics::BookProfileHoldingSettings()
 515 {
 516  // Book profile to hold all analysis settings.
 517
 518  TString analysisSettingsName = "fAnalysisSettings";
 519  fAnalysisSettings = new TProfile(analysisSettingsName.Data(),"Settings for analysis with mixed harmonics",5,0,5);
 520  fAnalysisSettings->GetXaxis()->SetBinLabel(1,"Integer n in cos(n(2#phi_{1}-#phi_{2}-#phi_{3}))");
 521  fAnalysisSettings->Fill(0.5,fCorrelatorInteger);
 522  fAnalysisSettings->GetXaxis()->SetBinLabel(2,"Corrected for detector effects?");
 523  fAnalysisSettings->Fill(1.5,(Int_t)fCorrectForDetectorEffects);
 524  fAnalysisSettings->GetXaxis()->SetBinLabel(3,"# of multiplicity bins");
 525  fAnalysisSettings->Fill(2.5,fNoOfMultipicityBins);
 526  fAnalysisSettings->GetXaxis()->SetBinLabel(4,"Width of multiplicity bins");
 527  fAnalysisSettings->Fill(3.5,fMultipicityBinWidth);
 528  fAnalysisSettings->GetXaxis()->SetBinLabel(5,"Minimal multiplicity");
 529  fAnalysisSettings->Fill(4.5,fMinMultiplicity);
 530  fHistList->Add(fAnalysisSettings);
 531
 532 } // end of void AliFlowAnalysisWithMixedHarmonics::BookProfileHoldingSettings()
 533
 534 //================================================================================================================
 535
 536 void AliFlowAnalysisWithMixedHarmonics::BookCommonHistograms()
 537 {
 538  // Book common control histograms and common histograms for final results.
 539
 540  TString commonHistsName = "AliFlowCommonHistMH";
 541  fCommonHists = new AliFlowCommonHist(commonHistsName.Data());
 542  fHistList->Add(fCommonHists);
 543
 544 } // end of void AliFlowAnalysisWithMixedHarmonics::BookCommonHistograms()
 545
 546 //================================================================================================================
 547
 548 void AliFlowAnalysisWithMixedHarmonics::BookAllEventByEventQuantities()
 549 {
 550  // Book all event-by-event quantitites.
 551
 552  // Q_{n,k} and S{p,k}:
 553  fReQnk = new TMatrixD(2,4); // to be improved(bound on k)
 554  fImQnk = new TMatrixD(2,4); // to be improved(bound on k)
 555  fSpk = new TMatrixD(4,4); // to be improved(bound on p and k)
 556  // 3-p correlator <cos[n(2phi1-phi2-phi3)]> for single event:
 557  TString s3pCorrelatorEBEName = "f3pCorrelatorEBE";
 558  f3pCorrelatorEBE = new TH1D(s3pCorrelatorEBEName.Data(),"<cos[n(2#phi_{1}-#phi_{2}-#phi_{3})]> for single event",1,0,1);
 559  // Correction terms for non-uniform acceptance to <cos[n(2phi1-phi2-phi3)]> for single event:
 560  TString nonIsotropicTermsEBEName = "fNonIsotropicTermsEBE";
 561  fNonIsotropicTermsEBE = new TH1D(nonIsotropicTermsEBEName.Data(),"Non-isotropic terms for single event",10,0,10);
 562
 563 } // end fo void AliFlowAnalysisWithMixedHarmonics::BookAllEventByEventQuantities()
 564 //================================================================================================================
 565
 566 void AliFlowAnalysisWithMixedHarmonics::BookAllAllEventQuantities()
 567 {
 568  // Book all all-event quantitites.
 569
 570  // a) Book quantites without multiplicity binning;
 571  // b) Book quantites with multiplicity binning.
 572
 573  // a) Book quantites without multiplicity binning:
 574  // 3-p correlator <<cos[n(2phi1-phi2-phi3)]>> for all events (with wrong error):
 575  TString s3pCorrelatorProName = "f3pCorrelatorPro";
 576  f3pCorrelatorPro = new TProfile(s3pCorrelatorProName.Data(),"<<cos[n(2#phi_{1}-#phi_{2}-#phi_{3})]>> for all events (errors are wrong here!)",1,0,1,"s");
 577  fProfileList->Add(f3pCorrelatorPro);
 578  // 3-p correlator <<cos[n(2phi1-phi2-phi3)]>> for all events (with correct error):
 579  TString s3pCorrelatorHistName = "f3pCorrelatorHist";
 580  f3pCorrelatorHist = new TH1D(s3pCorrelatorHistName.Data(),"<<cos[n(2#phi_{1}-#phi_{2}-#phi_{3})]>> for all events",1,0,1);
 581  fResultsList->Add(f3pCorrelatorHist);
 582  // Correction terms for non-uniform acceptance to <cos[n(2phi1-phi2-phi3)]> for all events (with wrong errors):
 583  TString nonIsotropicTermsProName = "fNonIsotropicTermsPro";
 584  fNonIsotropicTermsPro = new TProfile(nonIsotropicTermsProName.Data(),"Non-isotropic terms for all events (errors are wrong here!)",10,0,10,"s");
 585  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(1,"cos(n(#phi_{1}))");
 586  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(2,"sin(n(#phi_{1}))");
 587  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(3,"cos(2n(#phi_{1}))");
 588  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(4,"sin(2n(#phi_{1}))");
 589  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(5,"cos(n(#phi_{1}+#phi_{2}))");
 590  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(6,"sin(n(#phi_{1}+#phi_{2}))");
 591  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(7,"cos(n(2#phi_{1}-#phi_{2}))");
 592  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(8,"sin(n(2#phi_{1}-#phi_{2}))");
 593  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(9,"cos(n(#phi_{1}-#phi_{2}-#phi_{3}))");
 594  fNonIsotropicTermsPro->GetXaxis()->SetBinLabel(10,"sin(n(#phi_{1}-#phi_{2}-#phi_{3})");
 595  fProfileList->Add(fNonIsotropicTermsPro);
 596  // Correction terms for non-uniform acceptance to <cos[n(2phi1-phi2-phi3)]> for all events (with correct errors):
 597  TString nonIsotropicTermsHistName = "fNonIsotropicTermsHist";
 598  fNonIsotropicTermsHist = new TH1D(nonIsotropicTermsHistName.Data(),"Non-isotropic terms for all events",10,0,10);
 599  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(1,"cos(n(#phi_{1}))");
 600  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(2,"sin(n(#phi_{1}))");
 601  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(3,"cos(2n(#phi_{1}))");
 602  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(4,"sin(2n(#phi_{1}))");
 603  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(5,"cos(n(#phi_{1}+#phi_{2}))");
 604  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(6,"sin(n(#phi_{1}+#phi_{2}))");
 605  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(7,"cos(n(2#phi_{1}-#phi_{2}))");
 606  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(8,"sin(n(2#phi_{1}-#phi_{2}))");
 607  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(9,"cos(n(#phi_{1}-#phi_{2}-#phi_{3}))");
 608  fNonIsotropicTermsHist->GetXaxis()->SetBinLabel(10,"sin(n(#phi_{1}-#phi_{2}-#phi_{3})");
 609  fResultsList->Add(fNonIsotropicTermsHist);
 610  // Quantified bias comming from detector inefficiencies to 3-p correlator <<cos[n(2phi1-phi2-phi3)]>> (in %)
 611  TString detectorBiasHistName = "fDetectorBiasHist";
 612  fDetectorBiasHist = new TH1D(detectorBiasHistName.Data(),"Bias coming from detector inefficiences",1,0,1);
 613  fDetectorBiasHist->GetXaxis()->SetBinLabel(1,"#left|1 - #frac{corrected}{not-corrected}#right| (in %)");
 614  fResultsList->Add(fDetectorBiasHist);
 615
 616  // b) Book quantites with multiplicity binning.
 617  // 3-p correlator <<cos[n(2phi1-phi2-phi3)]>> for all events versus multiplicity (with wrong error):
 618  TString s3pCorrelatorVsMProName = "f3pCorrelatorVsMPro";
 619  f3pCorrelatorVsMPro = new TProfile(s3pCorrelatorVsMProName.Data(),"<<cos[n(2#phi_{1}-#phi_{2}-#phi_{3})]>> vs M for all events (errors are wrong here!)",fNoOfMultipicityBins+1,0,fNoOfMultipicityBins+1,"s");
 620  for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
 621  {
 622   f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+b*fMultipicityBinWidth)));
 623  }
 624  f3pCorrelatorVsMPro->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+1,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 625  fProfileList->Add(f3pCorrelatorVsMPro);
 626  // 3-p correlator <<cos[n(2phi1-phi2-phi3)]>> for all events (with correct error):
 627  TString s3pCorrelatorVsMHistName = "f3pCorrelatorVsMHist";
 628  f3pCorrelatorVsMHist = new TH1D(s3pCorrelatorVsMHistName.Data(),"<<cos[n(2#phi_{1}-#phi_{2}-#phi_{3})]>> vs M for all events",fNoOfMultipicityBins+1,0,fNoOfMultipicityBins+1);
 629  for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
 630  {
 631   f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+b*fMultipicityBinWidth)));
 632  }
 633  f3pCorrelatorVsMHist->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+1,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 634  fResultsList->Add(f3pCorrelatorVsMHist);
 635  // Correction terms for non-uniform acceptance to <cos[n(2phi1-phi2-phi3)]> for all events vs M (with wrong errors):
 636  TString nonIsotropicTermsVsMProName = "fNonIsotropicTermsVsMPro";
 637  fNonIsotropicTermsVsMPro = new TProfile2D(nonIsotropicTermsVsMProName.Data(),"Non-isotropic terms for all events vs M (errors are wrong here!)",10,0,10,fNoOfMultipicityBins+1,0,fNoOfMultipicityBins+1,"s");
 638  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(1,"cos(n(#phi_{1}))");
 639  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(2,"sin(n(#phi_{1}))");
 640  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(3,"cos(2n(#phi_{1}))");
 641  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(4,"sin(2n(#phi_{1}))");
 642  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(5,"cos(n(#phi_{1}+#phi_{2}))");
 643  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(6,"sin(n(#phi_{1}+#phi_{2}))");
 644  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(7,"cos(n(2#phi_{1}-#phi_{2}))");
 645  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(8,"sin(n(2#phi_{1}-#phi_{2}))");
 646  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(9,"cos(n(#phi_{1}-#phi_{2}-#phi_{3}))");
 647  fNonIsotropicTermsVsMPro->GetXaxis()->SetBinLabel(10,"sin(n(#phi_{1}-#phi_{2}-#phi_{3})");
 648  for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
 649  {
 650   fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+b*fMultipicityBinWidth)));
 651  }
 652  fNonIsotropicTermsVsMPro->GetYaxis()->SetBinLabel(fNoOfMultipicityBins+1,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 653  fProfileList->Add(fNonIsotropicTermsVsMPro);
 654  // Correction terms for non-uniform acceptance to <cos[n(2phi1-phi2-phi3)]> for all events (with correct errors):
 655  TString nonIsotropicTermsVsMHistName = "fNonIsotropicTermsVsMHist";
 656  fNonIsotropicTermsVsMHist = new TH2D(nonIsotropicTermsVsMHistName.Data(),"Non-isotropic terms for all events vs M ",10,0,10,fNoOfMultipicityBins+1,0,fNoOfMultipicityBins+1);
 657  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(1,"cos(n(#phi_{1}))");
 658  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(2,"sin(n(#phi_{1}))");
 659  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(3,"cos(2n(#phi_{1}))");
 660  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(4,"sin(2n(#phi_{1}))");
 661  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(5,"cos(n(#phi_{1}+#phi_{2}))");
 662  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(6,"sin(n(#phi_{1}+#phi_{2}))");
 663  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(7,"cos(n(2#phi_{1}-#phi_{2}))");
 664  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(8,"sin(n(2#phi_{1}-#phi_{2}))");
 665  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(9,"cos(n(#phi_{1}-#phi_{2}-#phi_{3}))");
 666  fNonIsotropicTermsVsMHist->GetXaxis()->SetBinLabel(10,"sin(n(#phi_{1}-#phi_{2}-#phi_{3})");
 667  for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
 668  {
 669   fNonIsotropicTermsVsMHist->GetYaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+b*fMultipicityBinWidth)));
 670  }
 671  fNonIsotropicTermsVsMHist->GetYaxis()->SetBinLabel(fNoOfMultipicityBins+1,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 672  fResultsList->Add(fNonIsotropicTermsVsMHist);
 673  // Quantified bias comming from detector inefficiencies to 3-p correlator <<cos[n(2phi1-phi2-phi3)]>> (in %)
 674  TString detectorBiasVsMHistName = "fDetectorBiasVsMHist";
 675  fDetectorBiasVsMHist = new TH1D(detectorBiasVsMHistName.Data(),"#left|1 - #frac{corrected}{not-corrected}#right| (in %)",fNoOfMultipicityBins+1,0,fNoOfMultipicityBins+1);
 676  for(Int_t b=1;b<=fNoOfMultipicityBins;b++)
 677  {
 678   fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(b,Form("%d #leq M < %d",(Int_t)(fMinMultiplicity+(b-1)*fMultipicityBinWidth),(Int_t)(fMinMultiplicity+b*fMultipicityBinWidth)));
 679  }
 680  fDetectorBiasVsMHist->GetXaxis()->SetBinLabel(fNoOfMultipicityBins+1,Form(" M #geq %d",(Int_t)(fMinMultiplicity+fNoOfMultipicityBins*fMultipicityBinWidth)));
 681  fResultsList->Add(fDetectorBiasVsMHist);
 682
 683 } // end of void AliFlowAnalysisWithMixedHarmonics::BookAllAllEventQuantities()
 684
 685 //================================================================================================================
 686
 687 void AliFlowAnalysisWithMixedHarmonics::AccessConstants()
 688 {
 689  // Access needed common constants from AliFlowCommonConstants.
 690
 691  fnBinsPhi = AliFlowCommonConstants::GetMaster()->GetNbinsPhi();
 692  fPhiMin = AliFlowCommonConstants::GetMaster()->GetPhiMin();
 693  fPhiMax = AliFlowCommonConstants::GetMaster()->GetPhiMax();
 694  if(fnBinsPhi) fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;
 695  fnBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 696  fPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 697  fPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 698  if(fnBinsPt) fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;
 699  fnBinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 700  fEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 701  fEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 702  if(fnBinsEta) fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;
 703
 704 } // end of void AliFlowAnalysisWithMixedHarmonics::AccessConstants()
 705
 706 //================================================================================================================
 707
 708 void AliFlowAnalysisWithMixedHarmonics::CrossCheckSettings()
 709 {
 710  // Cross-check if the user settings make sense.
 711
 712 } // end of void AliFlowAnalysisWithMixedHarmonics::CrossCheckSettings()
 713
 714 //================================================================================================================
 715
 716 void AliFlowAnalysisWithMixedHarmonics::BookAndFillWeightsHistograms()
 717 {
 718  // Book and fill (by accessing file "weights.root") histograms which hold phi, pt and eta weights.
 719
 720  if(!fWeightsList)
 721  {
 722   cout<<"WARNING: fWeightsList is NULL in AFAWMH::BAFWH() !!!!"<<endl;
 723   exit(0);
 724  }
 725  // Profile to hold flags for weights:
 726  TString fUseParticleWeightsName = "fUseParticleWeightsMH";
 727  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
 728  fUseParticleWeights->SetLabelSize(0.06);
 729  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
 730  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
 731  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
 732  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
 733  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
 734  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
 735  fWeightsList->Add(fUseParticleWeights);
 736  // Phi-weights:
 737  if(fUsePhiWeights)
 738  {
 739   if(fWeightsList->FindObject("phi_weights"))
 740   {
 741    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
 742    if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
 743    {
 744     cout<<endl;
 745     cout<<"WARNING (MH): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
 746     cout<<endl;
 747     exit(0);
 748    }
 749   } else
 750     {
 751      cout<<"WARNING (MH): fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWMH::BAFWH() !!!!"<<endl;
 752      exit(0);
 753     }
 754  } // end of if(fUsePhiWeights)
 755  // Pt-weights:
 756  if(fUsePtWeights)
 757  {
 758   if(fWeightsList->FindObject("pt_weights"))
 759   {
 760    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
 761    if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
 762    {
 763     cout<<endl;
 764     cout<<"WARNING (MH): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
 765     cout<<endl;
 766     exit(0);
 767    }
 768   } else
 769     {
 770      cout<<"WARNING (MH): fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWMH::BAFWH() !!!!"<<endl;
 771      exit(0);
 772     }
 773  } // end of if(fUsePtWeights)
 774  // Eta-weights:
 775  if(fUseEtaWeights)
 776  {
 777   if(fWeightsList->FindObject("eta_weights"))
 778   {
 779    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
 780    if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
 781    {
 782     cout<<endl;
 783     cout<<"WARNING (MH): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
 784     cout<<endl;
 785     exit(0);
 786    }
 787   } else
 788     {
 789      cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWMH::BAFWH() !!!!"<<endl;
 790      exit(0);
 791     }
 792  } // end of if(fUseEtaWeights)
 793
 794 } // end of AliFlowAnalysisWithMixedHarmonics::BookAndFillWeightsHistograms()
 795
 796 //================================================================================================================
 797
 798 void AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInMake()
 799 {
 800  // Check pointers used in method Make().
 801
 802  if(!fReQnk || !fImQnk || !fSpk || !f3pCorrelatorEBE)
 803  {
 804   cout<<endl;
 805   cout<<" WARNING (MH): (!fReQnk || !fImQnk || !fSpk || !f3pCorrelatorEBE) is NULL !!!!"<<endl;
 806   cout<<endl;
 807   exit(0);
 808  }
 809  if(!f3pCorrelatorPro)
 810  {
 811   cout<<endl;
 812   cout<<" WARNING (MH): (!f3pCorrelatorPro) is NULL !!!!"<<endl;
 813   cout<<endl;
 814   exit(0);
 815  }
 816
 817 } // end of AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInMake()
 818
 819 //================================================================================================================
 820
 821 void AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInFinish()
 822 {
 823  // Check pointers used in method Finish().
 824
 825  if(!fAnalysisSettings)
 826  {
 827   cout<<endl;
 828   cout<<" WARNING (MH): !fAnalysisSettings is NULL !!!!"<<endl;
 829   cout<<endl;
 830   exit(0);
 831  }
 832  if(!f3pCorrelatorPro)
 833  {
 834   cout<<endl;
 835   cout<<" WARNING (MH): !f3pCorrelatorPro is NULL !!!!"<<endl;
 836   cout<<endl;
 837   exit(0);
 838  }
 839  if(!fNonIsotropicTermsPro)
 840  {
 841   cout<<endl;
 842   cout<<" WARNING (MH): !fNonIsotropicTermsPro is NULL !!!!"<<endl;
 843   cout<<endl;
 844   exit(0);
 845  }
 846  if(!f3pCorrelatorVsMPro)
 847  {
 848   cout<<endl;
 849   cout<<" WARNING (MH): !f3pCorrelatorVsMPro is NULL !!!!"<<endl;
 850   cout<<endl;
 851   exit(0);
 852  }
 853  if(!fNonIsotropicTermsVsMPro)
 854  {
 855   cout<<endl;
 856   cout<<" WARNING (MH): !fNonIsotropicTermsVsMPro is NULL !!!!"<<endl;
 857   cout<<endl;
 858   exit(0);
 859  }
 860  if(!f3pCorrelatorHist)
 861  {
 862   cout<<endl;
 863   cout<<" WARNING (MH): !f3pCorrelatorHist is NULL !!!!"<<endl;
 864   cout<<endl;
 865   exit(0);
 866  }
 867  if(!fNonIsotropicTermsHist)
 868  {
 869   cout<<endl;
 870   cout<<" WARNING (MH): !fNonIsotropicTermsHist is NULL !!!!"<<endl;
 871   cout<<endl;
 872   exit(0);
 873  }
 874  if(!fDetectorBiasHist)
 875  {
 876   cout<<endl;
 877   cout<<" WARNING (MH): !fDetectorBiasHist is NULL !!!!"<<endl;
 878   cout<<endl;
 879   exit(0);
 880  }
 881  if(!f3pCorrelatorVsMHist)
 882  {
 883   cout<<endl;
 884   cout<<" WARNING (MH): !f3pCorrelatorVsMHist is NULL !!!!"<<endl;
 885   cout<<endl;
 886   exit(0);
 887  }
 888  if(!fNonIsotropicTermsVsMHist)
 889  {
 890   cout<<endl;
 891   cout<<" WARNING (MH): !fNonIsotropicTermsVsMHist is NULL !!!!"<<endl;
 892   cout<<endl;
 893   exit(0);
 894  }
 895  if(!fDetectorBiasVsMHist)
 896  {
 897   cout<<endl;
 898   cout<<" WARNING (MH): !fDetectorBiasVsMHist is NULL !!!!"<<endl;
 899   cout<<endl;
 900   exit(0);
 901  }
 902
 903 } // end of AliFlowAnalysisWithMixedHarmonics::CheckPointersUsedInFinish()
 904
 905 //================================================================================================================
 906
 907 void AliFlowAnalysisWithMixedHarmonics::AccessSettings()
 908 {
 909  // Access the settings for analysis with mixed harmonics.
 910
 911  fCorrelatorInteger = (Int_t)fAnalysisSettings->GetBinContent(1);
 912  fCorrectForDetectorEffects = (Bool_t)fAnalysisSettings->GetBinContent(2);
 913  fNoOfMultipicityBins = (Int_t)fAnalysisSettings->GetBinContent(3);
 914  fMultipicityBinWidth = (Double_t)fAnalysisSettings->GetBinContent(4);
 915  fMinMultiplicity = (Double_t)fAnalysisSettings->GetBinContent(5);
 916
 917 } // end of AliFlowAnalysisWithMixedHarmonics::AccessSettings()
 918
 919 //================================================================================================================
 920
 921 void AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffects()
 922 {
 923  // Correct measured 3-p correlator cos[n(2phi1-phi2-phi3)] for detector effects.
 924
 925  Double_t measured3pCorrelator = f3pCorrelatorPro->GetBinContent(1); // biased by detector effects
 926  Double_t corrected3pCorrelator = 0.; // corrected for detector effects
 927  Double_t nonIsotropicTerms[10] = {0.}; // there are 10 distinct non-isotropic terms
 928  for(Int_t nit=0;nit<10;nit++)
 929  {
 930   nonIsotropicTerms[nit] = fNonIsotropicTermsPro->GetBinContent(nit+1);
 931  }
 932  // Calculate corrected 3-p correlator:
 933  corrected3pCorrelator = measured3pCorrelator
 934                        - nonIsotropicTerms[2]*nonIsotropicTerms[4]
 935                        - nonIsotropicTerms[3]*nonIsotropicTerms[5]
 936                        - 2.*nonIsotropicTerms[0]*nonIsotropicTerms[6]
 937                        - 2.*nonIsotropicTerms[1]*nonIsotropicTerms[7]
 938                        + 2.*nonIsotropicTerms[2]*(pow(nonIsotropicTerms[0],2.)-pow(nonIsotropicTerms[1],2.))
 939                        + 4.*nonIsotropicTerms[3]*nonIsotropicTerms[0]*nonIsotropicTerms[1];
 940  // Store corrected 3-p correlator:
 941  f3pCorrelatorHist->SetBinContent(1,corrected3pCorrelator);
 942
 943 } // end of AliFlowAnalysisWithMixedHarmonics::CorrectForDetectorEffects()
 944
 945 //================================================================================================================
 946
 947 void AliFlowAnalysisWithMixedHarmonics::FinalizeNonIsotropicTerms()
 948 {
 949  // Calculate correctly the errors of non-isotropic terms.
 950
 951  // Transfer mean values from profiles to histograms: // to be improved
 952  for(Int_t mv=1;mv<=10;mv++)
 953  {
 954   fNonIsotropicTermsHist->SetBinContent(mv,fNonIsotropicTermsPro->GetBinContent(mv));
 955  }
 956
 957 } // end of void AliFlowAnalysisWithMixedHarmonics::FinalizeNonIsotropicTerms()
 958
 959 //================================================================================================================
 960
 961 void AliFlowAnalysisWithMixedHarmonics::QuantifyBiasFromDetectorEffects()
 962 {
 963  // Quantify bias from detector inefficiences to 3-p correlator.
 964
 965  // Remark: Bias is quantified as 1 - |corrected/nonCorrected| (in %) and stored in histogram fDetectorBias
 966
 967  Double_t bias = 0.;
 968  if(f3pCorrelatorPro->GetBinContent(1))
 969  {
 970   bias = 100.*TMath::Abs(1.-f3pCorrelatorHist->GetBinContent(1)/f3pCorrelatorPro->GetBinContent(1));
 971  }
 972  fDetectorBiasHist->SetBinContent(1,bias);
 973
 974 } // end of void AliFlowAnalysisWithMixedHarmonics::QuantifyBiasFromDetectorEffects()
 975
 976 //================================================================================================================
 977 void AliFlowAnalysisWithMixedHarmonics::ResetEventByEventQuantities()
 978 {
 979  // Reset all event by event quantities.
 980
 981  fReQnk->Zero();
 982  fImQnk->Zero();
 983  fSpk->Zero();
 984
 985  f3pCorrelatorEBE->Reset();
 986
 987 } // end of void AliFlowAnalysisWithMixedHarmonics::ResetEventByEventQuantities()
 988
 989 //================================================================================================================
 990
 991 void AliFlowAnalysisWithMixedHarmonics::Calculate3pCorrelator()
 992 {
 993  // Calculate 3-p azimuthal correlator cos[n(2phi1-phi2-phi3)] in terms of Q_{n,k} and S_{p,k}.
 994
 995  // a) Calculate 3-p correlator without using particle weights;
 996  // b) Calculate 3-p correlator with using particle weights.
 997
 998  // a) Calculate 3-p correlator without using particle weights:
 999  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1000  {
1001   // Multiplicity (number of RPs):
1002   Double_t dMult = (*fSpk)(0,0);
1003   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n and 2n:
1004   Double_t dReQ1n = (*fReQnk)(0,0);
1005   Double_t dReQ2n = (*fReQnk)(1,0);
1006   Double_t dImQ1n = (*fImQnk)(0,0);
1007   Double_t dImQ2n = (*fImQnk)(1,0);
1008   // 3-particle azimuthal correlator <cos(n*(2.*phi1-phi2-phi3))>:
1009   Double_t three2n1n1n = (pow(dReQ1n,2.)*dReQ2n + 2.*dReQ1n*dImQ1n*dImQ2n - pow(dImQ1n,2.)*dReQ2n
1010                        - 2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
1011                        - (pow(dReQ2n,2.)+pow(dImQ2n,2.))+2.*dMult)
1012                        / (dMult*(dMult-1.)*(dMult-2.));
1013   // Fill event-by-event histogram and all-events profile:
1014   f3pCorrelatorEBE->SetBinContent(1,three2n1n1n);
1015   f3pCorrelatorPro->Fill(0.5,three2n1n1n,dMult*(dMult-1.)*(dMult-2.));
1016
1017   // Correlator vs multiplicity: // to be improved
1018   f3pCorrelatorVsMPro->Fill(0.5+(Int_t)((dMult-fMinMultiplicity)/fMultipicityBinWidth),dMult*(dMult-1.)*(dMult-2.)); // to be improved (cross-checked)
1019  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1020
1021  // b) Calculate 3-p correlator without using particle weights:
1022  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1023  {
1024
1025  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1026
1027 } // end of void AliFlowAnalysisWithMixedHarmonics::Calculate3pCorrelator()
1028
1029 //================================================================================================================
1030
1031 void AliFlowAnalysisWithMixedHarmonics::CalculateNonIsotropicTerms()
1032 {
1033  // Calculate correction terms for non-uniform acceptance to 3-p correlator <cos[n(2phi1-phi2-phi3)]>.
1034
1035  // a) Calculate using particle weights;
1036  // b) Calculate without using particle weights.
1037
1038  // These non-isotropic terms are stored in fNonIsotropicTermsPro. The same results but with
1039  // correct errors are stored in fNonIsotropicTermsHist. Binning of fNonIsotropicTermsPro and
1040  // fNonIsotropicTermsHist is organized as follows:
1041  //  1st bin: <<cos(n*phi1)>>
1042  //  2nd bin: <<sin(n*phi1)>>
1043  //  3rd bin: <<cos(2n*phi1)>>
1044  //  4th bin: <<sin(2n*phi1)>>
1045  //  5th bin: <<cos(n*(phi1+phi2)>>
1046  //  6th bin: <<sin(n*(phi1+phi2)>>
1047  //  7th bin: <<cos(n*(2phi1-phi2)>>
1048  //  8th bin: <<sin(n*(2phi1-phi2)>>
1049  //  9th bin: <<cos(n*(phi1-phi2-phi3)>>
1050  // 10th bin: <<sin(n*(phi1-phi2-phi3)>>
1051
1052  // a) Calculate using particle weights:
1053  if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1054  {
1055   // Multiplicity (number of RPs):
1056   Double_t dMult = (*fSpk)(0,0);
1057   // Real and imaginary parts of non-weighted Q-vectors (Q_{n,0}) evaluated in harmonics n and 2n:
1058   Double_t dReQ1n = (*fReQnk)(0,0);
1059   Double_t dReQ2n = (*fReQnk)(1,0);
1060   Double_t dImQ1n = (*fImQnk)(0,0);
1061   Double_t dImQ2n = (*fImQnk)(1,0);
1062   // 1-particle terms:
1063   Double_t cosP1n = 0.; // <cos(n*(phi1))>
1064   Double_t sinP1n = 0.; // <sin(n*(phi1))>
1065   Double_t cosP2n = 0.; // <cos(2n*(phi1))>
1066   Double_t sinP2n = 0.; // <sin(2n*(phi1))>
1067   if(dMult>0)
1068   {
1069    cosP1n = dReQ1n/dMult;
1070    sinP1n = dImQ1n/dMult;
1071    cosP2n = dReQ2n/dMult;
1072    sinP2n = dImQ2n/dMult;
1073    // Single event avarages:
1074    fNonIsotropicTermsEBE->SetBinContent(1,cosP1n); // <cos(n*(phi1))>
1075    fNonIsotropicTermsEBE->SetBinContent(2,sinP1n); // <sin(n*(phi1))>
1076    fNonIsotropicTermsEBE->SetBinContent(3,cosP2n); // <cos(2n*(phi1))>
1077    fNonIsotropicTermsEBE->SetBinContent(4,sinP2n); // <sin(2n*(phi1))>
1078    // All-events avarages:
1079    fNonIsotropicTermsPro->Fill(0.5,cosP1n,dMult); // <<cos(n*(phi1))>>
1080    fNonIsotropicTermsPro->Fill(1.5,sinP1n,dMult); // <<sin(n*(phi1))>>
1081    fNonIsotropicTermsPro->Fill(2.5,cosP2n,dMult); // <<cos(2n*(phi1))>>
1082    fNonIsotropicTermsPro->Fill(3.5,sinP2n,dMult); // <<sin(2n*(phi1))>>
1083   }
1084   // 2-particle terms:
1085   Double_t cosP1nP1n = 0.; // <cos(n*(phi1+phi2))>
1086   Double_t sinP1nP1n = 0.; // <sin(n*(phi1+phi2))>
1087   Double_t cosP2nM1n = 0.; // <cos(n*(2phi1-phi2))>
1088   Double_t sinP2nM1n = 0.; // <sin(n*(2phi1-phi2))>
1089   if(dMult>1)
1090   {
1091    cosP1nP1n = (pow(dReQ1n,2)-pow(dImQ1n,2)-dReQ2n)/(dMult*(dMult-1));
1092    sinP1nP1n = (2.*dReQ1n*dImQ1n-dImQ2n)/(dMult*(dMult-1));
1093    cosP2nM1n = (dReQ2n*dReQ1n+dImQ2n*dImQ1n-dReQ1n)/(dMult*(dMult-1));
1094    sinP2nM1n = (dImQ2n*dReQ1n-dReQ2n*dImQ1n-dImQ1n)/(dMult*(dMult-1));
1095    // Single event avarages:
1096    fNonIsotropicTermsEBE->SetBinContent(5,cosP1nP1n); // <cos(n*(phi1+phi2))>
1097    fNonIsotropicTermsEBE->SetBinContent(6,sinP1nP1n); // <sin(n*(phi1+phi2))>
1098    fNonIsotropicTermsEBE->SetBinContent(7,cosP2nM1n); // <cos(n*(2phi1-phi2))>
1099    fNonIsotropicTermsEBE->SetBinContent(8,sinP2nM1n); // <sin(n*(2phi1-phi2))>
1100    // All-events avarages:
1101    fNonIsotropicTermsPro->Fill(4.5,cosP1nP1n,dMult*(dMult-1.)); // <<cos(n*(phi1+phi2))>>
1102    fNonIsotropicTermsPro->Fill(5.5,sinP1nP1n,dMult*(dMult-1.)); // <<sin(n*(phi1+phi2))>>
1103    fNonIsotropicTermsPro->Fill(6.5,cosP2nM1n,dMult*(dMult-1.)); // <<cos(n*(2phi1-phi2))>>
1104    fNonIsotropicTermsPro->Fill(7.5,sinP2nM1n,dMult*(dMult-1.)); // <<sin(n*(2phi1-phi2))>>
1105   }
1106
1107   // 3-particle:
1108   Double_t cosP1nM1nM1n = 0.; // <cos(n*(phi1-phi2-phi3))>
1109   Double_t sinP1nM1nM1n = 0.; // <sin(n*(phi1-phi2-phi3))>
1110   if(dMult>2)
1111   {
1112    cosP1nM1nM1n = (dReQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))-dReQ1n*dReQ2n-dImQ1n*dImQ2n-2.*(dMult-1)*dReQ1n)
1113                 / (dMult*(dMult-1)*(dMult-2));
1114    sinP1nM1nM1n = (-dImQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))+dReQ1n*dImQ2n-dImQ1n*dReQ2n+2.*(dMult-1)*dImQ1n)
1115                 / (dMult*(dMult-1)*(dMult-2));
1116    // Single event avarages:
1117    fNonIsotropicTermsEBE->SetBinContent(9,cosP1nM1nM1n); // <cos(n*(phi1-phi2-phi3))>
1118    fNonIsotropicTermsEBE->SetBinContent(10,sinP1nM1nM1n); // <sin(n*(phi1-phi2-phi3))>
1119    // All-events avarages:
1120    fNonIsotropicTermsPro->Fill(8.5,cosP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<cos(n*(phi1-phi2-phi3))>>
1121    fNonIsotropicTermsPro->Fill(9.5,sinP1nM1nM1n,dMult*(dMult-1.)*(dMult-2.)); // <<sin(n*(phi1-phi2-phi3))>>
1122   }
1123  } // end of if(!(fUsePhiWeights || fUsePtWeights || fUseEtaWeights))
1124
1125  // b) Calculate without using particle weights:
1126  if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1127  {
1128
1129  } // end of if(fUsePhiWeights || fUsePtWeights || fUseEtaWeights)
1130
1131 } // end of void AliFlowAnalysisWithMixedHarmonics::CalculateNonIsotropicTerms()
1132