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