PWG/FLOW/Base/AliFlowAnalysisWithMCEventPlane.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 // AliFlowAnalysisWithMCEventPlane:
  16 // Description: Maker to analyze Flow from the generated MC reaction plane.
  17 //              This class is used to get the real value of the flow
  18 //              to compare the other methods to when analysing simulated events
  19 // author: N. van der Kolk (kolk@nikhef.nl)
  20
  21 #define AliFlowAnalysisWithMCEventPlane_cxx
  22
  23 #include "Riostream.h"
  24 #include "TFile.h"
  25 #include "TProfile.h"
  26 #include "TProfile2D.h"
  27 #include "TList.h"
  28 #include "TH1F.h"
  29 #include "TMath.h"
  30 #include "TVector2.h"
  31
  32 #include "AliFlowCommonConstants.h"
  33 #include "AliFlowEventSimple.h"
  34 #include "AliFlowTrackSimple.h"
  35 #include "AliFlowCommonHist.h"
  36 #include "AliFlowCommonHistResults.h"
  37 #include "AliFlowAnalysisWithMCEventPlane.h"
  38 #include "AliFlowVector.h"
  39
  40 class AliFlowVector;
  41
  42 using std::endl;
  43 using std::cout;
  44 ClassImp(AliFlowAnalysisWithMCEventPlane)
  45
  46   //-----------------------------------------------------------------------
  47
  48  AliFlowAnalysisWithMCEventPlane::AliFlowAnalysisWithMCEventPlane():
  49    fQsum(NULL),
  50    fQ2sum(0),
  51    fEventNumber(0),
  52    fDebug(kFALSE),
  53    fHistList(NULL),
  54    fCommonHists(NULL),
  55    fCommonHistsRes(NULL),
  56    fHistRP(NULL),
  57    fHistProIntFlow(NULL),
  58    fHistProIntFlowVsM(NULL),
  59    fHistProDiffFlowPtEtaRP(NULL),
  60    fHistProDiffFlowPtRP(NULL),
  61    fHistProDiffFlowEtaRP(NULL),
  62    fHistProDiffFlowPtEtaPOI(NULL),
  63    fHistProDiffFlowPtPOI(NULL),
  64    fHistProDiffFlowEtaPOI(NULL),
  65    fHistSpreadOfFlow(NULL),
  66    fHarmonic(2),
  67    fMixedHarmonicsList(NULL),
  68    fEvaluateMixedHarmonics(kFALSE),
  69    fMixedHarmonicsSettings(NULL),
  70    fnBinsMult(10000),
  71    fMinMult(0.),
  72    fMaxMult(10000.),
  73    fNinCorrelator(2),
  74    fMinCorrelator(2),
  75    fXinPairAngle(0.5)
  76 {
  77
  78   // Constructor.
  79   fHistList = new TList();
  80
  81   fQsum = new TVector2;        // flow vector sum
  82
  83   fMixedHarmonicsList = new TList();
  84   this->InitalizeArraysForMixedHarmonics();
  85
  86 }
  87
  88 //-----------------------------------------------------------------------
  89
  90 AliFlowAnalysisWithMCEventPlane::~AliFlowAnalysisWithMCEventPlane()
  91 {
  92  //destructor
  93   delete fHistList;
  94   delete fQsum;
  95 }
  96
  97 //-----------------------------------------------------------------------
  98
  99 void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString* outputFileName)
 100 {
 101  //store the final results in output .root file
 102  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 103  //output->WriteObject(fHistList, "cobjMCEP","SingleKey");
 104  fHistList->SetName("cobjMCEP");
 105  fHistList->SetOwner(kTRUE);
 106  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 107  delete output;
 108 }
 109
 110 //-----------------------------------------------------------------------
 111
 112 void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString outputFileName)
 113 {
 114  //store the final results in output .root file
 115  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 116  //output->WriteObject(fHistList, "cobjMCEP","SingleKey");
 117  fHistList->SetName("cobjMCEP");
 118  fHistList->SetOwner(kTRUE);
 119  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 120  delete output;
 121 }
 122
 123 //-----------------------------------------------------------------------
 124
 125 void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TDirectoryFile *outputFileName)
 126 {
 127  //store the final results in output .root file
 128  fHistList->SetName("cobjMCEP");
 129  fHistList->SetOwner(kTRUE);
 130  outputFileName->Add(fHistList);
 131  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 132 }
 133
 134 //-----------------------------------------------------------------------
 135 void AliFlowAnalysisWithMCEventPlane::Init() {
 136
 137   //Define all histograms
 138   cout<<"---Analysis with the real MC Event Plane---"<<endl;
 139
 140   //save old value and prevent histograms from being added to directory
 141   //to avoid name clashes in case multiple analaysis objects are used
 142   //in an analysis
 143   Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 144   TH1::AddDirectory(kFALSE);
 145
 146   Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 147   Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 148   Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 149
 150   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 151   Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 152   Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 153
 154   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistMCEP");
 155   fHistList->Add(fCommonHists);
 156   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsMCEP","",fHarmonic);
 157   fHistList->Add(fCommonHistsRes);
 158
 159   // store harmonic in common control histogram:
 160   (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
 161
 162   fHistRP = new TH1F("Flow_RP_MCEP","Flow_RP_MCEP",100,0.,3.14);
 163   fHistRP->SetXTitle("Reaction Plane Angle");
 164   fHistRP->SetYTitle("Counts");
 165   fHistList->Add(fHistRP);
 166
 167   fHistProIntFlow = new TProfile("FlowPro_V_MCEP","FlowPro_V_MCEP",1,0.,1.);
 168   fHistProIntFlow->SetLabelSize(0.06);
 169   (fHistProIntFlow->GetXaxis())->SetBinLabel(1,"v_{n}{MCEP}");
 170   fHistProIntFlow->SetYTitle("");
 171   fHistList->Add(fHistProIntFlow);
 172
 173   fHistProIntFlowVsM = new TProfile("FlowPro_VsM_MCEP","FlowPro_VsM_MCEP",10000,0.,10000.); // to be improved - hardwired 10000
 174   //fHistProIntFlowVsM->SetLabelSize(0.06);
 175   (fHistProIntFlowVsM->GetXaxis())->SetTitle("M");
 176   fHistProIntFlowVsM->SetYTitle("");
 177   fHistList->Add(fHistProIntFlowVsM);
 178
 179   fHistProDiffFlowPtEtaRP = new TProfile2D("FlowPro_VPtEtaRP_MCEP","FlowPro_VPtEtaRP_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
 180   fHistProDiffFlowPtEtaRP->SetXTitle("P_{t}");
 181   fHistProDiffFlowPtEtaRP->SetYTitle("#eta");
 182   fHistList->Add(fHistProDiffFlowPtEtaRP);
 183
 184   fHistProDiffFlowPtRP = new TProfile("FlowPro_VPtRP_MCEP","FlowPro_VPtRP_MCEP",iNbinsPt,dPtMin,dPtMax);
 185   fHistProDiffFlowPtRP->SetXTitle("P_{t}");
 186   fHistProDiffFlowPtRP->SetYTitle("");
 187   fHistList->Add(fHistProDiffFlowPtRP);
 188
 189   fHistProDiffFlowEtaRP = new TProfile("FlowPro_VetaRP_MCEP","FlowPro_VetaRP_MCEP",iNbinsEta,dEtaMin,dEtaMax);
 190   fHistProDiffFlowEtaRP->SetXTitle("#eta");
 191   fHistProDiffFlowEtaRP->SetYTitle("");
 192   fHistList->Add(fHistProDiffFlowEtaRP);
 193
 194   fHistProDiffFlowPtEtaPOI = new TProfile2D("FlowPro_VPtEtaPOI_MCEP","FlowPro_VPtEtaPOI_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
 195   fHistProDiffFlowPtEtaPOI->SetXTitle("P_{t}");
 196   fHistProDiffFlowPtEtaPOI->SetYTitle("#eta");
 197   fHistList->Add(fHistProDiffFlowPtEtaPOI);
 198
 199   fHistProDiffFlowPtPOI = new TProfile("FlowPro_VPtPOI_MCEP","FlowPro_VPtPOI_MCEP",iNbinsPt,dPtMin,dPtMax);
 200   fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
 201   fHistProDiffFlowPtPOI->SetYTitle("");
 202   fHistList->Add(fHistProDiffFlowPtPOI);
 203
 204   fHistProDiffFlowEtaPOI = new TProfile("FlowPro_VetaPOI_MCEP","FlowPro_VetaPOI_MCEP",iNbinsEta,dEtaMin,dEtaMax);
 205   fHistProDiffFlowEtaPOI->SetXTitle("#eta");
 206   fHistProDiffFlowEtaPOI->SetYTitle("");
 207   fHistList->Add(fHistProDiffFlowEtaPOI);
 208
 209   fHistSpreadOfFlow = new TH1D("fHistSpreadOfFlow","fHistSpreadOfFlow",1000,-1,1);
 210   fHistSpreadOfFlow->SetXTitle("v_{2}");
 211   fHistSpreadOfFlow->SetYTitle("counts");
 212   fHistList->Add(fHistSpreadOfFlow);
 213
 214   fEventNumber = 0;  //set number of events to zero
 215
 216   if(fEvaluateMixedHarmonics) this->BookObjectsForMixedHarmonics();
 217
 218   TH1::AddDirectory(oldHistAddStatus);
 219 }
 220
 221 //-----------------------------------------------------------------------
 222
 223 void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent) {
 224
 225   //Calculate v2 from the MC reaction plane
 226   if (anEvent) {
 227
 228     // get the MC reaction plane angle
 229     Double_t aRP = anEvent->GetMCReactionPlaneAngle();
 230     //fill control histograms
 231     fCommonHists->FillControlHistograms(anEvent);
 232
 233     //get the Q vector from the FlowEvent
 234     AliFlowVector vQ = anEvent->GetQ(fHarmonic);
 235     //cout<<"vQ.Mod() = " << vQ.Mod() << endl;
 236     //for chi calculation:
 237     *fQsum += vQ;
 238     //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
 239     fQ2sum += vQ.Mod2();
 240     //cout<<"fQ2sum = "<<fQ2sum<<endl;
 241
 242     fHistRP->Fill(aRP);
 243
 244     Double_t dPhi = 0.;
 245     Double_t dv  = 0.;
 246     Double_t dPt  = 0.;
 247     Double_t dEta = 0.;
 248     //Double_t dPi = TMath::Pi();
 249
 250     // profile to calculate flow e-b-y:
 251     TProfile *flowEBE = new TProfile("flowEBE","flowEBE",1,0,1);
 252
 253     //calculate flow
 254     //loop over the tracks of the event
 255     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 256     Int_t iNumberOfRPs = anEvent->GetEventNSelTracksRP();
 257     for (Int_t i=0;i<iNumberOfTracks;i++)
 258       {
 259         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 260         if (pTrack){
 261           if (pTrack->InRPSelection()){
 262             dPhi = pTrack->Phi();
 263             dv  = TMath::Cos(fHarmonic*(dPhi-aRP));
 264                  dPt  = pTrack->Pt();
 265                  dEta = pTrack->Eta();
 266             //reference flow:
 267             fHistProIntFlow->Fill(0.,dv);
 268             //reference flow versus multiplicity:
 269             fHistProIntFlowVsM->Fill(iNumberOfRPs+0.5,dv);
 270             //reference flow e-b-e:
 271             flowEBE->Fill(0.,dv);
 272             //differential flow (Pt, Eta, RP):
 273             fHistProDiffFlowPtEtaRP->Fill(dPt,dEta,dv,1.);
 274             //differential flow (Pt, RP):
 275             fHistProDiffFlowPtRP->Fill(dPt,dv,1.);
 276             //differential flow (Eta, RP):
 277             fHistProDiffFlowEtaRP->Fill(dEta,dv,1.);
 278           }
 279           if (pTrack->InPOISelection()) {
 280             dPhi = pTrack->Phi();
 281             //if (dPhi<0.) dPhi+=2*TMath::Pi();
 282             //calculate flow v2:
 283             dv  = TMath::Cos(fHarmonic*(dPhi-aRP));
 284             dPt  = pTrack->Pt();
 285             dEta = pTrack->Eta();
 286             //differential flow (Pt, Eta, POI):
 287             fHistProDiffFlowPtEtaPOI->Fill(dPt,dEta,dv,1.);
 288             //differential flow (Pt, POI):
 289             fHistProDiffFlowPtPOI->Fill(dPt,dv,1.);
 290             //differential flow (Eta, POI):
 291             fHistProDiffFlowEtaPOI->Fill(dEta,dv,1.);
 292           }
 293         }//track selected
 294       }//loop over tracks
 295
 296     fEventNumber++;
 297     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 298
 299     // store flow value for this event:
 300     fHistSpreadOfFlow->Fill(flowEBE->GetBinContent(1),flowEBE->GetBinEntries(1));
 301     delete flowEBE;
 302
 303     if(fEvaluateMixedHarmonics) EvaluateMixedHarmonics(anEvent);
 304   }
 305 }
 306   //--------------------------------------------------------------------
 307
 308 void AliFlowAnalysisWithMCEventPlane::GetOutputHistograms(TList *outputListHistos) {
 309  // get the pointers to all output histograms before calling Finish()
 310  if (outputListHistos) {
 311     //Get the common histograms from the output list
 312     AliFlowCommonHist *pCommonHists = dynamic_cast<AliFlowCommonHist*>
 313       (outputListHistos->FindObject("AliFlowCommonHistMCEP"));
 314     AliFlowCommonHistResults *pCommonHistResults =
 315       dynamic_cast<AliFlowCommonHistResults*>
 316       (outputListHistos->FindObject("AliFlowCommonHistResultsMCEP"));
 317
 318     TProfile *pHistProIntFlow = dynamic_cast<TProfile*>
 319       (outputListHistos->FindObject("FlowPro_V_MCEP"));
 320
 321     TProfile2D *pHistProDiffFlowPtEtaRP = dynamic_cast<TProfile2D*>
 322       (outputListHistos->FindObject("FlowPro_VPtEtaRP_MCEP"));
 323
 324     TProfile *pHistProDiffFlowPtRP = dynamic_cast<TProfile*>
 325       (outputListHistos->FindObject("FlowPro_VPtRP_MCEP"));
 326
 327     TProfile *pHistProDiffFlowEtaRP = dynamic_cast<TProfile*>
 328       (outputListHistos->FindObject("FlowPro_VetaRP_MCEP"));
 329
 330     TProfile2D *pHistProDiffFlowPtEtaPOI = dynamic_cast<TProfile2D*>
 331       (outputListHistos->FindObject("FlowPro_VPtEtaPOI_MCEP"));
 332
 333     TProfile *pHistProDiffFlowPtPOI = dynamic_cast<TProfile*>
 334       (outputListHistos->FindObject("FlowPro_VPtPOI_MCEP"));
 335
 336     TProfile *pHistProDiffFlowEtaPOI = dynamic_cast<TProfile*>
 337       (outputListHistos->FindObject("FlowPro_VetaPOI_MCEP"));
 338
 339     if (pCommonHists && pCommonHistResults && pHistProIntFlow &&
 340         pHistProDiffFlowPtRP && pHistProDiffFlowEtaRP &&
 341         pHistProDiffFlowPtPOI && pHistProDiffFlowEtaPOI) {
 342       this->SetCommonHists(pCommonHists);
 343       this->SetCommonHistsRes(pCommonHistResults);
 344       this->SetHistProIntFlow(pHistProIntFlow);
 345       this->SetHistProDiffFlowPtEtaRP(pHistProDiffFlowPtEtaRP);
 346       this->SetHistProDiffFlowPtRP(pHistProDiffFlowPtRP);
 347       this->SetHistProDiffFlowEtaRP(pHistProDiffFlowEtaRP);
 348       this->SetHistProDiffFlowPtEtaPOI(pHistProDiffFlowPtEtaPOI);
 349       this->SetHistProDiffFlowPtPOI(pHistProDiffFlowPtPOI);
 350       this->SetHistProDiffFlowEtaPOI(pHistProDiffFlowEtaPOI);
 351     } else {
 352       cout<<"WARNING: Histograms needed to run Finish() are not accessible!"<<endl;  }
 353
 354     //fListHistos->Print();
 355
 356    TList *pMixedHarmonicsList = dynamic_cast<TList*>
 357      (outputListHistos->FindObject("Mixed Harmonics"));
 358    if(pMixedHarmonicsList) {this->GetOutputHistoramsForMixedHarmonics(pMixedHarmonicsList);}
 359
 360   } else { cout << "histogram list pointer is empty" << endl;}
 361
 362 }
 363
 364 //--------------------------------------------------------------------
 365
 366 void AliFlowAnalysisWithMCEventPlane::GetOutputHistoramsForMixedHarmonics(TList *mixedHarmonicsList)
 367 {
 368  // Get pointers to all objects relevant for mixed harmonics study.
 369
 370  if(mixedHarmonicsList)
 371  {
 372   // ...
 373  } else
 374    {
 375     cout<<endl;
 376     cout<<"WARNING (MCEP): mixedHarmonicsList in NULL in MCEP::GetOutputHistoramsForMixedHarmonics() !!!! "<<endl;
 377     cout<<endl;
 378    }
 379
 380 } // end of void AliFlowAnalysisWithMCEventPlane::GetOutputHistoramsForMixedHarmonics(TList *mixedHarmonicsList)
 381
 382 //--------------------------------------------------------------------
 383
 384 void AliFlowAnalysisWithMCEventPlane::Finish() {
 385
 386   //*************make histograms etc.
 387   if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
 388
 389   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 390   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 391
 392   // access harmonic:
 393   if(fCommonHists && fCommonHists->GetHarmonic())
 394   {
 395    fHarmonic = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1); // to be improved (moved somewhere else?)
 396   }
 397
 398   //reference flow :
 399   Double_t dV = fHistProIntFlow->GetBinContent(1);
 400   Double_t dErrV = fHistProIntFlow->GetBinError(1); // to be improved (treatment of errors for non-Gaussian distribution needed!)
 401   //fill reference flow:
 402   fCommonHistsRes->FillIntegratedFlow(dV,dErrV);
 403   cout<<"dV"<<fHarmonic<<"{MC} is       "<<dV<<" +- "<<dErrV<<endl;
 404
 405   //RP:
 406   TH1F* fHistPtRP = NULL;
 407   if(fCommonHists && fCommonHists->GetHistPtRP())
 408   {
 409    fHistPtRP = fCommonHists->GetHistPtRP();
 410   }
 411   Double_t dYieldPtRP = 0.;
 412   Double_t dVRP = 0.;
 413   Double_t dErrVRP = 0.;
 414   Double_t dSumRP = 0.;
 415   //differential flow (RP, Pt):
 416   Double_t dvPtRP = 0.;
 417   Double_t dErrvPtRP = 0.;
 418   for(Int_t b=1;b<=iNbinsPt;b++)
 419   {
 420    dvPtRP    = fHistProDiffFlowPtRP->GetBinContent(b);
 421    dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 422    fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
 423    if(fHistPtRP){
 424         //integrated flow (RP)
 425         dYieldPtRP = fHistPtRP->GetBinContent(b);
 426         dVRP += dvPtRP*dYieldPtRP;
 427         dSumRP += dYieldPtRP;
 428         //error on integrated flow
 429         dErrVRP += dYieldPtRP*dYieldPtRP*dErrvPtRP*dErrvPtRP;
 430       }
 431   }
 432
 433   if (!(TMath::AreEqualAbs(dSumRP, 0.0, 1e-10)) ) {
 434     dVRP /= dSumRP;  //because pt distribution should be normalised
 435     dErrVRP /= (dSumRP*dSumRP);
 436     dErrVRP = TMath::Sqrt(dErrVRP);
 437   }
 438   // fill integrated flow (RP):
 439   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
 440   cout<<"dV"<<fHarmonic<<"{MC} (RP) is  "<<dVRP<<" +- "<<dErrVRP<<endl;
 441
 442   //differential flow (RP, Eta):
 443   Double_t dvEtaRP = 0.;
 444   Double_t dErrvEtaRP = 0.;
 445   for(Int_t b=1;b<=iNbinsEta;b++)
 446   {
 447    dvEtaRP    = fHistProDiffFlowEtaRP->GetBinContent(b);
 448    dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 449    fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
 450   }
 451
 452   //POI:
 453   TH1F* fHistPtPOI = NULL;
 454   if(fCommonHists && fCommonHists->GetHistPtPOI())
 455   {
 456    fHistPtPOI = fCommonHists->GetHistPtPOI();
 457   }
 458   Double_t dYieldPtPOI = 0.;
 459   Double_t dVPOI = 0.;
 460   Double_t dErrVPOI = 0.;
 461   Double_t dSumPOI = 0.;
 462   Double_t dvproPtPOI = 0.;
 463   Double_t dErrdifcombPtPOI = 0.;
 464   Double_t dvproEtaPOI = 0.;
 465   Double_t dErrdifcombEtaPOI = 0.;
 466   //Pt:
 467   if(fHistProDiffFlowPtPOI) {
 468     for(Int_t b=1;b<=iNbinsPt;b++){
 469       dvproPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
 470       dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 471       //fill TH1D
 472       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dvproPtPOI, dErrdifcombPtPOI);
 473       if (fHistPtPOI){
 474         //integrated flow (POI)
 475         dYieldPtPOI = fHistPtPOI->GetBinContent(b);
 476         dVPOI += dvproPtPOI*dYieldPtPOI;
 477         dSumPOI += dYieldPtPOI;
 478         //error on integrated flow
 479         dErrVPOI += dYieldPtPOI*dYieldPtPOI*dErrdifcombPtPOI*dErrdifcombPtPOI;
 480       }
 481     }//end of for(Int_t b=0;b<iNbinsPt;b++)
 482   } else { cout<<"fHistProFlow is NULL"<<endl; }
 483
 484   if (!(TMath::AreEqualAbs(dSumPOI, 0.0, 1e-10)) ) {
 485     dVPOI /= dSumPOI;  //because pt distribution should be normalised
 486     dErrVPOI /= (dSumPOI*dSumPOI);
 487     dErrVPOI = TMath::Sqrt(dErrVPOI);
 488   }
 489   cout<<"dV"<<fHarmonic<<"{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;
 490
 491   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
 492
 493   //Eta:
 494   if(fHistProDiffFlowEtaPOI)
 495   {
 496    for(Int_t b=1;b<=iNbinsEta;b++)
 497    {
 498     dvproEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
 499     dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 500     //fill common hist results:
 501     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dvproEtaPOI, dErrdifcombEtaPOI);
 502    }
 503   }
 504
 505   cout<<endl;
 506   //cout<<".....finished"<<endl;
 507 }
 508
 509 //-----------------------------------------------------------------------
 510
 511 void AliFlowAnalysisWithMCEventPlane::InitalizeArraysForMixedHarmonics()
 512 {
 513  // Iinitialize all arrays for mixed harmonics.
 514
 515  for(Int_t cs=0;cs<2;cs++) // cos/sin
 516  {
 517   fPairCorrelator[cs] = NULL;
 518   fPairCorrelatorVsM[cs] = NULL;
 519   for(Int_t sd=0;sd<2;sd++) // pt sum/difference
 520   {
 521    fPairCorrelatorVsPtSumDiff[cs][sd] = NULL;
 522   }
 523  }
 524
 525 } // end of void InitalizeArraysForMixedHarmonics()
 526
 527 //-----------------------------------------------------------------------
 528
 529 void AliFlowAnalysisWithMCEventPlane::BookObjectsForMixedHarmonics()
 530 {
 531  // Book all objects needed for mixed harmonics.
 532
 533  // List holding all objects relevant for mixed harmonics:
 534  fMixedHarmonicsList->SetName("Mixed Harmonics");
 535  fMixedHarmonicsList->SetOwner(kTRUE);
 536  fHistList->Add(fMixedHarmonicsList);
 537
 538  // Profile holding settings relevant for mixed harmonics:
 539  TString mixedHarmonicsSettingsName = "fMixedHarmonicsSettings";
 540  fMixedHarmonicsSettings = new TProfile(mixedHarmonicsSettingsName.Data(),"Settings for Mixed Harmonics",4,0,4);
 541  //fMixedHarmonicsSettings->GetXaxis()->SetLabelSize(0.025);
 542  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(1,"fEvaluateMixedHarmonics");
 543  fMixedHarmonicsSettings->Fill(0.5,(Int_t)fEvaluateMixedHarmonics);
 544  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(2,"fNinCorrelator");
 545  fMixedHarmonicsSettings->Fill(1.5,(Int_t)fNinCorrelator);
 546  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(3,"fMinCorrelator");
 547  fMixedHarmonicsSettings->Fill(2.5,(Int_t)fMinCorrelator);
 548  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(4,"x in #phi_{pair}"); // phi_{pair} = x*phi1+(1-x)*phi2
 549  fMixedHarmonicsSettings->Fill(3.5,fXinPairAngle);
 550  fMixedHarmonicsList->Add(fMixedHarmonicsSettings);
 551
 552  // Profiles used to calculate <cos[m*phi_{pair}-n*RP]> and <sin[m*phi_{pair}-n*RP]>, where phi_{pair} = x*phi1+(1-x)*phi2:
 553  TString cosSinFlag[2] = {"Cos","Sin"};
 554  TString cosSinTitleFlag[2] = {"cos[m#phi_{pair}-n#psi_{RP}]","sin[m#phi_{pair}-n#psi_{RP}]"};
 555  if(fNinCorrelator == 2 && fMinCorrelator == 2 && TMath::Abs(fXinPairAngle-0.5)<1.e-44) // default values
 556  {
 557   cosSinTitleFlag[0] = "cos[#phi_{1}+#phi_{2}-2#psi_{RP})]";
 558   cosSinTitleFlag[1] = "sin[#phi_{1}+#phi_{2}-2#psi_{RP})]";
 559  }
 560  TString psdFlag[2] = {"Sum","Diff"};
 561  TString psdTitleFlag[2] = {"(p_{t,1}+p_{t,2})/2","#left|p_{t,1}-p_{t,2}#right|"};
 562  TString pairCorrelatorName = "fPairCorrelator";
 563  TString pairCorrelatorVsMName = "fPairCorrelatorVsM";
 564  TString pairCorrelatorVsPtSumDiffName = "fPairCorrelatorVsPt";
 565  Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 566  Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 567  Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 568  for(Int_t cs=0;cs<2;cs++)
 569  {
 570   fPairCorrelator[cs] = new TProfile(Form("%s, %s",pairCorrelatorName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),1,0.,1.);
 571   fPairCorrelator[cs]->GetXaxis()->SetBinLabel(1,cosSinTitleFlag[cs].Data());
 572   fMixedHarmonicsList->Add(fPairCorrelator[cs]);
 573
 574   fPairCorrelatorVsM[cs] = new TProfile(Form("%s, %s",pairCorrelatorVsMName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),fnBinsMult,fMinMult,fMaxMult);
 575   fPairCorrelatorVsM[cs]->GetXaxis()->SetTitle("# of RPs");
 576   fMixedHarmonicsList->Add(fPairCorrelatorVsM[cs]);
 577
 578   for(Int_t sd=0;sd<2;sd++)
 579   {
 580    fPairCorrelatorVsPtSumDiff[cs][sd] = new TProfile(Form("%s%s, %s",pairCorrelatorVsPtSumDiffName.Data(),psdFlag[sd].Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),iNbinsPt,dPtMin,dPtMax);
 581    fPairCorrelatorVsPtSumDiff[cs][sd]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
 582    fMixedHarmonicsList->Add(fPairCorrelatorVsPtSumDiff[cs][sd]);
 583   } // end of for(Int_t sd=0;sd<2;sd++)
 584  } // end of for(Int_t cs=0;cs<2;cs++)
 585
 586 } // end of void AliFlowAnalysisWithMCEventPlane::BookObjectsForMixedHarmonics()
 587
 588 //-----------------------------------------------------------------------
 589
 590 void AliFlowAnalysisWithMCEventPlane::EvaluateMixedHarmonics(AliFlowEventSimple* anEvent)
 591 {
 592  // Evaluate correlators relevant for the mixed harmonics.
 593
 594  // Get the MC reaction plane angle:
 595  Double_t dReactionPlane = anEvent->GetMCReactionPlaneAngle();
 596  // Get the number of tracks:
 597  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 598  Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of Reference Particles
 599  AliFlowTrackSimple *pTrack = NULL;
 600  Double_t dPhi1 = 0.;
 601  Double_t dPhi2 = 0.;
 602  Double_t dPt1 = 0.;
 603  Double_t dPt2 = 0.;
 604  Double_t n = fNinCorrelator; // shortcut
 605  Double_t m = fMinCorrelator; // shortcut
 606  Double_t x = fXinPairAngle; // shortcut
 607  for(Int_t i=0;i<iNumberOfTracks;i++)
 608  {
 609   pTrack = anEvent->GetTrack(i);
 610   if(pTrack && pTrack->InRPSelection())
 611   {
 612    dPhi1 = pTrack->Phi();
 613    dPt1 = pTrack->Pt();
 614   }
 615   for(Int_t j=0;j<iNumberOfTracks;j++)
 616   {
 617    if(j==i) continue;
 618    pTrack = anEvent->GetTrack(j);
 619    if(pTrack && pTrack->InPOISelection())
 620    {
 621     dPhi2 = pTrack->Phi();
 622     dPt2 = pTrack->Pt();
 623    }
 624    Double_t dPhiPair = x*dPhi1+(1.-x)*dPhi2;
 625    Double_t dPtSum = 0.5*(dPt1+dPt2);
 626    Double_t dPtDiff = TMath::Abs(dPt1-dPt2);
 627    fPairCorrelator[0]->Fill(0.5,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 628    fPairCorrelator[1]->Fill(0.5,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 629    fPairCorrelatorVsM[0]->Fill(nRP+0.5,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 630    fPairCorrelatorVsM[1]->Fill(nRP+0.5,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 631    fPairCorrelatorVsPtSumDiff[0][0]->Fill(dPtSum,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 632    fPairCorrelatorVsPtSumDiff[1][0]->Fill(dPtSum,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 633    fPairCorrelatorVsPtSumDiff[0][1]->Fill(dPtDiff,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 634    fPairCorrelatorVsPtSumDiff[1][1]->Fill(dPtDiff,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 635   } // end of for(Int_t j=i+1;j<iNumberOfTracks;j++)
 636  } // end of for(Int_t i=0;i<iNumberOfTracks;i++)
 637
 638 } // end of void AliFlowAnalysisWithMCEventPlane::EvaluateMixedHarmonics()
 639
 640
 641