git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
changing default of SetAcceptKinkDaughters to kTRUE as it should be
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithLYZEventPlane.cxx
CommitLineData
f1d945a1 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"
28ca24ad 25#include "TList.h"
f1d945a1 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
9d062fe3 30#include "TVector2.h"
f1d945a1 31
32class 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
43class 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
52ClassImp(AliFlowAnalysisWithLYZEventPlane)
53
54 //-----------------------------------------------------------------------
55
56 AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
9d062fe3 57 fHistList(NULL),
58 fSecondRunList(NULL),
59 fSecondReDtheta(NULL),
60 fSecondImDtheta(NULL),
61 fFirstr0theta(NULL),
80e93fef 62 fHistProVetaRP(NULL),
63 fHistProVetaPOI(NULL),
64 fHistProVPtRP(NULL),
65 fHistProVPtPOI(NULL),
9d062fe3 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)
f1d945a1 81{
82
83 // Constructor.
0092f3c2 84 fQsum = new TVector2(); // flow vector sum
28ca24ad 85
86 fHistList = new TList();
9d062fe3 87 fSecondRunList = new TList();
f1d945a1 88}
89
90
91
92 //-----------------------------------------------------------------------
93
94
95 AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
96 {
97 //destructor
0092f3c2 98 delete fQsum;
28ca24ad 99 delete fHistList;
9d062fe3 100 delete fSecondRunList;
f1d945a1 101 }
102
103
1dfa3c16 104//-----------------------------------------------------------------------
105
106void 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
b0fda271 115//-----------------------------------------------------------------------
116
117void 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
f1d945a1 126//-----------------------------------------------------------------------
127void AliFlowAnalysisWithLYZEventPlane::Init() {
128
129 //Initialise all histograms
130 cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
131
132 //input histograms
9d062fe3 133 if (fSecondRunList) {
134 fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
28ca24ad 135 fHistList->Add(fSecondReDtheta);
9d062fe3 136
137 fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
28ca24ad 138 fHistList->Add(fSecondImDtheta);
9d062fe3 139
140 fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
28ca24ad 141 fHistList->Add(fFirstr0theta);
f1d945a1 142
9d062fe3 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; }
f1d945a1 147
9d062fe3 148 }
28ca24ad 149
9d062fe3 150 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
151 fHistList->Add(fCommonHists);
f1d945a1 152
9d062fe3 153 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
154 fHistList->Add(fCommonHistsRes);
155
882ffd6a 156 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
80e93fef 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);
f1d945a1 182
f1d945a1 183 fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
184 fHistProWr->SetXTitle("Q");
185 fHistProWr->SetYTitle("Wr");
28ca24ad 186 fHistList->Add(fHistProWr);
f1d945a1 187
9d062fe3 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
f1d945a1 193 fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
194 fHistDeltaPhi->SetXTitle("DeltaPhi");
195 fHistDeltaPhi->SetYTitle("Counts");
28ca24ad 196 fHistList->Add(fHistDeltaPhi);
f1d945a1 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");
28ca24ad 201 fHistList->Add(fHistPhiLYZ);
f1d945a1 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");
28ca24ad 206 fHistList->Add(fHistPhiEP);
f1d945a1 207
208 fEventNumber = 0; //set number of events to zero
f1d945a1 209
210}
211
212//-----------------------------------------------------------------------
213
882ffd6a 214void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
215
f1d945a1 216 //Get the event plane and weight for each event
8232a5ec 217 if (anEvent) {
f1d945a1 218
219 //fill control histograms
8232a5ec 220 fCommonHists->FillControlHistograms(anEvent);
f1d945a1 221
222 //get the Q vector from the FlowEvent
882ffd6a 223 AliFlowVector vQ = anEvent->GetQ();
9d062fe3 224 if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
f1d945a1 225 //Weight with the multiplicity
882ffd6a 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);
9d062fe3 233 //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
28ca24ad 234
f1d945a1 235 //for chi calculation:
882ffd6a 236 *fQsum += vQ;
9d062fe3 237 fHistQsumforChi->SetBinContent(1,fQsum->X());
238 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 239 fQ2sum += vQ.Mod2();
9d062fe3 240 fHistQsumforChi->SetBinContent(3,fQ2sum);
241 //cout<<"fQ2sum = "<<fQ2sum<<endl;
f1d945a1 242
243 //call AliFlowLYZEventPlane::CalculateRPandW() here!
882ffd6a 244 aLYZEP->CalculateRPandW(vQ);
f1d945a1 245
882ffd6a 246 Double_t dWR = aLYZEP->GetWR();
247 Double_t dRP = aLYZEP->GetPsi();
f1d945a1 248
882ffd6a 249 //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
250 fHistPhiLYZ->Fill(dRP);
f1d945a1 251
252 //plot difference between event plane from EP-method and LYZ-method
882ffd6a 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);
f1d945a1 257
882ffd6a 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);
f1d945a1 262
263 //Flip sign of WR
882ffd6a 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;
f1d945a1 270 }
882ffd6a 271 fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
f1d945a1 272
80e93fef 273 //calculate flow for RP and POI selections
f1d945a1 274 //loop over the tracks of the event
882ffd6a 275 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
276 for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1 277 {
882ffd6a 278 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
279 if (pTrack){
80e93fef 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 }
882ffd6a 291 if (pTrack->UseForDifferentialFlow()) {
80e93fef 292 //fill histograms for POI selection
293 fHistProVetaPOI -> Fill(dEta,dv2);
294 fHistProVPtPOI -> Fill(dPt,dv2);
f1d945a1 295 }
80e93fef 296 }//track
f1d945a1 297 }//loop over tracks
298
299 fEventNumber++;
deebd72f 300 // cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
f1d945a1 301 }
302}
303
304 //--------------------------------------------------------------------
305void AliFlowAnalysisWithLYZEventPlane::Finish() {
306
307 //plot histograms etc.
bc92c0cb 308 cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
9d062fe3 309
310 //constants:
882ffd6a 311 Double_t dJ01 = 2.405;
80e93fef 312 Int_t iNtheta = AliFlowLYZConstants::kTheta;
313 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
314 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
9d062fe3 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
882ffd6a 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;
f1d945a1 330 }
882ffd6a 331 dV /= iNtheta;
f1d945a1 332
882ffd6a 333 //calculate dChi
334 Double_t dSigma2 = 0;
335 Double_t dChi= 0;
f1d945a1 336 if (fEventNumber!=0) {
0092f3c2 337 *fQsum /= fEventNumber;
bc92c0cb 338 //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
339 //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
f1d945a1 340 fQ2sum /= fEventNumber;
bc92c0cb 341 //cerr<<"fEventNumber = "<<fEventNumber<<endl;
342 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a 343 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
bc92c0cb 344 //cerr<<"dSigma2"<<dSigma2<<endl;
882ffd6a 345 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
346 else dChi = -1.;
80e93fef 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
f1d945a1 380 }
381
80e93fef 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.) {
882ffd6a 392 for (Int_t theta=0;theta<iNtheta;theta++) {
393 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
80e93fef 394 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 395 TMath::Cos(dTheta));
80e93fef 396 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 397 TMath::Cos(dTheta));
80e93fef 398 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
882ffd6a 399 TMath::BesselJ1(dJ01)))*
80e93fef 400 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
401 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
f1d945a1 402 } //loop over theta
80e93fef 403 }
f1d945a1 404
80e93fef 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
f1d945a1 444
80e93fef 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;
f1d945a1 470 }
80e93fef 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);
f1d945a1 486
80e93fef 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);
f1d945a1 533
80e93fef 534 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
535 cout<<endl;
536 cout<<"*************************************"<<endl;
537 cout<<"*************************************"<<endl;
538
539 //cout<<".....finished"<<endl;
f1d945a1 540 }
882ffd6a 541