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