POI's and RP's for LeeYang Zeroes eventplane
[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
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