PWG2/FLOW/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 "TComplex.h"   //needed as include
  22 #include "TList.h"
  23
  24 class TH1F;
  25
  26 #include "AliFlowCommonConstants.h"    //needed as include
  27 #include "AliFlowEventSimple.h"
  28 #include "AliFlowTrackSimple.h"
  29 #include "AliFlowCommonHist.h"
  30 #include "AliFlowCommonHistResults.h"
  31 #include "AliFlowAnalysisWithMCEventPlane.h"
  32
  33 class AliFlowVector;
  34
  35 // AliFlowAnalysisWithMCEventPlane:
  36 // Description: Maker to analyze Flow from the generated MC reaction plane.
  37 //              This class is used to get the real value of the flow
  38 //              to compare the other methods to when analysing simulated events
  39 // author: N. van der Kolk (kolk@nikhef.nl)
  40
  41 ClassImp(AliFlowAnalysisWithMCEventPlane)
  42
  43   //-----------------------------------------------------------------------
  44
  45  AliFlowAnalysisWithMCEventPlane::AliFlowAnalysisWithMCEventPlane():
  46    fQsum(NULL),
  47    fQ2sum(0),
  48    fEventNumber(0),
  49    fDebug(kFALSE),
  50    fHistList(NULL),
  51    fCommonHists(NULL),
  52    fCommonHistsRes(NULL),
  53    fHistProFlow(NULL),
  54    fHistRP(NULL)
  55
  56 {
  57
  58   // Constructor.
  59   fHistList = new TList();
  60
  61   fQsum = new TVector2;        // flow vector sum
  62 }
  63
  64
  65  //-----------------------------------------------------------------------
  66
  67
  68  AliFlowAnalysisWithMCEventPlane::~AliFlowAnalysisWithMCEventPlane()
  69  {
  70    //destructor
  71    delete fHistList;
  72    delete fQsum;
  73  }
  74
  75
  76 //-----------------------------------------------------------------------
  77 void AliFlowAnalysisWithMCEventPlane::Init() {
  78
  79   //Define all histograms
  80   cout<<"---Analysis with the real MC Event Plane---"<<endl;
  81
  82   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
  83   Double_t  dPtMin = AliFlowCommonConstants::GetPtMin();
  84   Double_t  dPtMax = AliFlowCommonConstants::GetPtMax();
  85
  86   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistMCEP");
  87   fHistList->Add(fCommonHists);
  88   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsMCEP");
  89   fHistList->Add(fCommonHistsRes);
  90
  91   fHistProFlow = new TProfile("FlowPro_VPt_MCEP","FlowPro_VPt_MCEP",iNbinsPt,dPtMin,dPtMax);
  92   fHistProFlow->SetXTitle("Pt");
  93   fHistProFlow->SetYTitle("v2 (%)");
  94   fHistList->Add(fHistProFlow);
  95
  96   fHistRP = new TH1F("Flow_RP_MCEP","Flow_RP_MCEP",100,0.,3.14);
  97   fHistRP->SetXTitle("Reaction Plane Angle");
  98   fHistRP->SetYTitle("Counts");
  99   fHistList->Add(fHistRP);
 100
 101
 102   fEventNumber = 0;  //set number of events to zero
 103
 104 }
 105
 106 //-----------------------------------------------------------------------
 107
 108 void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent, Double_t aRP) {
 109
 110   //Calculate v2 from the MC reaction plane
 111   if (anEvent) {
 112
 113     //fill control histograms
 114     fCommonHists->FillControlHistograms(anEvent);
 115
 116     //get the Q vector from the FlowEvent
 117     AliFlowVector vQ = anEvent->GetQ();
 118     //cout<<"vQ.Mod() = " << vQ.Mod() << endl;
 119     //for chi calculation:
 120     *fQsum += vQ;
 121     //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
 122     fQ2sum += vQ.Mod2();
 123     //cout<<"fQ2sum = "<<fQ2sum<<endl;
 124
 125     fHistRP->Fill(aRP);
 126
 127     //calculate flow
 128     //loop over the tracks of the event
 129     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 130     for (Int_t i=0;i<iNumberOfTracks;i++)
 131       {
 132         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 133         if (pTrack){
 134           if (pTrack->UseForDifferentialFlow()) {
 135             Double_t dPhi = pTrack->Phi();
 136             //if (dPhi<0.) dPhi+=2*TMath::Pi();
 137             //calculate flow v2:
 138             Double_t dv2 = TMath::Cos(2*(dPhi-aRP));
 139             Double_t dPt = pTrack->Pt();
 140             //fill histogram
 141             fHistProFlow->Fill(dPt,100*dv2);
 142           }
 143         }//track selected
 144       }//loop over tracks
 145
 146     fEventNumber++;
 147     cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 148   }
 149 }
 150
 151   //--------------------------------------------------------------------
 152 void AliFlowAnalysisWithMCEventPlane::Finish() {
 153
 154   //*************make histograms etc.
 155   if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
 156
 157   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 158
 159
 160   TH1F* fHistPtDiff = fCommonHists->GetHistPtDiff();
 161   Double_t dV = 0.;
 162   Double_t dErrV = 0.;
 163   Double_t dSum = 0.;
 164   for(Int_t b=0;b<iNbinsPt;b++){
 165     Double_t dv2pro = 0.;
 166     Double_t dErrdifcomb = 0.;
 167     if(fHistProFlow) {
 168       dv2pro = fHistProFlow->GetBinContent(b);
 169       dErrdifcomb = fHistProFlow->GetBinError(b); //in case error from profile is correct
 170       //fill TH1D
 171       fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
 172       if (fHistPtDiff){
 173         //integrated flow
 174         Double_t dYield = fHistPtDiff->GetBinContent(b);
 175         dV += dv2pro/100*dYield ;
 176         dSum += dYield;
 177         //error on integrated flow
 178         dErrV += dYield*dYield*(dErrdifcomb/100)*(dErrdifcomb/100);
 179       }
 180     } else { cout<<"fHistProFlow is NULL"<<endl; }
 181   }
 182   if (dSum != 0. ) {
 183     dV /= dSum;  //because pt distribution should be normalised
 184     dErrV /= dSum*dSum;
 185     dErrV = TMath::Sqrt(dErrV); }
 186   cout<<"dV is "<<dV<<" +- "<<dErrV<<endl;
 187   fCommonHistsRes->FillIntegratedFlow(dV,dErrV);
 188
 189   cout<<".....finished"<<endl;
 190 }
 191
 192
 193