fix coding violations
[u/mrichter/AliRoot.git] / 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