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