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
  16 #define AliFlowAnalysisWithMCEventPlane_cxx
  17
  18 #include "Riostream.h"  //needed as include
  19 #include "TFile.h"      //needed as include
  20 #include "TProfile.h"   //needed as include
  21 #include "TProfile2D.h"
  22 #include "TComplex.h"   //needed as include
  23 #include "TList.h"
  24
  25 class TH1F;
  26
  27 #include "AliFlowCommonConstants.h"    //needed as include
  28 #include "AliFlowEventSimple.h"
  29 #include "AliFlowTrackSimple.h"
  30 #include "AliFlowCommonHist.h"
  31 #include "AliFlowCommonHistResults.h"
  32 #include "AliFlowAnalysisWithMCEventPlane.h"
  33
  34 class AliFlowVector;
  35
  36 // AliFlowAnalysisWithMCEventPlane:
  37 // Description: Maker to analyze Flow from the generated MC reaction plane.
  38 //              This class is used to get the real value of the flow
  39 //              to compare the other methods to when analysing simulated events
  40 // author: N. van der Kolk (kolk@nikhef.nl)
  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)\r
  64 {
  65
  66   // Constructor.
  67   fHistList = new TList();
  68
  69   fQsum = new TVector2;        // flow vector sum
  70 }
  71
  72
  73  //-----------------------------------------------------------------------
  74
  75
  76  AliFlowAnalysisWithMCEventPlane::~AliFlowAnalysisWithMCEventPlane()
  77  {
  78    //destructor
  79    delete fHistList;
  80    delete fQsum;
  81  }
  82
  83 //-----------------------------------------------------------------------
  84
  85 void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString* outputFileName)
  86 {
  87  //store the final results in output .root file
  88  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  89  //output->WriteObject(fHistList, "cobjMCEP","SingleKey");
  90  fHistList->SetName("cobjMCEP");
  91  fHistList->SetOwner(kTRUE);
  92  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  93  delete output;
  94 }
  95
  96 //-----------------------------------------------------------------------
  97
  98 void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString outputFileName)
  99 {
 100  //store the final results in output .root file
 101  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 102  //output->WriteObject(fHistList, "cobjMCEP","SingleKey");
 103  fHistList->SetName("cobjMCEP");
 104  fHistList->SetOwner(kTRUE);
 105  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 106  delete output;
 107 }
 108
 109 //-----------------------------------------------------------------------
 110 void AliFlowAnalysisWithMCEventPlane::Init() {
 111
 112   //Define all histograms
 113   cout<<"---Analysis with the real MC Event Plane---"<<endl;
 114
 115   Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 116   Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 117   Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 118
 119   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 120   Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 121   Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 122
 123   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistMCEP");
 124   fHistList->Add(fCommonHists);
 125   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsMCEP");
 126   fHistList->Add(fCommonHistsRes);
 127
 128   // store harmonic in common control histogram:
 129   (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
 130
 131   fHistRP = new TH1F("Flow_RP_MCEP","Flow_RP_MCEP",100,0.,3.14);
 132   fHistRP->SetXTitle("Reaction Plane Angle");
 133   fHistRP->SetYTitle("Counts");
 134   fHistList->Add(fHistRP);
 135
 136   fHistProIntFlow = new TProfile("FlowPro_V_MCEP","FlowPro_V_MCEP",1,0.,1.);
 137   fHistProIntFlow->SetLabelSize(0.06);
 138   (fHistProIntFlow->GetXaxis())->SetBinLabel(1,"v_{n}{2}");
 139   fHistProIntFlow->SetYTitle("");
 140   fHistList->Add(fHistProIntFlow);
 141
 142   fHistProDiffFlowPtEtaRP = new TProfile2D("FlowPro_VPtEtaRP_MCEP","FlowPro_VPtEtaRP_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
 143   fHistProDiffFlowPtEtaRP->SetXTitle("P_{t}");
 144   fHistProDiffFlowPtEtaRP->SetYTitle("#eta");
 145   fHistList->Add(fHistProDiffFlowPtEtaRP);
 146
 147   fHistProDiffFlowPtRP = new TProfile("FlowPro_VPtRP_MCEP","FlowPro_VPtRP_MCEP",iNbinsPt,dPtMin,dPtMax);
 148   fHistProDiffFlowPtRP->SetXTitle("P_{t}");
 149   fHistProDiffFlowPtRP->SetYTitle("");
 150   fHistList->Add(fHistProDiffFlowPtRP);
 151
 152   fHistProDiffFlowEtaRP = new TProfile("FlowPro_VetaRP_MCEP","FlowPro_VetaRP_MCEP",iNbinsEta,dEtaMin,dEtaMax);
 153   fHistProDiffFlowEtaRP->SetXTitle("#eta");
 154   fHistProDiffFlowEtaRP->SetYTitle("");
 155   fHistList->Add(fHistProDiffFlowEtaRP);
 156
 157   fHistProDiffFlowPtEtaPOI = new TProfile2D("FlowPro_VPtEtaPOI_MCEP","FlowPro_VPtEtaPOI_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
 158   fHistProDiffFlowPtEtaPOI->SetXTitle("P_{t}");
 159   fHistProDiffFlowPtEtaPOI->SetYTitle("#eta");
 160   fHistList->Add(fHistProDiffFlowPtEtaPOI);
 161
 162   fHistProDiffFlowPtPOI = new TProfile("FlowPro_VPtPOI_MCEP","FlowPro_VPtPOI_MCEP",iNbinsPt,dPtMin,dPtMax);
 163   fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
 164   fHistProDiffFlowPtPOI->SetYTitle("");
 165   fHistList->Add(fHistProDiffFlowPtPOI);
 166
 167   fHistProDiffFlowEtaPOI = new TProfile("FlowPro_VetaPOI_MCEP","FlowPro_VetaPOI_MCEP",iNbinsEta,dEtaMin,dEtaMax);
 168   fHistProDiffFlowEtaPOI->SetXTitle("#eta");
 169   fHistProDiffFlowEtaPOI->SetYTitle("");
 170   fHistList->Add(fHistProDiffFlowEtaPOI);
 171
 172   fHistSpreadOfFlow = new TH1D("fHistSpreadOfFlow","fHistSpreadOfFlow",1000,-1,1);
 173   fHistSpreadOfFlow->SetXTitle("v_{2}");
 174   fHistSpreadOfFlow->SetYTitle("counts");
 175   fHistList->Add(fHistSpreadOfFlow);
 176
 177   fEventNumber = 0;  //set number of events to zero
 178
 179 }
 180
 181 //-----------------------------------------------------------------------
 182
 183 void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent) {
 184
 185   //Calculate v2 from the MC reaction plane
 186   if (anEvent) {
 187
 188     // get the MC reaction plane angle
 189     Double_t aRP = anEvent->GetMCReactionPlaneAngle();
 190     //fill control histograms
 191     fCommonHists->FillControlHistograms(anEvent);
 192
 193     //get the Q vector from the FlowEvent
 194     AliFlowVector vQ = anEvent->GetQ(fHarmonic);
 195     //cout<<"vQ.Mod() = " << vQ.Mod() << endl;
 196     //for chi calculation:
 197     *fQsum += vQ;
 198     //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
 199     fQ2sum += vQ.Mod2();
 200     //cout<<"fQ2sum = "<<fQ2sum<<endl;
 201
 202     fHistRP->Fill(aRP);
 203
 204     Double_t dPhi = 0.;
 205     Double_t dv  = 0.;
 206     Double_t dPt  = 0.;
 207     Double_t dEta = 0.;
 208     //Double_t dPi = TMath::Pi();
 209
 210     // profile to calculate flow e-b-y:
 211     TProfile *flowEBE = new TProfile("flowEBE","flowEBE",1,0,1);
 212
 213     //calculate flow
 214     //loop over the tracks of the event
 215     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 216     for (Int_t i=0;i<iNumberOfTracks;i++)
 217       {
 218         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 219         if (pTrack){
 220           if (pTrack->InRPSelection()){
 221             dPhi = pTrack->Phi();
 222             dv  = TMath::Cos(fHarmonic*(dPhi-aRP));
 223                  dPt  = pTrack->Pt();
 224                  dEta = pTrack->Eta();
 225             //no-name int. flow (to be improved = name needed!):
 226             fHistProIntFlow->Fill(0.,dv);
 227             //no-name int. flow e-b-e (to be improved = name needed!):
 228             flowEBE->Fill(0.,dv);
 229             //differential flow (Pt, Eta, RP):
 230             fHistProDiffFlowPtEtaRP->Fill(dPt,dEta,dv,1.);
 231             //differential flow (Pt, RP):
 232             fHistProDiffFlowPtRP->Fill(dPt,dv,1.);
 233             //differential flow (Eta, RP):
 234             fHistProDiffFlowEtaRP->Fill(dEta,dv,1.);
 235           }
 236           if (pTrack->InPOISelection()) {
 237             dPhi = pTrack->Phi();
 238             //if (dPhi<0.) dPhi+=2*TMath::Pi();
 239             //calculate flow v2:
 240             dv  = TMath::Cos(fHarmonic*(dPhi-aRP));
 241             dPt  = pTrack->Pt();
 242             dEta = pTrack->Eta();
 243             //differential flow (Pt, Eta, POI):
 244             fHistProDiffFlowPtEtaPOI->Fill(dPt,dEta,dv,1.);
 245             //differential flow (Pt, POI):
 246             fHistProDiffFlowPtPOI->Fill(dPt,dv,1.);
 247             //differential flow (Eta, POI):
 248             fHistProDiffFlowEtaPOI->Fill(dEta,dv,1.);
 249           }
 250         }//track selected
 251       }//loop over tracks
 252
 253     fEventNumber++;
 254     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 255
 256     // store flow value for this event:
 257     fHistSpreadOfFlow->Fill(flowEBE->GetBinContent(1),flowEBE->GetBinEntries(1));
 258     delete flowEBE;
 259   }
 260 }
 261   //--------------------------------------------------------------------
 262
 263 void AliFlowAnalysisWithMCEventPlane::GetOutputHistograms(TList *outputListHistos) {
 264  // get the pointers to all output histograms before calling Finish()
 265  if (outputListHistos) {
 266     //Get the common histograms from the output list
 267     AliFlowCommonHist *pCommonHists = dynamic_cast<AliFlowCommonHist*>
 268       (outputListHistos->FindObject("AliFlowCommonHistMCEP"));
 269     AliFlowCommonHistResults *pCommonHistResults =
 270       dynamic_cast<AliFlowCommonHistResults*>
 271       (outputListHistos->FindObject("AliFlowCommonHistResultsMCEP"));
 272
 273     TProfile *pHistProIntFlow = dynamic_cast<TProfile*>
 274       (outputListHistos->FindObject("FlowPro_V_MCEP"));
 275
 276     TProfile2D *pHistProDiffFlowPtEtaRP = dynamic_cast<TProfile2D*>
 277       (outputListHistos->FindObject("FlowPro_VPtEtaRP_MCEP"));
 278
 279     TProfile *pHistProDiffFlowPtRP = dynamic_cast<TProfile*>
 280       (outputListHistos->FindObject("FlowPro_VPtRP_MCEP"));
 281
 282     TProfile *pHistProDiffFlowEtaRP = dynamic_cast<TProfile*>
 283       (outputListHistos->FindObject("FlowPro_VetaRP_MCEP"));
 284
 285     TProfile2D *pHistProDiffFlowPtEtaPOI = dynamic_cast<TProfile2D*>
 286       (outputListHistos->FindObject("FlowPro_VPtEtaPOI_MCEP"));
 287
 288     TProfile *pHistProDiffFlowPtPOI = dynamic_cast<TProfile*>
 289       (outputListHistos->FindObject("FlowPro_VPtPOI_MCEP"));
 290
 291     TProfile *pHistProDiffFlowEtaPOI = dynamic_cast<TProfile*>
 292       (outputListHistos->FindObject("FlowPro_VetaPOI_MCEP"));
 293
 294     if (pCommonHists && pCommonHistResults && pHistProIntFlow &&
 295         pHistProDiffFlowPtRP && pHistProDiffFlowEtaRP &&
 296         pHistProDiffFlowPtPOI && pHistProDiffFlowEtaPOI) {
 297       this->SetCommonHists(pCommonHists);
 298       this->SetCommonHistsRes(pCommonHistResults);
 299       this->SetHistProIntFlow(pHistProIntFlow);
 300       this->SetHistProDiffFlowPtEtaRP(pHistProDiffFlowPtEtaRP);
 301       this->SetHistProDiffFlowPtRP(pHistProDiffFlowPtRP);
 302       this->SetHistProDiffFlowEtaRP(pHistProDiffFlowEtaRP);
 303       this->SetHistProDiffFlowPtEtaPOI(pHistProDiffFlowPtEtaPOI);
 304       this->SetHistProDiffFlowPtPOI(pHistProDiffFlowPtPOI);
 305       this->SetHistProDiffFlowEtaPOI(pHistProDiffFlowEtaPOI);
 306     } else {
 307       cout<<"WARNING: Histograms needed to run Finish() are not accessible!"<<endl;  }
 308
 309     //fListHistos->Print();
 310   } else { cout << "histogram list pointer is empty" << endl;}
 311
 312 }
 313
 314   //--------------------------------------------------------------------
 315 void AliFlowAnalysisWithMCEventPlane::Finish() {
 316
 317   //*************make histograms etc.
 318   if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
 319
 320   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 321   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 322
 323   // access harmonic:
 324   if(fCommonHists && fCommonHists->GetHarmonic())\r
 325   {
 326    fHarmonic = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1); // to be improved (moved somewhere else?)
 327   }
 328
 329   // no-name int. flow (to be improved):
 330   Double_t dV = fHistProIntFlow->GetBinContent(1);
 331   Double_t dErrV = fHistProIntFlow->GetBinError(1); // to be improved (treatment of errors for non-Gaussian distribution needed!)
 332   // fill no-name int. flow (to be improved):
 333   fCommonHistsRes->FillIntegratedFlow(dV,dErrV);
 334   cout<<"dV"<<fHarmonic<<"{MC} is       "<<dV<<" +- "<<dErrV<<endl;
 335
 336   //RP:
 337   TH1F* fHistPtRP = fCommonHists->GetHistPtRP();
 338   Double_t dYieldPtRP = 0.;
 339   Double_t dVRP = 0.;
 340   Double_t dErrVRP = 0.;
 341   Double_t dSumRP = 0.;
 342   //differential flow (RP, Pt):
 343   Double_t dvPtRP = 0.;
 344   Double_t dErrvPtRP = 0.;
 345   for(Int_t b=1;b<=iNbinsPt;b++)
 346   {
 347    dvPtRP    = fHistProDiffFlowPtRP->GetBinContent(b);
 348    dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 349    fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
 350    if(fHistPtRP){
 351         //integrated flow (RP)
 352         dYieldPtRP = fHistPtRP->GetBinContent(b);
 353         dVRP += dvPtRP*dYieldPtRP;
 354         dSumRP += dYieldPtRP;
 355         //error on integrated flow
 356         dErrVRP += dYieldPtRP*dYieldPtRP*dErrvPtRP*dErrvPtRP;
 357       }
 358   }
 359   if (dSumRP != 0. ) {
 360     dVRP /= dSumRP;  //because pt distribution should be normalised
 361     dErrVRP /= (dSumRP*dSumRP);
 362     dErrVRP = TMath::Sqrt(dErrVRP);
 363   }
 364   // fill integrated flow (RP):
 365   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
 366   cout<<"dV"<<fHarmonic<<"{MC} (RP) is  "<<dVRP<<" +- "<<dErrVRP<<endl;
 367
 368   //differential flow (RP, Eta):
 369   Double_t dvEtaRP = 0.;
 370   Double_t dErrvEtaRP = 0.;
 371   for(Int_t b=1;b<=iNbinsEta;b++)
 372   {
 373    dvEtaRP    = fHistProDiffFlowEtaRP->GetBinContent(b);
 374    dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 375    fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
 376   }
 377
 378   //POI:
 379   TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI();
 380   Double_t dYieldPtPOI = 0.;
 381   Double_t dVPOI = 0.;
 382   Double_t dErrVPOI = 0.;
 383   Double_t dSumPOI = 0.;
 384   Double_t dvproPtPOI = 0.;
 385   Double_t dErrdifcombPtPOI = 0.;
 386   Double_t dvproEtaPOI = 0.;
 387   Double_t dErrdifcombEtaPOI = 0.;
 388   //Pt:
 389   if(fHistProDiffFlowPtPOI) {
 390     for(Int_t b=1;b<=iNbinsPt;b++){
 391       dvproPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
 392       dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 393       //fill TH1D
 394       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dvproPtPOI, dErrdifcombPtPOI);
 395       if (fHistPtPOI){
 396         //integrated flow (POI)
 397         dYieldPtPOI = fHistPtPOI->GetBinContent(b);
 398         dVPOI += dvproPtPOI*dYieldPtPOI;
 399         dSumPOI += dYieldPtPOI;
 400         //error on integrated flow
 401         dErrVPOI += dYieldPtPOI*dYieldPtPOI*dErrdifcombPtPOI*dErrdifcombPtPOI;
 402       }
 403     }//end of for(Int_t b=0;b<iNbinsPt;b++)
 404   } else { cout<<"fHistProFlow is NULL"<<endl; }
 405   if (dSumPOI != 0. ) {
 406     dVPOI /= dSumPOI;  //because pt distribution should be normalised
 407     dErrVPOI /= (dSumPOI*dSumPOI);
 408     dErrVPOI = TMath::Sqrt(dErrVPOI);
 409   }
 410   cout<<"dV"<<fHarmonic<<"{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;
 411
 412   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
 413
 414   //Eta:
 415   if(fHistProDiffFlowEtaPOI)
 416   {
 417    for(Int_t b=1;b<=iNbinsEta;b++)
 418    {
 419     dvproEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
 420     dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
 421     //fill common hist results:
 422     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dvproEtaPOI, dErrdifcombEtaPOI);
 423    }
 424   }
 425
 426   cout<<endl;
 427   //cout<<".....finished"<<endl;
 428 }
 429
 430
 431