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