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