PWG2/FLOW/AliFlowCommon/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(kTRUE),
  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");
 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 = fCommonHists->GetHistPtRP();
 405   Double_t dYieldPtRP = 0.;
 406   Double_t dVRP = 0.;
 407   Double_t dErrVRP = 0.;
 408   Double_t dSumRP = 0.;
 409   //differential flow (RP, Pt):
 410   Double_t dvPtRP = 0.;
 411   Double_t dErrvPtRP = 0.;
 412   for(Int_t b=1;b<=iNbinsPt;b++)
 413   {
 414    dvPtRP    = fHistProDiffFlowPtRP->GetBinContent(b);
 415    dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 416    fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
 417    if(fHistPtRP){
 418         //integrated flow (RP)
 419         dYieldPtRP = fHistPtRP->GetBinContent(b);
 420         dVRP += dvPtRP*dYieldPtRP;
 421         dSumRP += dYieldPtRP;
 422         //error on integrated flow
 423         dErrVRP += dYieldPtRP*dYieldPtRP*dErrvPtRP*dErrvPtRP;
 424       }
 425   }
 426
 427   if (!(TMath::AreEqualAbs(dSumRP, 0.0, 1e-10)) ) {
 428     dVRP /= dSumRP;  //because pt distribution should be normalised
 429     dErrVRP /= (dSumRP*dSumRP);
 430     dErrVRP = TMath::Sqrt(dErrVRP);
 431   }
 432   // fill integrated flow (RP):
 433   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
 434   cout<<"dV"<<fHarmonic<<"{MC} (RP) is  "<<dVRP<<" +- "<<dErrVRP<<endl;
 435
 436   //differential flow (RP, Eta):
 437   Double_t dvEtaRP = 0.;
 438   Double_t dErrvEtaRP = 0.;
 439   for(Int_t b=1;b<=iNbinsEta;b++)
 440   {
 441    dvEtaRP    = fHistProDiffFlowEtaRP->GetBinContent(b);
 442    dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 443    fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
 444   }
 445
 446   //POI:
 447   TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI();
 448   Double_t dYieldPtPOI = 0.;
 449   Double_t dVPOI = 0.;
 450   Double_t dErrVPOI = 0.;
 451   Double_t dSumPOI = 0.;
 452   Double_t dvproPtPOI = 0.;
 453   Double_t dErrdifcombPtPOI = 0.;
 454   Double_t dvproEtaPOI = 0.;
 455   Double_t dErrdifcombEtaPOI = 0.;
 456   //Pt:
 457   if(fHistProDiffFlowPtPOI) {
 458     for(Int_t b=1;b<=iNbinsPt;b++){
 459       dvproPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
 460       dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 461       //fill TH1D
 462       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dvproPtPOI, dErrdifcombPtPOI);
 463       if (fHistPtPOI){
 464         //integrated flow (POI)
 465         dYieldPtPOI = fHistPtPOI->GetBinContent(b);
 466         dVPOI += dvproPtPOI*dYieldPtPOI;
 467         dSumPOI += dYieldPtPOI;
 468         //error on integrated flow
 469         dErrVPOI += dYieldPtPOI*dYieldPtPOI*dErrdifcombPtPOI*dErrdifcombPtPOI;
 470       }
 471     }//end of for(Int_t b=0;b<iNbinsPt;b++)
 472   } else { cout<<"fHistProFlow is NULL"<<endl; }
 473
 474   if (!(TMath::AreEqualAbs(dSumPOI, 0.0, 1e-10)) ) {
 475     dVPOI /= dSumPOI;  //because pt distribution should be normalised
 476     dErrVPOI /= (dSumPOI*dSumPOI);
 477     dErrVPOI = TMath::Sqrt(dErrVPOI);
 478   }
 479   cout<<"dV"<<fHarmonic<<"{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;
 480
 481   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
 482
 483   //Eta:
 484   if(fHistProDiffFlowEtaPOI)
 485   {
 486    for(Int_t b=1;b<=iNbinsEta;b++)
 487    {
 488     dvproEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
 489     dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 490     //fill common hist results:
 491     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dvproEtaPOI, dErrdifcombEtaPOI);
 492    }
 493   }
 494
 495   cout<<endl;
 496   //cout<<".....finished"<<endl;
 497 }
 498
 499 //-----------------------------------------------------------------------
 500
 501 void AliFlowAnalysisWithMCEventPlane::InitalizeArraysForMixedHarmonics()
 502 {
 503  // Iinitialize all arrays for mixed harmonics.
 504
 505  for(Int_t cs=0;cs<2;cs++) // cos/sin
 506  {
 507   fPairCorrelator[cs] = NULL;
 508   fPairCorrelatorVsM[cs] = NULL;
 509   for(Int_t sd=0;sd<2;sd++) // pt sum/difference
 510   {
 511    fPairCorrelatorVsPtSumDiff[cs][sd] = NULL;
 512   }
 513  }
 514
 515 } // end of void InitalizeArraysForMixedHarmonics()
 516
 517 //-----------------------------------------------------------------------
 518
 519 void AliFlowAnalysisWithMCEventPlane::BookObjectsForMixedHarmonics()
 520 {
 521  // Book all objects needed for mixed harmonics.
 522
 523  // List holding all objects relevant for mixed harmonics:
 524  fMixedHarmonicsList->SetName("Mixed Harmonics");
 525  fMixedHarmonicsList->SetOwner(kTRUE);
 526  fHistList->Add(fMixedHarmonicsList);
 527
 528  // Profile holding settings relevant for mixed harmonics:
 529  TString mixedHarmonicsSettingsName = "fMixedHarmonicsSettings";
 530  fMixedHarmonicsSettings = new TProfile(mixedHarmonicsSettingsName.Data(),"Settings for Mixed Harmonics",4,0,4);
 531  //fMixedHarmonicsSettings->GetXaxis()->SetLabelSize(0.025);
 532  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(1,"fEvaluateMixedHarmonics");
 533  fMixedHarmonicsSettings->Fill(0.5,(Int_t)fEvaluateMixedHarmonics);
 534  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(2,"fNinCorrelator");
 535  fMixedHarmonicsSettings->Fill(1.5,(Int_t)fNinCorrelator);
 536  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(3,"fMinCorrelator");
 537  fMixedHarmonicsSettings->Fill(2.5,(Int_t)fMinCorrelator);
 538  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(4,"x in #phi_{pair}"); // phi_{pair} = x*phi1+(1-x)*phi2
 539  fMixedHarmonicsSettings->Fill(3.5,fXinPairAngle);
 540  fMixedHarmonicsList->Add(fMixedHarmonicsSettings);
 541
 542  // 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:
 543  TString cosSinFlag[2] = {"Cos","Sin"};
 544  TString cosSinTitleFlag[2] = {"cos[m#phi_{pair}-n#psi_{RP}]","sin[m#phi_{pair}-n#psi_{RP}]"};
 545  if(fNinCorrelator == 2 && fMinCorrelator == 2 && TMath::Abs(fXinPairAngle-0.5)<1.e-44) // default values
 546  {
 547   cosSinTitleFlag[0] = "cos[#phi_{1}+#phi_{2}-2#psi_{RP})]";
 548   cosSinTitleFlag[1] = "sin[#phi_{1}+#phi_{2}-2#psi_{RP})]";
 549  }
 550  TString psdFlag[2] = {"Sum","Diff"};
 551  TString psdTitleFlag[2] = {"(p_{t,1}+p_{t,2})/2","#left|p_{t,1}-p_{t,2}#right|"};
 552  TString pairCorrelatorName = "fPairCorrelator";
 553  TString pairCorrelatorVsMName = "fPairCorrelatorVsM";
 554  TString pairCorrelatorVsPtSumDiffName = "fPairCorrelatorVsPt";
 555  Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 556  Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 557  Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 558  for(Int_t cs=0;cs<2;cs++)
 559  {
 560   fPairCorrelator[cs] = new TProfile(Form("%s, %s",pairCorrelatorName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),1,0.,1.);
 561   fPairCorrelator[cs]->GetXaxis()->SetBinLabel(1,cosSinTitleFlag[cs].Data());
 562   fMixedHarmonicsList->Add(fPairCorrelator[cs]);
 563
 564   fPairCorrelatorVsM[cs] = new TProfile(Form("%s, %s",pairCorrelatorVsMName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),fnBinsMult,fMinMult,fMaxMult);
 565   fPairCorrelatorVsM[cs]->GetXaxis()->SetTitle("# of RPs");
 566   fMixedHarmonicsList->Add(fPairCorrelatorVsM[cs]);
 567
 568   for(Int_t sd=0;sd<2;sd++)
 569   {
 570    fPairCorrelatorVsPtSumDiff[cs][sd] = new TProfile(Form("%s%s, %s",pairCorrelatorVsPtSumDiffName.Data(),psdFlag[sd].Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),iNbinsPt,dPtMin,dPtMax);
 571    fPairCorrelatorVsPtSumDiff[cs][sd]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
 572    fMixedHarmonicsList->Add(fPairCorrelatorVsPtSumDiff[cs][sd]);
 573   } // end of for(Int_t sd=0;sd<2;sd++)
 574  } // end of for(Int_t cs=0;cs<2;cs++)
 575
 576 } // end of void AliFlowAnalysisWithMCEventPlane::BookObjectsForMixedHarmonics()
 577
 578 //-----------------------------------------------------------------------
 579
 580 void AliFlowAnalysisWithMCEventPlane::EvaluateMixedHarmonics(AliFlowEventSimple* anEvent)
 581 {
 582  // Evaluate correlators relevant for the mixed harmonics.
 583
 584  // Get the MC reaction plane angle:
 585  Double_t dReactionPlane = anEvent->GetMCReactionPlaneAngle();
 586  // Get the number of tracks:
 587  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 588  Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of Reference Particles
 589  AliFlowTrackSimple *pTrack = NULL;
 590  Double_t dPhi1 = 0.;
 591  Double_t dPhi2 = 0.;
 592  Double_t dPt1 = 0.;
 593  Double_t dPt2 = 0.;
 594  Double_t n = fNinCorrelator; // shortcut
 595  Double_t m = fMinCorrelator; // shortcut
 596  Double_t x = fXinPairAngle; // shortcut
 597  for(Int_t i=0;i<iNumberOfTracks;i++)
 598  {
 599   pTrack = anEvent->GetTrack(i);
 600   if(pTrack && pTrack->InRPSelection())
 601   {
 602    dPhi1 = pTrack->Phi();
 603    dPt1 = pTrack->Pt();
 604   }
 605   for(Int_t j=0;j<iNumberOfTracks;j++)
 606   {
 607    if(j==i) continue;
 608    pTrack = anEvent->GetTrack(j);
 609    if(pTrack && pTrack->InPOISelection())
 610    {
 611     dPhi2 = pTrack->Phi();
 612     dPt2 = pTrack->Pt();
 613    }
 614    Double_t dPhiPair = x*dPhi1+(1.-x)*dPhi2;
 615    Double_t dPtSum = 0.5*(dPt1+dPt2);
 616    Double_t dPtDiff = TMath::Abs(dPt1-dPt2);
 617    fPairCorrelator[0]->Fill(0.5,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 618    fPairCorrelator[1]->Fill(0.5,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 619    fPairCorrelatorVsM[0]->Fill(nRP+0.5,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 620    fPairCorrelatorVsM[1]->Fill(nRP+0.5,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 621    fPairCorrelatorVsPtSumDiff[0][0]->Fill(dPtSum,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 622    fPairCorrelatorVsPtSumDiff[1][0]->Fill(dPtSum,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 623    fPairCorrelatorVsPtSumDiff[0][1]->Fill(dPtDiff,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
 624    fPairCorrelatorVsPtSumDiff[1][1]->Fill(dPtDiff,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
 625   } // end of for(Int_t j=i+1;j<iNumberOfTracks;j++)
 626  } // end of for(Int_t i=0;i<iNumberOfTracks;i++)
 627
 628 } // end of void AliFlowAnalysisWithMCEventPlane::EvaluateMixedHarmonics()
 629
 630
 631