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    fHistProDiffFlowPtEtaRP(NULL),
  57    fHistProDiffFlowPtRP(NULL),
  58    fHistProDiffFlowEtaRP(NULL),
  59    fHistProDiffFlowPtEtaPOI(NULL),
  60    fHistProDiffFlowPtPOI(NULL),
  61    fHistProDiffFlowEtaPOI(NULL),
  62    fHistSpreadOfFlow(NULL),
  63    fHarmonic(2),
  64    fResonanceList(NULL),
  65    fFlowOfResonances(kFALSE),
  66    fResonanceSettings(NULL),
  67    fXinPairAngle(0.5)
  68 {
  69
  70   // Constructor.
  71   fHistList = new TList();
  72
  73   fQsum = new TVector2;        // flow vector sum
  74
  75   fResonanceList = new TList();
  76   // Initialize arrays needed for the flow of resonances:
  77   for(Int_t cs=0;cs<2;cs++) // cos/sin
  78   {
  79    fPairCorrelator[cs] = NULL;
  80   }
  81
  82 }
  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}{2}");
 168   fHistProIntFlow->SetYTitle("");
 169   fHistList->Add(fHistProIntFlow);
 170
 171   fHistProDiffFlowPtEtaRP = new TProfile2D("FlowPro_VPtEtaRP_MCEP","FlowPro_VPtEtaRP_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
 172   fHistProDiffFlowPtEtaRP->SetXTitle("P_{t}");
 173   fHistProDiffFlowPtEtaRP->SetYTitle("#eta");
 174   fHistList->Add(fHistProDiffFlowPtEtaRP);
 175
 176   fHistProDiffFlowPtRP = new TProfile("FlowPro_VPtRP_MCEP","FlowPro_VPtRP_MCEP",iNbinsPt,dPtMin,dPtMax);
 177   fHistProDiffFlowPtRP->SetXTitle("P_{t}");
 178   fHistProDiffFlowPtRP->SetYTitle("");
 179   fHistList->Add(fHistProDiffFlowPtRP);
 180
 181   fHistProDiffFlowEtaRP = new TProfile("FlowPro_VetaRP_MCEP","FlowPro_VetaRP_MCEP",iNbinsEta,dEtaMin,dEtaMax);
 182   fHistProDiffFlowEtaRP->SetXTitle("#eta");
 183   fHistProDiffFlowEtaRP->SetYTitle("");
 184   fHistList->Add(fHistProDiffFlowEtaRP);
 185
 186   fHistProDiffFlowPtEtaPOI = new TProfile2D("FlowPro_VPtEtaPOI_MCEP","FlowPro_VPtEtaPOI_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
 187   fHistProDiffFlowPtEtaPOI->SetXTitle("P_{t}");
 188   fHistProDiffFlowPtEtaPOI->SetYTitle("#eta");
 189   fHistList->Add(fHistProDiffFlowPtEtaPOI);
 190
 191   fHistProDiffFlowPtPOI = new TProfile("FlowPro_VPtPOI_MCEP","FlowPro_VPtPOI_MCEP",iNbinsPt,dPtMin,dPtMax);
 192   fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
 193   fHistProDiffFlowPtPOI->SetYTitle("");
 194   fHistList->Add(fHistProDiffFlowPtPOI);
 195
 196   fHistProDiffFlowEtaPOI = new TProfile("FlowPro_VetaPOI_MCEP","FlowPro_VetaPOI_MCEP",iNbinsEta,dEtaMin,dEtaMax);
 197   fHistProDiffFlowEtaPOI->SetXTitle("#eta");
 198   fHistProDiffFlowEtaPOI->SetYTitle("");
 199   fHistList->Add(fHistProDiffFlowEtaPOI);
 200
 201   fHistSpreadOfFlow = new TH1D("fHistSpreadOfFlow","fHistSpreadOfFlow",1000,-1,1);
 202   fHistSpreadOfFlow->SetXTitle("v_{2}");
 203   fHistSpreadOfFlow->SetYTitle("counts");
 204   fHistList->Add(fHistSpreadOfFlow);
 205
 206   fEventNumber = 0;  //set number of events to zero
 207
 208   if(fFlowOfResonances) this->BookObjectsForFlowOfResonances();
 209
 210   TH1::AddDirectory(oldHistAddStatus);
 211 }
 212
 213 //-----------------------------------------------------------------------
 214
 215 void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent) {
 216
 217   //Calculate v2 from the MC reaction plane
 218   if (anEvent) {
 219
 220     // get the MC reaction plane angle
 221     Double_t aRP = anEvent->GetMCReactionPlaneAngle();
 222     //fill control histograms
 223     fCommonHists->FillControlHistograms(anEvent);
 224
 225     //get the Q vector from the FlowEvent
 226     AliFlowVector vQ = anEvent->GetQ(fHarmonic);
 227     //cout<<"vQ.Mod() = " << vQ.Mod() << endl;
 228     //for chi calculation:
 229     *fQsum += vQ;
 230     //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
 231     fQ2sum += vQ.Mod2();
 232     //cout<<"fQ2sum = "<<fQ2sum<<endl;
 233
 234     fHistRP->Fill(aRP);
 235
 236     Double_t dPhi = 0.;
 237     Double_t dv  = 0.;
 238     Double_t dPt  = 0.;
 239     Double_t dEta = 0.;
 240     //Double_t dPi = TMath::Pi();
 241
 242     // profile to calculate flow e-b-y:
 243     TProfile *flowEBE = new TProfile("flowEBE","flowEBE",1,0,1);
 244
 245     //calculate flow
 246     //loop over the tracks of the event
 247     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 248     for (Int_t i=0;i<iNumberOfTracks;i++)
 249       {
 250         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 251         if (pTrack){
 252           if (pTrack->InRPSelection()){
 253             dPhi = pTrack->Phi();
 254             dv  = TMath::Cos(fHarmonic*(dPhi-aRP));
 255                  dPt  = pTrack->Pt();
 256                  dEta = pTrack->Eta();
 257             //no-name int. flow (to be improved = name needed!):
 258             fHistProIntFlow->Fill(0.,dv);
 259             //no-name int. flow e-b-e (to be improved = name needed!):
 260             flowEBE->Fill(0.,dv);
 261             //differential flow (Pt, Eta, RP):
 262             fHistProDiffFlowPtEtaRP->Fill(dPt,dEta,dv,1.);
 263             //differential flow (Pt, RP):
 264             fHistProDiffFlowPtRP->Fill(dPt,dv,1.);
 265             //differential flow (Eta, RP):
 266             fHistProDiffFlowEtaRP->Fill(dEta,dv,1.);
 267           }
 268           if (pTrack->InPOISelection()) {
 269             dPhi = pTrack->Phi();
 270             //if (dPhi<0.) dPhi+=2*TMath::Pi();
 271             //calculate flow v2:
 272             dv  = TMath::Cos(fHarmonic*(dPhi-aRP));
 273             dPt  = pTrack->Pt();
 274             dEta = pTrack->Eta();
 275             //differential flow (Pt, Eta, POI):
 276             fHistProDiffFlowPtEtaPOI->Fill(dPt,dEta,dv,1.);
 277             //differential flow (Pt, POI):
 278             fHistProDiffFlowPtPOI->Fill(dPt,dv,1.);
 279             //differential flow (Eta, POI):
 280             fHistProDiffFlowEtaPOI->Fill(dEta,dv,1.);
 281           }
 282         }//track selected
 283       }//loop over tracks
 284
 285     fEventNumber++;
 286     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 287
 288     // store flow value for this event:
 289     fHistSpreadOfFlow->Fill(flowEBE->GetBinContent(1),flowEBE->GetBinEntries(1));
 290     delete flowEBE;
 291
 292     if(fFlowOfResonances) FlowOfResonances(anEvent);
 293   }
 294 }
 295   //--------------------------------------------------------------------
 296
 297 void AliFlowAnalysisWithMCEventPlane::GetOutputHistograms(TList *outputListHistos) {
 298  // get the pointers to all output histograms before calling Finish()
 299  if (outputListHistos) {
 300     //Get the common histograms from the output list
 301     AliFlowCommonHist *pCommonHists = dynamic_cast<AliFlowCommonHist*>
 302       (outputListHistos->FindObject("AliFlowCommonHistMCEP"));
 303     AliFlowCommonHistResults *pCommonHistResults =
 304       dynamic_cast<AliFlowCommonHistResults*>
 305       (outputListHistos->FindObject("AliFlowCommonHistResultsMCEP"));
 306
 307     TProfile *pHistProIntFlow = dynamic_cast<TProfile*>
 308       (outputListHistos->FindObject("FlowPro_V_MCEP"));
 309
 310     TProfile2D *pHistProDiffFlowPtEtaRP = dynamic_cast<TProfile2D*>
 311       (outputListHistos->FindObject("FlowPro_VPtEtaRP_MCEP"));
 312
 313     TProfile *pHistProDiffFlowPtRP = dynamic_cast<TProfile*>
 314       (outputListHistos->FindObject("FlowPro_VPtRP_MCEP"));
 315
 316     TProfile *pHistProDiffFlowEtaRP = dynamic_cast<TProfile*>
 317       (outputListHistos->FindObject("FlowPro_VetaRP_MCEP"));
 318
 319     TProfile2D *pHistProDiffFlowPtEtaPOI = dynamic_cast<TProfile2D*>
 320       (outputListHistos->FindObject("FlowPro_VPtEtaPOI_MCEP"));
 321
 322     TProfile *pHistProDiffFlowPtPOI = dynamic_cast<TProfile*>
 323       (outputListHistos->FindObject("FlowPro_VPtPOI_MCEP"));
 324
 325     TProfile *pHistProDiffFlowEtaPOI = dynamic_cast<TProfile*>
 326       (outputListHistos->FindObject("FlowPro_VetaPOI_MCEP"));
 327
 328     if (pCommonHists && pCommonHistResults && pHistProIntFlow &&
 329         pHistProDiffFlowPtRP && pHistProDiffFlowEtaRP &&
 330         pHistProDiffFlowPtPOI && pHistProDiffFlowEtaPOI) {
 331       this->SetCommonHists(pCommonHists);
 332       this->SetCommonHistsRes(pCommonHistResults);
 333       this->SetHistProIntFlow(pHistProIntFlow);
 334       this->SetHistProDiffFlowPtEtaRP(pHistProDiffFlowPtEtaRP);
 335       this->SetHistProDiffFlowPtRP(pHistProDiffFlowPtRP);
 336       this->SetHistProDiffFlowEtaRP(pHistProDiffFlowEtaRP);
 337       this->SetHistProDiffFlowPtEtaPOI(pHistProDiffFlowPtEtaPOI);
 338       this->SetHistProDiffFlowPtPOI(pHistProDiffFlowPtPOI);
 339       this->SetHistProDiffFlowEtaPOI(pHistProDiffFlowEtaPOI);
 340     } else {
 341       cout<<"WARNING: Histograms needed to run Finish() are not accessible!"<<endl;  }
 342
 343     //fListHistos->Print();
 344   } else { cout << "histogram list pointer is empty" << endl;}
 345
 346 }
 347
 348   //--------------------------------------------------------------------
 349 void AliFlowAnalysisWithMCEventPlane::Finish() {
 350
 351   //*************make histograms etc.
 352   if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
 353
 354   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 355   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 356
 357   // access harmonic:
 358   if(fCommonHists && fCommonHists->GetHarmonic())
 359   {
 360    fHarmonic = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1); // to be improved (moved somewhere else?)
 361   }
 362
 363   // no-name int. flow (to be improved):
 364   Double_t dV = fHistProIntFlow->GetBinContent(1);
 365   Double_t dErrV = fHistProIntFlow->GetBinError(1); // to be improved (treatment of errors for non-Gaussian distribution needed!)
 366   // fill no-name int. flow (to be improved):
 367   fCommonHistsRes->FillIntegratedFlow(dV,dErrV);
 368   cout<<"dV"<<fHarmonic<<"{MC} is       "<<dV<<" +- "<<dErrV<<endl;
 369
 370   //RP:
 371   TH1F* fHistPtRP = fCommonHists->GetHistPtRP();
 372   Double_t dYieldPtRP = 0.;
 373   Double_t dVRP = 0.;
 374   Double_t dErrVRP = 0.;
 375   Double_t dSumRP = 0.;
 376   //differential flow (RP, Pt):
 377   Double_t dvPtRP = 0.;
 378   Double_t dErrvPtRP = 0.;
 379   for(Int_t b=1;b<=iNbinsPt;b++)
 380   {
 381    dvPtRP    = fHistProDiffFlowPtRP->GetBinContent(b);
 382    dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 383    fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
 384    if(fHistPtRP){
 385         //integrated flow (RP)
 386         dYieldPtRP = fHistPtRP->GetBinContent(b);
 387         dVRP += dvPtRP*dYieldPtRP;
 388         dSumRP += dYieldPtRP;
 389         //error on integrated flow
 390         dErrVRP += dYieldPtRP*dYieldPtRP*dErrvPtRP*dErrvPtRP;
 391       }
 392   }
 393
 394   if (!(TMath::AreEqualAbs(dSumRP, 0.0, 1e-10)) ) {
 395     dVRP /= dSumRP;  //because pt distribution should be normalised
 396     dErrVRP /= (dSumRP*dSumRP);
 397     dErrVRP = TMath::Sqrt(dErrVRP);
 398   }
 399   // fill integrated flow (RP):
 400   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
 401   cout<<"dV"<<fHarmonic<<"{MC} (RP) is  "<<dVRP<<" +- "<<dErrVRP<<endl;
 402
 403   //differential flow (RP, Eta):
 404   Double_t dvEtaRP = 0.;
 405   Double_t dErrvEtaRP = 0.;
 406   for(Int_t b=1;b<=iNbinsEta;b++)
 407   {
 408    dvEtaRP    = fHistProDiffFlowEtaRP->GetBinContent(b);
 409    dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 410    fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
 411   }
 412
 413   //POI:
 414   TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI();
 415   Double_t dYieldPtPOI = 0.;
 416   Double_t dVPOI = 0.;
 417   Double_t dErrVPOI = 0.;
 418   Double_t dSumPOI = 0.;
 419   Double_t dvproPtPOI = 0.;
 420   Double_t dErrdifcombPtPOI = 0.;
 421   Double_t dvproEtaPOI = 0.;
 422   Double_t dErrdifcombEtaPOI = 0.;
 423   //Pt:
 424   if(fHistProDiffFlowPtPOI) {
 425     for(Int_t b=1;b<=iNbinsPt;b++){
 426       dvproPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
 427       dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 428       //fill TH1D
 429       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dvproPtPOI, dErrdifcombPtPOI);
 430       if (fHistPtPOI){
 431         //integrated flow (POI)
 432         dYieldPtPOI = fHistPtPOI->GetBinContent(b);
 433         dVPOI += dvproPtPOI*dYieldPtPOI;
 434         dSumPOI += dYieldPtPOI;
 435         //error on integrated flow
 436         dErrVPOI += dYieldPtPOI*dYieldPtPOI*dErrdifcombPtPOI*dErrdifcombPtPOI;
 437       }
 438     }//end of for(Int_t b=0;b<iNbinsPt;b++)
 439   } else { cout<<"fHistProFlow is NULL"<<endl; }
 440
 441   if (!(TMath::AreEqualAbs(dSumPOI, 0.0, 1e-10)) ) {
 442     dVPOI /= dSumPOI;  //because pt distribution should be normalised
 443     dErrVPOI /= (dSumPOI*dSumPOI);
 444     dErrVPOI = TMath::Sqrt(dErrVPOI);
 445   }
 446   cout<<"dV"<<fHarmonic<<"{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;
 447
 448   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
 449
 450   //Eta:
 451   if(fHistProDiffFlowEtaPOI)
 452   {
 453    for(Int_t b=1;b<=iNbinsEta;b++)
 454    {
 455     dvproEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
 456     dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 457     //fill common hist results:
 458     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dvproEtaPOI, dErrdifcombEtaPOI);
 459    }
 460   }
 461
 462   cout<<endl;
 463   //cout<<".....finished"<<endl;
 464 }
 465
 466 void AliFlowAnalysisWithMCEventPlane::BookObjectsForFlowOfResonances()
 467 {
 468  // Book all objecsts needed to study flow of resonances.
 469
 470  // List:
 471  fResonanceList->SetName("Resonances");
 472  fResonanceList->SetOwner(kTRUE);
 473  fHistList->Add(fResonanceList);
 474
 475  // Profile holding settings for flow of resoances:
 476  TString resonanceSettingsName = "fResonanceSettings";
 477  fResonanceSettings = new TProfile(resonanceSettingsName.Data(),"Settings for flow of resonances",2,0,2);
 478  //fResonanceSettings->GetXaxis()->SetLabelSize(0.025);
 479  fResonanceSettings->GetXaxis()->SetBinLabel(1,"fFlowOfResonances");
 480  fResonanceSettings->Fill(0.5,(Int_t)fFlowOfResonances);
 481  fResonanceSettings->GetXaxis()->SetBinLabel(2,"x in #phi_{pair}"); // phi_{pair} = x*phi1+(1-x)*phi2
 482  fResonanceSettings->Fill(1.5,fXinPairAngle);
 483  fResonanceList->Add(fResonanceSettings);
 484
 485  // Profiles used to calculate <cos[n(phi_{pair}-RP)]> and <sin[n(phi_{pair}-RP)]> (0 = cos, 1 = sin), where phi_{pair} = x*phi1+(1-x)*phi2:
 486  TString cosSinFlag[2] = {"Cos","Sin"};
 487  TString cosSinTitleFlag[2] = {"cos[n(#phi_{pair}-#psi_{RP})]","sin[n(#phi_{pair}-#psi_{RP})]"};
 488  TString pairCorrelatorName = "fPairCorrelator";
 489  for(Int_t cs=0;cs<2;cs++)
 490  {
 491   fPairCorrelator[cs] = new TProfile(Form("%s%s",pairCorrelatorName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),1,0.,1.);
 492   fPairCorrelator[cs]->GetXaxis()->SetBinLabel(1,cosSinTitleFlag[cs].Data());
 493   fResonanceList->Add(fPairCorrelator[cs]);
 494  }
 495
 496 } // end of void AliFlowAnalysisWithMCEventPlane::BookObjectsForFlowOfResonances()
 497
 498 void AliFlowAnalysisWithMCEventPlane::FlowOfResonances(AliFlowEventSimple* anEvent)
 499 {
 500  // Evaluate correlators relevant for the flow of resonances.
 501
 502  // Get the MC reaction plane angle:
 503  Double_t RP = anEvent->GetMCReactionPlaneAngle();
 504  // Get the number of tracks:
 505  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 506  AliFlowTrackSimple *pTrack = NULL;
 507  Double_t dPhi1 = 0.;
 508  Double_t dPhi2 = 0.;
 509  Double_t x = fXinPairAngle; // shortcut
 510  Double_t n = fHarmonic; // shortcut
 511  for(Int_t i=0;i<iNumberOfTracks;i++)
 512  {
 513   pTrack = anEvent->GetTrack(i);
 514   if(pTrack && pTrack->InRPSelection())
 515   {
 516    dPhi1 = pTrack->Phi();
 517   }
 518   for(Int_t j=0;j<iNumberOfTracks;j++)
 519   {
 520    if(j==i) continue;
 521    pTrack = anEvent->GetTrack(j);
 522    if(pTrack && pTrack->InRPSelection())
 523    {
 524     dPhi2 = pTrack->Phi();
 525    }
 526    Double_t dPhiPair = x*dPhi1+(1.-x)*dPhi2;
 527    fPairCorrelator[0]->Fill(0.5,TMath::Cos(n*(dPhiPair-RP)),1.);
 528    fPairCorrelator[1]->Fill(0.5,TMath::Sin(n*(dPhiPair-RP)),1.);
 529   } // end of for(Int_t j=i+1;j<iNumberOfTracks;j++)
 530  } // end of for(Int_t i=0;i<iNumberOfTracks;i++)
 531
 532 } // end of void AliFlowAnalysisWithMCEventPlane::FlowOfResonances()
 533
 534