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