PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.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 // AliFlowAnalysisWithLYZEventPlane:
  17 //
  18 // Class to do flow analysis with the event plane from the LYZ method
  19 //
  20 // author: N. van der Kolk (kolk@nikhef.nl)
  21
  22 /*
  23 $Log$
  24 */
  25
  26 //#define AliFlowAnalysisWithLYZEventPlane_cxx
  27
  28 #include "Riostream.h"  //needed as include
  29 #include "TMath.h"   //needed as include
  30 #include "TProfile.h"   //needed as include
  31 #include "TH1F.h"
  32 #include "TFile.h"
  33 #include "TList.h"
  34 #include "TVector2.h"
  35
  36 #include "AliFlowLYZConstants.h"    //needed as include
  37 #include "AliFlowCommonConstants.h" //needed as include
  38 #include "AliFlowEventSimple.h"
  39 #include "AliFlowTrackSimple.h"
  40 #include "AliFlowCommonHist.h"
  41 #include "AliFlowCommonHistResults.h"
  42 #include "AliFlowLYZEventPlane.h"
  43 #include "AliFlowAnalysisWithLYZEventPlane.h"
  44 #include "AliFlowVector.h"
  45
  46
  47 ClassImp(AliFlowAnalysisWithLYZEventPlane)
  48
  49   //-----------------------------------------------------------------------
  50
  51  AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
  52    fHistList(NULL),
  53    fSecondRunList(NULL),
  54    fSecondReDtheta(NULL),
  55    fSecondImDtheta(NULL),
  56    fFirstr0theta(NULL),
  57    fHistProVetaRP(NULL),
  58    fHistProVetaPOI(NULL),
  59    fHistProVPtRP(NULL),
  60    fHistProVPtPOI(NULL),
  61    fHistProWr(NULL),
  62    fHistProWrCorr(NULL),
  63    fHistQsumforChi(NULL),
  64    fHistDeltaPhi(NULL),
  65    fHistDeltaPhi2(NULL),
  66    fHistDeltaPhihere(NULL),
  67    fHistPhiEP(NULL),
  68    fHistPhiEPhere(NULL),
  69    fHistPhiLYZ(NULL),
  70    fHistPhiLYZ2(NULL),
  71    fCommonHists(NULL),
  72    fCommonHistsRes(NULL),
  73    fEventNumber(0),
  74    fQsum(NULL),
  75    fQ2sum(0)
  76 {
  77
  78   // Constructor.
  79   fQsum = new TVector2();        // flow vector sum
  80
  81   fHistList = new TList();
  82   fSecondRunList = new TList();
  83 }
  84
  85
  86
  87  //-----------------------------------------------------------------------
  88
  89
  90  AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
  91  {
  92    //destructor
  93    delete fQsum;
  94    delete fHistList;
  95    delete fSecondRunList;
  96  }
  97
  98
  99 //-----------------------------------------------------------------------
 100
 101 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString* outputFileName)
 102 {
 103  //store the final results in output .root file
 104
 105   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 106   //output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
 107   fHistList->SetName("cobjLYZEP");
 108   fHistList->SetOwner(kTRUE);
 109   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 110   delete output;
 111 }
 112
 113 //-----------------------------------------------------------------------
 114
 115 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString outputFileName)
 116 {
 117  //store the final results in output .root file
 118
 119   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 120   //output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
 121   fHistList->SetName("cobjLYZEP");
 122   fHistList->SetOwner(kTRUE);
 123   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 124   delete output;
 125 }
 126
 127 //-----------------------------------------------------------------------
 128
 129 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TDirectoryFile *outputFileName)
 130 {
 131  //store the final results in output .root file
 132  fHistList->SetName("cobjLYZEP");
 133  fHistList->SetOwner(kTRUE);
 134  outputFileName->Add(fHistList);
 135  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 136 }
 137
 138 //-----------------------------------------------------------------------
 139
 140 void AliFlowAnalysisWithLYZEventPlane::Init() {
 141
 142   //Initialise all histograms
 143   cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
 144
 145   //save old value and prevent histograms from being added to directory
 146   //to avoid name clashes in case multiple analaysis objects are used
 147   //in an analysis
 148   Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 149   TH1::AddDirectory(kFALSE);
 150
 151   //input histograms
 152   if (fSecondRunList) {
 153     fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZSUM");
 154     fHistList->Add(fSecondReDtheta);
 155
 156     fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZSUM");
 157     fHistList->Add(fSecondImDtheta);
 158
 159     fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZSUM");
 160     fHistList->Add(fFirstr0theta);
 161
 162     //warnings
 163     if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
 164     if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
 165     if (!fFirstr0theta)   {cout<<"fFirstr0theta is NULL!"<<endl; }
 166
 167   }
 168
 169   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
 170   fHistList->Add(fCommonHists);
 171
 172   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
 173   fHistList->Add(fCommonHistsRes);
 174
 175   Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 176   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 177   Double_t  dPtMin  = AliFlowCommonConstants::GetMaster()->GetPtMin();
 178   Double_t  dPtMax  = AliFlowCommonConstants::GetMaster()->GetPtMax();
 179   Double_t  dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 180   Double_t  dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 181
 182   fHistProVetaRP = new TProfile("FlowPro_VetaRP_LYZEP","FlowPro_VetaRP_LYZEP",iNbinsEta,dEtaMin,dEtaMax);
 183   fHistProVetaRP->SetXTitle("rapidity");
 184   fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
 185   fHistList->Add(fHistProVetaRP);
 186
 187   fHistProVetaPOI = new TProfile("FlowPro_VetaPOI_LYZEP","FlowPro_VetaPOI_LYZEP",iNbinsEta,dEtaMin,dEtaMax);
 188   fHistProVetaPOI->SetXTitle("rapidity");
 189   fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
 190   fHistList->Add(fHistProVetaPOI);
 191
 192   fHistProVPtRP = new TProfile("FlowPro_VPtRP_LYZEP","FlowPro_VPtRP_LYZEP",iNbinsPt,dPtMin,dPtMax);
 193   fHistProVPtRP->SetXTitle("Pt");
 194   fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
 195   fHistList->Add(fHistProVPtRP);
 196
 197   fHistProVPtPOI = new TProfile("FlowPro_VPtPOI_LYZEP","FlowPro_VPtPOI_LYZEP",iNbinsPt,dPtMin,dPtMax);
 198   fHistProVPtPOI->SetXTitle("p_{T}");
 199   fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
 200   fHistList->Add(fHistProVPtPOI);
 201
 202   fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
 203   fHistProWr->SetXTitle("Q");
 204   fHistProWr->SetYTitle("Wr");
 205   fHistList->Add(fHistProWr);
 206
 207   fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
 208   fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
 209   fHistQsumforChi->SetYTitle("value");
 210   fHistList->Add(fHistQsumforChi);
 211
 212   fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
 213   fHistDeltaPhi->SetXTitle("DeltaPhi");
 214   fHistDeltaPhi->SetYTitle("Counts");
 215   fHistList->Add(fHistDeltaPhi);
 216
 217   fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
 218   fHistPhiLYZ->SetXTitle("Phi from LYZ");
 219   fHistPhiLYZ->SetYTitle("Counts");
 220   fHistList->Add(fHistPhiLYZ);
 221
 222   fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
 223   fHistPhiEP->SetXTitle("Phi from EP");
 224   fHistPhiEP->SetYTitle("Counts");
 225   fHistList->Add(fHistPhiEP);
 226
 227   fEventNumber = 0;  //set number of events to zero
 228
 229   //restore old status
 230   TH1::AddDirectory(oldHistAddStatus);
 231 }
 232
 233 //-----------------------------------------------------------------------
 234
 235 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
 236
 237   //Get the event plane and weight for each event
 238   if (anEvent) {
 239
 240     //fill control histograms
 241     fCommonHists->FillControlHistograms(anEvent);
 242
 243     //get the Q vector from the FlowEvent
 244     AliFlowVector vQ = anEvent->GetQ();
 245     //if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; } //coding violation
 246     //Weight with the multiplicity
 247     Double_t dQX = 0.;
 248     Double_t dQY = 0.;
 249     if (TMath::AreEqualAbs(vQ.GetMult(),0.0,1e-10)) {
 250       dQX = vQ.X()/vQ.GetMult();
 251       dQY = vQ.Y()/vQ.GetMult();
 252     } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
 253     vQ.Set(dQX,dQY);
 254     //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
 255
 256     //for chi calculation:
 257     *fQsum += vQ;
 258     fHistQsumforChi->SetBinContent(1,fQsum->X());
 259     fHistQsumforChi->SetBinContent(2,fQsum->Y());
 260     fQ2sum += vQ.Mod2();
 261     fHistQsumforChi->SetBinContent(3,fQ2sum);
 262     //cout<<"fQ2sum = "<<fQ2sum<<endl;
 263
 264     //call AliFlowLYZEventPlane::CalculateRPandW() here!
 265     aLYZEP->CalculateRPandW(vQ);
 266
 267     Double_t dWR = aLYZEP->GetWR();
 268     Double_t dRP = aLYZEP->GetPsi();
 269
 270     //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
 271     fHistPhiLYZ->Fill(dRP);
 272
 273     //plot difference between event plane from EP-method and LYZ-method
 274     Double_t dRPEP = vQ.Phi()/2;                              //gives distribution from (0 to pi)
 275     //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X());       //gives distribution from (-pi/2 to pi/2)
 276     //cout<<"dRPEP = "<< dRPEP <<endl;
 277     fHistPhiEP->Fill(dRPEP);
 278
 279     Double_t dDeltaPhi = dRPEP - dRP;
 280     if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); }        //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
 281     //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
 282     fHistDeltaPhi->Fill(dDeltaPhi);
 283
 284     //Flip sign of WR
 285     Double_t dLow = TMath::Pi()/4.;
 286     Double_t dHigh = 3.*(TMath::Pi()/4.);
 287     if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
 288       dRP -= (TMath::Pi()/2);
 289       dWR = -dWR;
 290       cerr<<"*** dRP modified ***"<<endl;
 291     }
 292     fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
 293
 294     //calculate flow for RP and POI selections
 295     //loop over the tracks of the event
 296     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 297     for (Int_t i=0;i<iNumberOfTracks;i++)
 298       {
 299         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 300         if (pTrack){
 301           Double_t dPhi = pTrack->Phi();
 302           //if (dPhi<0.) fPhi+=2*TMath::Pi();
 303           Double_t dPt  = pTrack->Pt();
 304           Double_t dEta = pTrack->Eta();
 305           //calculate flow v2:
 306           Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
 307           if (pTrack->InRPSelection()) {
 308             //fill histograms for RP selection
 309             fHistProVetaRP -> Fill(dEta,dv2);
 310             fHistProVPtRP  -> Fill(dPt,dv2);
 311           }
 312           if (pTrack->InPOISelection()) {
 313             //fill histograms for POI selection
 314             fHistProVetaPOI -> Fill(dEta,dv2);
 315             fHistProVPtPOI  -> Fill(dPt,dv2);
 316           }
 317         }//track
 318       }//loop over tracks
 319
 320     fEventNumber++;
 321     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 322   }
 323 }
 324
 325   //--------------------------------------------------------------------
 326 void AliFlowAnalysisWithLYZEventPlane::GetOutputHistograms(TList *outputListHistos){
 327  //get pointers to all output histograms (called before Finish())
 328  if (outputListHistos) {
 329     //Get the common histograms from the output list
 330     AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
 331       (outputListHistos->FindObject("AliFlowCommonHistLYZEP"));
 332     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 333       (outputListHistos->FindObject("AliFlowCommonHistResultsLYZEP"));
 334
 335     TProfile* pHistProR0theta = dynamic_cast<TProfile*>
 336       (outputListHistos->FindObject("First_FlowPro_r0theta_LYZSUM"));
 337
 338     TProfile* pHistProVetaRP = dynamic_cast<TProfile*>
 339       (outputListHistos->FindObject("FlowPro_VetaRP_LYZEP"));
 340     TProfile* pHistProVetaPOI = dynamic_cast<TProfile*>
 341       (outputListHistos->FindObject("FlowPro_VetaPOI_LYZEP"));
 342     TProfile* pHistProVPtRP = dynamic_cast<TProfile*>
 343       (outputListHistos->FindObject("FlowPro_VPtRP_LYZEP"));
 344     TProfile* pHistProVPtPOI = dynamic_cast<TProfile*>
 345       (outputListHistos->FindObject("FlowPro_VPtPOI_LYZEP"));
 346
 347     TH1F* pHistQsumforChi = dynamic_cast<TH1F*>
 348       (outputListHistos->FindObject("Flow_QsumforChi_LYZEP"));
 349
 350     if (pCommonHist && pCommonHistResults && pHistProR0theta &&
 351         pHistProVetaRP && pHistProVetaPOI && pHistProVPtRP &&
 352         pHistProVPtPOI && pHistQsumforChi ) {
 353       this -> SetCommonHists(pCommonHist);
 354       this -> SetCommonHistsRes(pCommonHistResults);
 355       this -> SetFirstr0theta(pHistProR0theta);
 356       this -> SetHistProVetaRP(pHistProVetaRP);
 357       this -> SetHistProVetaPOI(pHistProVetaPOI);
 358       this -> SetHistProVPtRP(pHistProVPtRP);
 359       this -> SetHistProVPtPOI(pHistProVPtPOI);
 360       this -> SetHistQsumforChi(pHistQsumforChi);
 361      }
 362   } else {
 363       cout<<"WARNING: Histograms needed to run Finish() are not accessable!"<<endl;
 364     }
 365 }
 366
 367   //--------------------------------------------------------------------
 368 void AliFlowAnalysisWithLYZEventPlane::Finish() {
 369
 370   //plot histograms etc.
 371   cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
 372
 373   //constants:
 374   Double_t  dJ01 = 2.405;
 375   Int_t iNtheta   = AliFlowLYZConstants::GetMaster()->GetNtheta();
 376   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 377   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 378   //set the event number
 379   if (fCommonHists) {
 380   SetEventNumber((int)fCommonHists->GetHistQ()->GetEntries());
 381   //cout<<"number of events processed is "<<fEventNumber<<endl;
 382   }
 383
 384   //set the sum of Q vectors
 385   fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
 386   SetQ2sum(fHistQsumforChi->GetBinContent(3));
 387
 388   //calculate dV the mean of dVtheta
 389   Double_t  dVtheta = 0;
 390   Double_t  dV = 0;
 391   for (Int_t theta=0;theta<iNtheta;theta++)     {
 392     Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
 393     if (TMath::AreEqualAbs(dR0,0.0,1e-10)) { dVtheta = dJ01/dR0 ;}
 394     dV += dVtheta;
 395   }
 396   dV /= iNtheta;
 397
 398   //calculate dChi
 399   Double_t  dSigma2 = 0;
 400   Double_t  dChi= 0;
 401   if (fEventNumber!=0) {
 402     *fQsum /= fEventNumber;
 403     //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
 404     //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
 405     fQ2sum /= fEventNumber;
 406     //cerr<<"fEventNumber = "<<fEventNumber<<endl;
 407     //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 408     dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 409     //cerr<<"dSigma2"<<dSigma2<<endl;
 410     if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 411     else dChi = -1.;
 412     fCommonHistsRes->FillChiRP(dChi);
 413
 414     // recalculate statistical errors on integrated flow
 415     //combining 5 theta angles to 1 relative error BP eq. 89
 416     Double_t dRelErr2comb = 0.;
 417     Int_t iEvts = fEventNumber;
 418     if (iEvts!=0) {
 419       for (Int_t theta=0;theta<iNtheta;theta++){
 420         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 421         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 422                                        TMath::Cos(dTheta));
 423         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 424                                       TMath::Cos(dTheta));
 425         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 426                           TMath::BesselJ1(dJ01)))*
 427           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 428            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 429       }
 430       dRelErr2comb /= iNtheta;
 431     }
 432     Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 433     Double_t dVErr = dV*dRelErrcomb ;
 434     fCommonHistsRes->FillIntegratedFlow(dV, dVErr);
 435
 436     cout<<"*************************************"<<endl;
 437     cout<<"*************************************"<<endl;
 438     cout<<"      Integrated flow from           "<<endl;
 439     cout<<"  Lee-Yang Zeroes Event Plane        "<<endl;
 440     cout<<endl;
 441     cout<<"dChi = "<<dChi<<endl;
 442     cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
 443     cout<<endl;
 444
 445   }
 446
 447   //copy content of profile into TH1D, add error and fill the AliFlowCommonHistResults
 448
 449   //v as a function of eta for RP selection
 450   for(Int_t b=0;b<iNbinsEta;b++) {
 451     Double_t dv2pro  = fHistProVetaRP->GetBinContent(b);
 452     Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
 453     Double_t dErrdifcomb = 0.;  //set error to zero
 454     Double_t dErr2difcomb = 0.; //set error to zero
 455     //calculate error
 456     if (TMath::AreEqualAbs(dNprime,0.0,1e-10)) {
 457       for (Int_t theta=0;theta<iNtheta;theta++) {
 458         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 459         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 460                                            TMath::Cos(dTheta));
 461         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 462                                           TMath::Cos(dTheta));
 463         dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 464                                                  TMath::BesselJ1(dJ01)))*
 465           ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 466            (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 467       } //loop over theta
 468     }
 469
 470     if (TMath::AreEqualAbs(dErr2difcomb, 0.0, 1e-10)) {
 471       dErr2difcomb /= iNtheta;
 472       dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 473     }
 474     else {dErrdifcomb = 0.;}
 475     //fill TH1D
 476     fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
 477   } //loop over bins b
 478
 479
 480   //v as a function of eta for POI selection
 481   for(Int_t b=0;b<iNbinsEta;b++) {
 482     Double_t dv2pro  = fHistProVetaPOI->GetBinContent(b);
 483     Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
 484     Double_t dErrdifcomb = 0.;  //set error to zero
 485     Double_t dErr2difcomb = 0.; //set error to zero
 486     //calculate error
 487     if (dNprime!=0.) {
 488       for (Int_t theta=0;theta<iNtheta;theta++) {
 489         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 490         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 491                                          TMath::Cos(dTheta));
 492         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 493                                         TMath::Cos(dTheta));
 494         dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 495                                              TMath::BesselJ1(dJ01)))*
 496           ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 497            (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 498       } //loop over theta
 499     }
 500
 501     if (dErr2difcomb!=0.) {
 502       dErr2difcomb /= iNtheta;
 503       dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 504     }
 505     else {dErrdifcomb = 0.;}
 506     //fill TH1D
 507     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
 508   } //loop over bins b
 509
 510   //v as a function of Pt for RP selection
 511   TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 512   Double_t dVRP = 0.;
 513   Double_t dSum = 0.;
 514   Double_t dErrV =0.;
 515
 516   for(Int_t b=0;b<iNbinsPt;b++) {
 517     Double_t dv2pro  = fHistProVPtRP->GetBinContent(b);
 518     Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
 519     Double_t dErrdifcomb = 0.;  //set error to zero
 520     Double_t dErr2difcomb = 0.; //set error to zero
 521     //calculate error
 522     if (dNprime!=0.) {
 523       for (Int_t theta=0;theta<iNtheta;theta++) {
 524         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 525         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 526                                            TMath::Cos(dTheta));
 527         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 528                                           TMath::Cos(dTheta));
 529         dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 530                                                TMath::BesselJ1(dJ01)))*
 531           ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 532            (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 533       } //loop over theta
 534     }
 535
 536     if (dErr2difcomb!=0.) {
 537       dErr2difcomb /= iNtheta;
 538       dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 539       //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 540     }
 541     else {dErrdifcomb = 0.;}
 542
 543     //fill TH1D
 544     fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
 545     //calculate integrated flow for RP selection
 546     if (fHistPtRP){
 547       Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 548       dVRP += dv2pro*dYieldPt;
 549       dSum +=dYieldPt;
 550       dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
 551     } else { cout<<"fHistPtRP is NULL"<<endl; }
 552
 553   } //loop over bins b
 554
 555   if (TMath::AreEqualAbs(dSum, 0.0, 1e-10)) {
 556     dVRP /= dSum; //the pt distribution should be normalised
 557     dErrV /= (dSum*dSum);
 558     dErrV = TMath::Sqrt(dErrV);
 559   }
 560   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
 561
 562   cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
 563   //cout<<endl;
 564
 565
 566   //v as a function of Pt for POI selection
 567   TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 568   Double_t dVPOI = 0.;
 569   dSum = 0.;
 570   dErrV =0.;
 571
 572   for(Int_t b=0;b<iNbinsPt;b++) {
 573     Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
 574     Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
 575
 576     //cerr<<"dNprime = "<<dNprime<<endl;
 577     //Int_t iNprime = TMath::Nint(fHistProVPtPOI->GetBinEntries(b));
 578     //cerr<<"iNprime = "<<iNprime<<endl;
 579
 580     Double_t dErrdifcomb = 0.;  //set error to zero
 581     Double_t dErr2difcomb = 0.; //set error to zero
 582     //calculate error
 583     if (dNprime!=0.) {
 584       for (Int_t theta=0;theta<iNtheta;theta++) {
 585         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 586         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 587                                            TMath::Cos(dTheta));
 588         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 589                                           TMath::Cos(dTheta));
 590         dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 591                                                  TMath::BesselJ1(dJ01)))*
 592           ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 593            (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 594       } //loop over theta
 595     }
 596
 597     if (dErr2difcomb!=0.) {
 598       dErr2difcomb /= iNtheta;
 599       dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 600       //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 601     }
 602     else {dErrdifcomb = 0.;}
 603
 604     //fill TH1D
 605     fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
 606
 607     //calculate integrated flow for POI selection
 608     if (fHistPtPOI){
 609       Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 610       dVPOI += dv2pro*dYieldPt;
 611       dSum +=dYieldPt;
 612       dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
 613     } else { cout<<"fHistPtPOI is NULL"<<endl; }
 614   } //loop over bins b
 615
 616   if (dSum != 0.) {
 617     dVPOI /= dSum; //the pt distribution should be normalised
 618     dErrV /= (dSum*dSum);
 619     dErrV = TMath::Sqrt(dErrV);
 620   }
 621   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
 622
 623   cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
 624   cout<<endl;
 625   cout<<"*************************************"<<endl;
 626   cout<<"*************************************"<<endl;
 627
 628   //cout<<".....finished"<<endl;
 629  }
 630