]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
select only primaries
[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
f1d945a1 23#include "TComplex.h" //needed as include
11d3944d 24#include "TProfile.h" //needed as include
f1d945a1 25
26class TH1F;
11d3944d 27class TFile;
28class TList;
29class TVector2;
f1d945a1 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
40class 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
49ClassImp(AliFlowAnalysisWithLYZEventPlane)
50
51 //-----------------------------------------------------------------------
52
53 AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
9d062fe3 54 fHistList(NULL),
55 fSecondRunList(NULL),
56 fSecondReDtheta(NULL),
57 fSecondImDtheta(NULL),
58 fFirstr0theta(NULL),
80e93fef 59 fHistProVetaRP(NULL),
60 fHistProVetaPOI(NULL),
61 fHistProVPtRP(NULL),
62 fHistProVPtPOI(NULL),
9d062fe3 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)
f1d945a1 78{
79
80 // Constructor.
0092f3c2 81 fQsum = new TVector2(); // flow vector sum
28ca24ad 82
83 fHistList = new TList();
9d062fe3 84 fSecondRunList = new TList();
f1d945a1 85}
86
87
88
89 //-----------------------------------------------------------------------
90
91
92 AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
93 {
94 //destructor
0092f3c2 95 delete fQsum;
28ca24ad 96 delete fHistList;
9d062fe3 97 delete fSecondRunList;
f1d945a1 98 }
99
100
1dfa3c16 101//-----------------------------------------------------------------------
102
103void 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
b0fda271 112//-----------------------------------------------------------------------
113
114void 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
f1d945a1 123//-----------------------------------------------------------------------
124void AliFlowAnalysisWithLYZEventPlane::Init() {
125
126 //Initialise all histograms
127 cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
128
129 //input histograms
9d062fe3 130 if (fSecondRunList) {
131 fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
28ca24ad 132 fHistList->Add(fSecondReDtheta);
9d062fe3 133
134 fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
28ca24ad 135 fHistList->Add(fSecondImDtheta);
9d062fe3 136
137 fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
28ca24ad 138 fHistList->Add(fFirstr0theta);
f1d945a1 139
9d062fe3 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; }
f1d945a1 144
9d062fe3 145 }
28ca24ad 146
9d062fe3 147 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
148 fHistList->Add(fCommonHists);
f1d945a1 149
9d062fe3 150 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
151 fHistList->Add(fCommonHistsRes);
152
882ffd6a 153 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
80e93fef 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);
f1d945a1 179
f1d945a1 180 fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
181 fHistProWr->SetXTitle("Q");
182 fHistProWr->SetYTitle("Wr");
28ca24ad 183 fHistList->Add(fHistProWr);
f1d945a1 184
9d062fe3 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
f1d945a1 190 fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
191 fHistDeltaPhi->SetXTitle("DeltaPhi");
192 fHistDeltaPhi->SetYTitle("Counts");
28ca24ad 193 fHistList->Add(fHistDeltaPhi);
f1d945a1 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");
28ca24ad 198 fHistList->Add(fHistPhiLYZ);
f1d945a1 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");
28ca24ad 203 fHistList->Add(fHistPhiEP);
f1d945a1 204
205 fEventNumber = 0; //set number of events to zero
f1d945a1 206
207}
208
209//-----------------------------------------------------------------------
210
882ffd6a 211void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
212
f1d945a1 213 //Get the event plane and weight for each event
8232a5ec 214 if (anEvent) {
f1d945a1 215
216 //fill control histograms
8232a5ec 217 fCommonHists->FillControlHistograms(anEvent);
f1d945a1 218
219 //get the Q vector from the FlowEvent
882ffd6a 220 AliFlowVector vQ = anEvent->GetQ();
9d062fe3 221 if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
f1d945a1 222 //Weight with the multiplicity
882ffd6a 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);
9d062fe3 230 //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
28ca24ad 231
f1d945a1 232 //for chi calculation:
882ffd6a 233 *fQsum += vQ;
9d062fe3 234 fHistQsumforChi->SetBinContent(1,fQsum->X());
235 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 236 fQ2sum += vQ.Mod2();
9d062fe3 237 fHistQsumforChi->SetBinContent(3,fQ2sum);
238 //cout<<"fQ2sum = "<<fQ2sum<<endl;
f1d945a1 239
240 //call AliFlowLYZEventPlane::CalculateRPandW() here!
882ffd6a 241 aLYZEP->CalculateRPandW(vQ);
f1d945a1 242
882ffd6a 243 Double_t dWR = aLYZEP->GetWR();
244 Double_t dRP = aLYZEP->GetPsi();
f1d945a1 245
882ffd6a 246 //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
247 fHistPhiLYZ->Fill(dRP);
f1d945a1 248
249 //plot difference between event plane from EP-method and LYZ-method
882ffd6a 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);
f1d945a1 254
882ffd6a 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);
f1d945a1 259
260 //Flip sign of WR
882ffd6a 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;
f1d945a1 267 }
882ffd6a 268 fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
f1d945a1 269
80e93fef 270 //calculate flow for RP and POI selections
f1d945a1 271 //loop over the tracks of the event
882ffd6a 272 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
273 for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1 274 {
882ffd6a 275 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
276 if (pTrack){
80e93fef 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));
b7cb54d5 283 if (pTrack->InRPSelection()) {
80e93fef 284 //fill histograms for RP selection
285 fHistProVetaRP -> Fill(dEta,dv2);
286 fHistProVPtRP -> Fill(dPt,dv2);
287 }
b7cb54d5 288 if (pTrack->InPOISelection()) {
80e93fef 289 //fill histograms for POI selection
290 fHistProVetaPOI -> Fill(dEta,dv2);
291 fHistProVPtPOI -> Fill(dPt,dv2);
f1d945a1 292 }
80e93fef 293 }//track
f1d945a1 294 }//loop over tracks
295
296 fEventNumber++;
deebd72f 297 // cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
f1d945a1 298 }
299}
300
301 //--------------------------------------------------------------------
302void AliFlowAnalysisWithLYZEventPlane::Finish() {
303
304 //plot histograms etc.
bc92c0cb 305 cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
9d062fe3 306
307 //constants:
882ffd6a 308 Double_t dJ01 = 2.405;
80e93fef 309 Int_t iNtheta = AliFlowLYZConstants::kTheta;
310 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
311 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
9d062fe3 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
882ffd6a 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;
f1d945a1 327 }
882ffd6a 328 dV /= iNtheta;
f1d945a1 329
882ffd6a 330 //calculate dChi
331 Double_t dSigma2 = 0;
332 Double_t dChi= 0;
f1d945a1 333 if (fEventNumber!=0) {
0092f3c2 334 *fQsum /= fEventNumber;
bc92c0cb 335 //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
336 //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
f1d945a1 337 fQ2sum /= fEventNumber;
bc92c0cb 338 //cerr<<"fEventNumber = "<<fEventNumber<<endl;
339 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a 340 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
bc92c0cb 341 //cerr<<"dSigma2"<<dSigma2<<endl;
882ffd6a 342 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
343 else dChi = -1.;
80e93fef 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 ;
81bbfdbc 366 fCommonHistsRes->FillIntegratedFlow(dV, dVErr);
80e93fef 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;
81bbfdbc 374 cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
80e93fef 375 cout<<endl;
376
f1d945a1 377 }
378
80e93fef 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.) {
882ffd6a 389 for (Int_t theta=0;theta<iNtheta;theta++) {
390 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
80e93fef 391 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 392 TMath::Cos(dTheta));
80e93fef 393 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 394 TMath::Cos(dTheta));
80e93fef 395 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
882ffd6a 396 TMath::BesselJ1(dJ01)))*
80e93fef 397 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
398 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
f1d945a1 399 } //loop over theta
80e93fef 400 }
f1d945a1 401
80e93fef 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
f1d945a1 441
80e93fef 442 //v as a function of Pt for RP selection
b7cb54d5 443 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
81bbfdbc 444 Double_t dVRP = 0.;
445 Double_t dSum = 0.;
446 Double_t dErrV =0.;
447
80e93fef 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;
f1d945a1 472 }
80e93fef 473 else {dErrdifcomb = 0.;}
474
475 //fill TH1D
81bbfdbc 476 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
477 //calculate integrated flow for RP selection
b7cb54d5 478 if (fHistPtRP){
479 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
81bbfdbc 480 dVRP += dv2pro*dYieldPt;
481 dSum +=dYieldPt;
482 dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
b7cb54d5 483 } else { cout<<"fHistPtRP is NULL"<<endl; }
81bbfdbc 484
80e93fef 485 } //loop over bins b
81bbfdbc 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
80e93fef 497
498 //v as a function of Pt for POI selection
b7cb54d5 499 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
80e93fef 500 Double_t dVPOI = 0.;
81bbfdbc 501 dSum = 0.;
502 dErrV =0.;
80e93fef 503
504 for(Int_t b=0;b<iNbinsPt;b++) {
505 Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
506 Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
f1d945a1 507
80e93fef 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
b7cb54d5 540 if (fHistPtPOI){
541 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
80e93fef 542 dVPOI += dv2pro*dYieldPt;
543 dSum +=dYieldPt;
544 dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
b7cb54d5 545 } else { cout<<"fHistPtPOI is NULL"<<endl; }
80e93fef 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);
f1d945a1 554
80e93fef 555 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
556 cout<<endl;
557 cout<<"*************************************"<<endl;
558 cout<<"*************************************"<<endl;
559
560 //cout<<".....finished"<<endl;
f1d945a1 561 }
882ffd6a 562