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