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