]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Base/AliFlowAnalysisWithLeeYangZeros.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- HF can take now all kind of histograms
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / AliFlowAnalysisWithLeeYangZeros.cxx
CommitLineData
f1d945a1 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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
80f72ed6 16/////////////////////////////////////////////////////////////////////////////////////////
17//Description: Maker to analyze Flow using the LeeYangZeros method
18// Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
19// Practical Guide (PG): J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
20// Adapted from StFlowLeeYangZerosMaker.cxx
21// by Markus Oldenberg and Art Poskanzer, LBNL
22// with advice from Jean-Yves Ollitrault and Nicolas Borghini
23//
24//Author: Naomi van der Kolk (kolk@nikhef.nl)
25/////////////////////////////////////////////////////////////////////////////////////////
f1d945a1 26
e77f0120 27#include "Riostream.h"
28#include "TObject.h"
29#include "TMath.h"
30#include "TFile.h"
88e00a8a 31#include "TList.h"
e77f0120 32#include "TVector2.h"
33#include "TH1F.h"
34#include "TComplex.h"
35#include "TProfile.h"
36#include "TDirectoryFile.h"
37
38#include "AliFlowCommonConstants.h"
39#include "AliFlowLYZConstants.h"
40#include "AliFlowAnalysisWithLeeYangZeros.h"
41#include "AliFlowLYZHist1.h"
42#include "AliFlowLYZHist2.h"
43#include "AliFlowCommonHist.h"
44#include "AliFlowCommonHistResults.h"
45#include "AliFlowEventSimple.h"
46#include "AliFlowTrackSimple.h"
47#include "AliFlowVector.h"
f1d945a1 48
3a7af7bd 49using std::endl;
50using std::cout;
51using std::cerr;
f1d945a1 52ClassImp(AliFlowAnalysisWithLeeYangZeros)
53
54 //-----------------------------------------------------------------------
55
56 AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
882ffd6a 57 fQsum(NULL),
58 fQ2sum(0),
59 fEventNumber(0),
60 fFirstRun(kTRUE),
61 fUseSum(kTRUE),
62 fDoubleLoop(kFALSE),
63 fDebug(kFALSE),
88e00a8a 64 fHistList(NULL),
9d062fe3 65 fFirstRunList(NULL),
b9ac26e9 66 fHistVtheta(NULL),
e8c3ff94 67 fHistProVetaRP(NULL),
68 fHistProVetaPOI(NULL),
69 fHistProVPtRP(NULL),
70 fHistProVPtPOI(NULL),
b9ac26e9 71 fHistR0theta(NULL),
882ffd6a 72 fHistProReDenom(NULL),
73 fHistProImDenom(NULL),
b9ac26e9 74 fHistReDtheta(NULL),
75 fHistImDtheta(NULL),
9d062fe3 76 fHistQsumforChi(NULL),
882ffd6a 77 fCommonHists(NULL),
78 fCommonHistsRes(NULL)
f1d945a1 79
80{
81 //default constructor
82 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
83
88e00a8a 84 fHistList = new TList();
9d062fe3 85 fFirstRunList = new TList();
f1d945a1 86
87 for(Int_t i = 0;i<5;i++)
88 {
89 fHist1[i]=0;
e8c3ff94 90 fHist2RP[i]=0;
91 fHist2POI[i]=0;
f1d945a1 92 }
93
0092f3c2 94 fQsum = new TVector2();
f1d945a1 95
96}
97
f1d945a1 98//-----------------------------------------------------------------------
99
100
101 AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
102 {
103 //default destructor
104 if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
0092f3c2 105 delete fQsum;
88e00a8a 106 delete fHistList;
9d062fe3 107 delete fFirstRunList;
108
f1d945a1 109 }
110
7ebf0358 111//-----------------------------------------------------------------------
112
113void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString* outputFileName)
114{
115 //store the final results in output .root file
ebaded83 116
117 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
118 if (GetFirstRun()) {
0fe80f88 119 //output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
c741f5d0 120 if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
121 else {fHistList->SetName("cobjLYZ1PROD");}
9455e15e 122 fHistList->SetOwner(kTRUE);
0fe80f88 123 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
ebaded83 124 }
6a209f0c 125 else {
0fe80f88 126 //output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
c741f5d0 127 if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
128 else { fHistList->SetName("cobjLYZ2PROD"); }
9455e15e 129 fHistList->SetOwner(kTRUE);
0fe80f88 130 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
6a209f0c 131 }
132 delete output;
133}
134
135//-----------------------------------------------------------------------
136
137void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString outputFileName)
138{
139 //store the final results in output .root file
140
141 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
142 if (GetFirstRun()) {
0fe80f88 143 //output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
c741f5d0 144 if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
145 else {fHistList->SetName("cobjLYZ1PROD");}
9455e15e 146 fHistList->SetOwner(kTRUE);
0fe80f88 147 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
6a209f0c 148 }
ebaded83 149 else {
0fe80f88 150 //output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
c741f5d0 151 if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
152 else { fHistList->SetName("cobjLYZ2PROD"); }
9455e15e 153 fHistList->SetOwner(kTRUE);
0fe80f88 154 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
ebaded83 155 }
156 delete output;
7ebf0358 157}
158
159//-----------------------------------------------------------------------
f1d945a1 160
ad87ae62 161void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TDirectoryFile *outputFileName)
162{
163 //store the final results in output .root file
164 if (GetFirstRun()) {
165 if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
166 else {fHistList->SetName("cobjLYZ1PROD");}
167 fHistList->SetOwner(kTRUE);
168 outputFileName->Add(fHistList);
169 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
170 }
171 else {
172 if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
173 else { fHistList->SetName("cobjLYZ2PROD"); }
174 fHistList->SetOwner(kTRUE);
175 outputFileName->Add(fHistList);
176 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
177 }
178}
179
180//-----------------------------------------------------------------------
181
f1d945a1 182Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
183{
184 //init method
185 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
186
489d5531 187 //save old value and prevent histograms from being added to directory
188 //to avoid name clashes in case multiple analaysis objects are used
189 //in an analysis
190 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
191 TH1::AddDirectory(kFALSE);
192
f1d945a1 193 // Book histograms
bee2efdc 194 Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
195 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
196 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
197
198 Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
199 Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
200 Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
201 Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
9d062fe3 202
203 //for control histograms
0c380f17 204 if (fFirstRun){
205 if (fUseSum) {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1SUM");}
206 else {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1PROD");}
207 fHistList->Add(fCommonHists);
208 if (fUseSum) {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1SUM");}
209 else {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1PROD");}
210 fHistList->Add(fCommonHistsRes);
7ebf0358 211 }
0c380f17 212 else {
213 if (fUseSum) {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2SUM");}
214 else {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2PROD");}
215 fHistList->Add(fCommonHists);
216 if (fUseSum) {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2SUM");}
217 else {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2PROD");}
218 fHistList->Add(fCommonHistsRes);
7ebf0358 219 }
0c380f17 220
221 TString nameChiHist;
222 if (fUseSum) {nameChiHist = "Flow_QsumforChi_LYZSUM";}
223 else {nameChiHist = "Flow_QsumforChi_LYZPROD";}
224 fHistQsumforChi = new TH1F(nameChiHist.Data(),nameChiHist.Data(),3,-1.,2.);
9d062fe3 225 fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
226 fHistQsumforChi->SetYTitle("value");
227 fHistList->Add(fHistQsumforChi);
f1d945a1 228
229 //for first loop over events
230 if (fFirstRun){
0c380f17 231 TString nameR0Hist;
b9ac26e9 232 if (fUseSum) {nameR0Hist = "First_Flow_r0theta_LYZSUM";}
233 else {nameR0Hist = "First_Flow_r0theta_LYZPROD";}
234 fHistR0theta = new TH1D(nameR0Hist.Data(),nameR0Hist.Data(),iNtheta,-0.5,iNtheta-0.5);
235 fHistR0theta->SetXTitle("#theta");
236 fHistR0theta->SetYTitle("r_{0}^{#theta}");
237 fHistList->Add(fHistR0theta);
f1d945a1 238
0c380f17 239 TString nameVHist;
b9ac26e9 240 if (fUseSum) {nameVHist = "First_Flow_Vtheta_LYZSUM";}
241 else {nameVHist = "First_Flow_Vtheta_LYZPROD";}
242 fHistVtheta = new TH1D(nameVHist.Data(),nameVHist.Data(),iNtheta,-0.5,iNtheta-0.5);
243 fHistVtheta->SetXTitle("#theta");
244 fHistVtheta->SetYTitle("V_{n}^{#theta}");
245 fHistList->Add(fHistVtheta);
246
247 //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta
882ffd6a 248 for (Int_t theta=0;theta<iNtheta;theta++) {
0c380f17 249 TString nameHist1;
250 if (fUseSum) {nameHist1 = "AliFlowLYZHist1_";}
251 else {nameHist1 = "AliFlowLYZHist1_";}
252 nameHist1 += theta;
253 fHist1[theta]=new AliFlowLYZHist1(theta, nameHist1, fUseSum);
9d062fe3 254 fHistList->Add(fHist1[theta]);
f1d945a1 255 }
9d062fe3 256
f1d945a1 257 }
258 //for second loop over events
259 else {
0c380f17 260 TString nameReDenomHist;
261 if (fUseSum) {nameReDenomHist = "Second_FlowPro_ReDenom_LYZSUM";}
262 else {nameReDenomHist = "Second_FlowPro_ReDenom_LYZPROD";}
263 fHistProReDenom = new TProfile(nameReDenomHist.Data(),nameReDenomHist.Data(), iNtheta, -0.5, iNtheta-0.5);
f1d945a1 264 fHistProReDenom->SetXTitle("#theta");
265 fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
88e00a8a 266 fHistList->Add(fHistProReDenom);
f1d945a1 267
0c380f17 268 TString nameImDenomHist;
269 if (fUseSum) {nameImDenomHist = "Second_FlowPro_ImDenom_LYZSUM";}
270 else {nameImDenomHist = "Second_FlowPro_ImDenom_LYZPROD";}
271 fHistProImDenom = new TProfile(nameImDenomHist.Data(),nameImDenomHist.Data(), iNtheta, -0.5, iNtheta-0.5);
f1d945a1 272 fHistProImDenom->SetXTitle("#theta");
273 fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
88e00a8a 274 fHistList->Add(fHistProImDenom);
f1d945a1 275
0c380f17 276 TString nameVetaRPHist;
277 if (fUseSum) {nameVetaRPHist = "Second_FlowPro_VetaRP_LYZSUM";}
278 else {nameVetaRPHist = "Second_FlowPro_VetaRP_LYZPROD";}
279 fHistProVetaRP = new TProfile(nameVetaRPHist.Data(),nameVetaRPHist.Data(),iNbinsEta,dEtaMin,dEtaMax);
e8c3ff94 280 fHistProVetaRP->SetXTitle("rapidity");
281 fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
282 fHistList->Add(fHistProVetaRP);
283
0c380f17 284 TString nameVetaPOIHist;
285 if (fUseSum) {nameVetaPOIHist = "Second_FlowPro_VetaPOI_LYZSUM";}
286 else {nameVetaPOIHist = "Second_FlowPro_VetaPOI_LYZPROD";}
287 fHistProVetaPOI = new TProfile(nameVetaPOIHist.Data(),nameVetaPOIHist.Data(),iNbinsEta,dEtaMin,dEtaMax);
e8c3ff94 288 fHistProVetaPOI->SetXTitle("rapidity");
289 fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
290 fHistList->Add(fHistProVetaPOI);
291
0c380f17 292 TString nameVPtRPHist;
293 if (fUseSum) {nameVPtRPHist = "Second_FlowPro_VPtRP_LYZSUM";}
294 else {nameVPtRPHist = "Second_FlowPro_VPtRP_LYZPROD";}
295 fHistProVPtRP = new TProfile(nameVPtRPHist.Data(),nameVPtRPHist.Data(),iNbinsPt,dPtMin,dPtMax);
e8c3ff94 296 fHistProVPtRP->SetXTitle("Pt");
297 fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
298 fHistList->Add(fHistProVPtRP);
f1d945a1 299
0c380f17 300 TString nameVPtPOIHist;
301 if (fUseSum) {nameVPtPOIHist = "Second_FlowPro_VPtPOI_LYZSUM";}
302 else {nameVPtPOIHist = "Second_FlowPro_VPtPOI_LYZPROD";}
303 fHistProVPtPOI = new TProfile(nameVPtPOIHist.Data(),nameVPtPOIHist.Data(),iNbinsPt,dPtMin,dPtMax);
e8c3ff94 304 fHistProVPtPOI->SetXTitle("p_{T}");
305 fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
306 fHistList->Add(fHistProVPtPOI);
b9ac26e9 307
0c380f17 308 if (fUseSum){
b9ac26e9 309 fHistReDtheta = new TH1D("Second_Flow_ReDtheta_LYZSUM","Second_Flow_ReDtheta_LYZSUM",iNtheta, -0.5, (double)iNtheta-0.5);
310 fHistReDtheta->SetXTitle("#theta");
311 fHistReDtheta->SetYTitle("Re(D^{#theta})");
312 fHistList->Add(fHistReDtheta);
313
314 fHistImDtheta = new TH1D("Second_Flow_ImDtheta_LYZSUM","Second_Flow_ImDtheta_LYZSUM",iNtheta, -0.5, (double)iNtheta-0.5);
315 fHistImDtheta->SetXTitle("#theta");
316 fHistImDtheta->SetYTitle("Im(D^{#theta})");
317 fHistList->Add(fHistImDtheta);
0c380f17 318 }
f1d945a1 319
320 //class AliFlowLYZHist2 defines the histograms:
882ffd6a 321 for (Int_t theta=0;theta<iNtheta;theta++) {
e8c3ff94 322 TString nameRP = "AliFlowLYZHist2RP_";
323 nameRP += theta;
0c380f17 324 fHist2RP[theta]=new AliFlowLYZHist2(theta, "RP", nameRP, fUseSum);
e8c3ff94 325 fHistList->Add(fHist2RP[theta]);
326
327 TString namePOI = "AliFlowLYZHist2POI_";
328 namePOI += theta;
0c380f17 329 fHist2POI[theta]=new AliFlowLYZHist2(theta, "POI", namePOI, fUseSum);
e8c3ff94 330 fHistList->Add(fHist2POI[theta]);
f1d945a1 331 }
9d062fe3 332
b9ac26e9 333 //read histogram fHistR0theta from the first run list
9d062fe3 334 if (fFirstRunList) {
b9ac26e9 335 if (fUseSum) { fHistR0theta = (TH1D*)fFirstRunList->FindObject("First_Flow_r0theta_LYZSUM");}
336 else{ fHistR0theta = (TH1D*)fFirstRunList->FindObject("First_Flow_r0theta_LYZPROD");}
337 if (!fHistR0theta) {cout<<"fHistR0theta has a NULL pointer!"<<endl;}
338 fHistList->Add(fHistR0theta);
9d062fe3 339 } else { cout<<"list is NULL pointer!"<<endl; }
340
f1d945a1 341 }
342
343
344 if (fDebug) cout<<"****Histograms initialised****"<<endl;
345
346 fEventNumber = 0; //set event counter to zero
489d5531 347
348 //resore old stuff
349 TH1::AddDirectory(oldHistAddStatus);
350
f1d945a1 351 return kTRUE;
352}
353
354 //-----------------------------------------------------------------------
355
0092f3c2 356Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* anEvent)
f1d945a1 357{
358 //make method
359 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
360
361 //get tracks from event
0092f3c2 362 if (anEvent) {
f1d945a1 363 if (fFirstRun){
0092f3c2 364 fCommonHists->FillControlHistograms(anEvent);
365 FillFromFlowEvent(anEvent);
f1d945a1 366 }
367 else {
0092f3c2 368 fCommonHists->FillControlHistograms(anEvent);
369 SecondFillFromFlowEvent(anEvent);
f1d945a1 370 }
371 }
372
373 else {
374 cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
375 return kFALSE;
376 }
951817d5 377 // cout<<"^^^^read event "<<fEventNumber<<endl;
f1d945a1 378 fEventNumber++;
379
380
381 return kTRUE;
382}
383
48385911 384 //-----------------------------------------------------------------------
385void AliFlowAnalysisWithLeeYangZeros::GetOutputHistograms(TList *outputListHistos) {
386 // get the pointers to all output histograms before calling Finish()
387
bee2efdc 388 const Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
48385911 389
390 if (outputListHistos) {
391
392 //define histograms for first and second run
393 AliFlowCommonHist *pCommonHist = NULL;
394 AliFlowCommonHistResults *pCommonHistResults = NULL;
b9ac26e9 395 TH1D* pHistR0theta = NULL;
396 TH1D* pHistVtheta = NULL;
48385911 397 TProfile* pHistProReDenom = NULL;
398 TProfile* pHistProImDenom = NULL;
0c380f17 399 TProfile* pHistProVetaRP = NULL;
48385911 400 TProfile* pHistProVetaPOI = NULL;
0c380f17 401 TProfile* pHistProVPtRP = NULL;
48385911 402 TProfile* pHistProVPtPOI = NULL;
0c380f17 403 TH1F* pHistQsumforChi = NULL;
bee2efdc 404 //AliFlowLYZHist1 *pLYZHist1[iNtheta] = {NULL}; //array of pointers to AliFlowLYZHist1
405 //AliFlowLYZHist2 *pLYZHist2RP[iNtheta] = {NULL}; //array of pointers to AliFlowLYZHist2
406 //AliFlowLYZHist2 *pLYZHist2POI[iNtheta] = {NULL}; //array of pointers to AliFlowLYZHist2
407 AliFlowLYZHist1 **pLYZHist1 = new AliFlowLYZHist1*[iNtheta]; //array of pointers to AliFlowLYZHist1
408 AliFlowLYZHist2 **pLYZHist2RP = new AliFlowLYZHist2*[iNtheta]; //array of pointers to AliFlowLYZHist2
409 AliFlowLYZHist2 **pLYZHist2POI = new AliFlowLYZHist2*[iNtheta]; //array of pointers to AliFlowLYZHist2
410 for (Int_t i=0; i<iNtheta; i++)
411 {
412 pLYZHist1[i] = NULL;
413 pLYZHist2RP[i] = NULL;
414 pLYZHist2POI[i] = NULL;
415 }
48385911 416
417 if (GetFirstRun()) { //first run
418 //Get the common histograms from the output list
0c380f17 419 if (GetUseSum()){
420 pCommonHist = dynamic_cast<AliFlowCommonHist*>
421 (outputListHistos->FindObject("AliFlowCommonHistLYZ1SUM"));
422 pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
423 (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ1SUM"));
424 //Get the histograms from the output list
425 for(Int_t theta = 0;theta<iNtheta;theta++){
426 TString name = "AliFlowLYZHist1_";
427 name += theta;
428 pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
429 (outputListHistos->FindObject(name));
430 }
b9ac26e9 431 pHistVtheta = dynamic_cast<TH1D*>
432 (outputListHistos->FindObject("First_Flow_Vtheta_LYZSUM"));
0c380f17 433
b9ac26e9 434 pHistR0theta = dynamic_cast<TH1D*>
435 (outputListHistos->FindObject("First_Flow_r0theta_LYZSUM"));
0c380f17 436
437 pHistQsumforChi = dynamic_cast<TH1F*>
438 (outputListHistos->FindObject("Flow_QsumforChi_LYZSUM"));
439
440 //Set the histogram pointers and call Finish()
441 if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
b9ac26e9 442 pHistVtheta && pHistR0theta && pHistQsumforChi ) {
0c380f17 443 this->SetCommonHists(pCommonHist);
444 this->SetCommonHistsRes(pCommonHistResults);
445 this->SetHist1(pLYZHist1);
b9ac26e9 446 this->SetHistVtheta(pHistVtheta);
447 this->SetHistR0theta(pHistR0theta);
0c380f17 448 this->SetHistQsumforChi(pHistQsumforChi);
449 }
450 else {
451 cout<<"WARNING: Histograms needed to run Finish() firstrun (SUM) are not accessable!"<<endl;
452 }
453 }
454 else {
455 pCommonHist = dynamic_cast<AliFlowCommonHist*>
456 (outputListHistos->FindObject("AliFlowCommonHistLYZ1PROD"));
457 pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
458 (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ1PROD"));
459 //Get the histograms from the output list
460 for(Int_t theta = 0;theta<iNtheta;theta++){
461 TString name = "AliFlowLYZHist1_";
462 name += theta;
463 pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
464 (outputListHistos->FindObject(name));
465 }
b9ac26e9 466 pHistVtheta = dynamic_cast<TH1D*>
467 (outputListHistos->FindObject("First_Flow_Vtheta_LYZPROD"));
0c380f17 468
b9ac26e9 469 pHistR0theta = dynamic_cast<TH1D*>
470 (outputListHistos->FindObject("First_Flow_r0theta_LYZPROD"));
0c380f17 471
472 pHistQsumforChi = dynamic_cast<TH1F*>
473 (outputListHistos->FindObject("Flow_QsumforChi_LYZPROD"));
474
475 //Set the histogram pointers and call Finish()
476 if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
b9ac26e9 477 pHistVtheta && pHistR0theta && pHistQsumforChi ) {
0c380f17 478 this->SetCommonHists(pCommonHist);
479 this->SetCommonHistsRes(pCommonHistResults);
480 this->SetHist1(pLYZHist1);
b9ac26e9 481 this->SetHistVtheta(pHistVtheta);
482 this->SetHistR0theta(pHistR0theta);
0c380f17 483 this->SetHistQsumforChi(pHistQsumforChi);
484 } else {
485 cout<<"WARNING: Histograms needed to run Finish() firstrun (PROD) are not accessable!"<<endl;
486 }
487 }
48385911 488 }
489 else { //second run
490 //Get the common histograms from the output list
0c380f17 491 if (GetUseSum()){
492 pCommonHist = dynamic_cast<AliFlowCommonHist*>
493 (outputListHistos->FindObject("AliFlowCommonHistLYZ2SUM"));
494 pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
495 (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ2SUM"));
496
b9ac26e9 497 pHistR0theta = dynamic_cast<TH1D*>
498 (outputListHistos->FindObject("First_Flow_r0theta_LYZSUM"));
0c380f17 499
500 pHistQsumforChi = dynamic_cast<TH1F*>
501 (outputListHistos->FindObject("Flow_QsumforChi_LYZSUM"));
502
503 //Get the histograms from the output list
504 for(Int_t theta = 0;theta<iNtheta;theta++){
505 TString nameRP = "AliFlowLYZHist2RP_";
506 nameRP += theta;
507 pLYZHist2RP[theta] = dynamic_cast<AliFlowLYZHist2*>
508 (outputListHistos->FindObject(nameRP));
509 TString namePOI = "AliFlowLYZHist2POI_";
510 namePOI += theta;
511 pLYZHist2POI[theta] = dynamic_cast<AliFlowLYZHist2*>
512 (outputListHistos->FindObject(namePOI));
513 }
514 pHistProReDenom = dynamic_cast<TProfile*>
515 (outputListHistos->FindObject("Second_FlowPro_ReDenom_LYZSUM"));
516 pHistProImDenom = dynamic_cast<TProfile*>
517 (outputListHistos->FindObject("Second_FlowPro_ImDenom_LYZSUM"));
518
b9ac26e9 519 TH1D* pHistReDtheta = dynamic_cast<TH1D*>
520 (outputListHistos->FindObject("Second_Flow_ReDtheta_LYZSUM"));
521 TH1D* pHistImDtheta = dynamic_cast<TH1D*>
522 (outputListHistos->FindObject("Second_Flow_ImDtheta_LYZSUM"));
0c380f17 523
524 pHistProVetaRP = dynamic_cast<TProfile*>
525 (outputListHistos->FindObject("Second_FlowPro_VetaRP_LYZSUM"));
526 pHistProVetaPOI = dynamic_cast<TProfile*>
527 (outputListHistos->FindObject("Second_FlowPro_VetaPOI_LYZSUM"));
528 pHistProVPtRP = dynamic_cast<TProfile*>
529 (outputListHistos->FindObject("Second_FlowPro_VPtRP_LYZSUM"));
530 pHistProVPtPOI = dynamic_cast<TProfile*>
531 (outputListHistos->FindObject("Second_FlowPro_VPtPOI_LYZSUM"));
532
533
534 //Set the histogram pointers and call Finish()
b9ac26e9 535 if (pCommonHist && pCommonHistResults &&
536 pLYZHist2RP[0] && pLYZHist2POI[0] &&
537 pHistR0theta &&
538 pHistProReDenom && pHistProImDenom &&
539 pHistReDtheta && pHistImDtheta &&
540 pHistProVetaRP && pHistProVetaPOI &&
541 pHistProVPtRP && pHistProVPtPOI &&
542 pHistQsumforChi) {
0c380f17 543 this->SetCommonHists(pCommonHist);
544 this->SetCommonHistsRes(pCommonHistResults);
545 this->SetHist2RP(pLYZHist2RP);
546 this->SetHist2POI(pLYZHist2POI);
b9ac26e9 547 this->SetHistR0theta(pHistR0theta);
0c380f17 548 this->SetHistProReDenom(pHistProReDenom);
549 this->SetHistProImDenom(pHistProImDenom);
b9ac26e9 550 this->SetHistReDtheta(pHistReDtheta);
551 this->SetHistImDtheta(pHistImDtheta);
0c380f17 552 this->SetHistProVetaRP(pHistProVetaRP);
553 this->SetHistProVetaPOI(pHistProVetaPOI);
554 this->SetHistProVPtRP(pHistProVPtRP);
555 this->SetHistProVPtPOI(pHistProVPtPOI);
556 this->SetHistQsumforChi(pHistQsumforChi);
557 }
558 else {
559 cout<<"WARNING: Histograms needed to run Finish() secondrun (SUM) are not accessable!"<<endl;
560 }
48385911 561 }
0c380f17 562 else {
563 pCommonHist = dynamic_cast<AliFlowCommonHist*>
564 (outputListHistos->FindObject("AliFlowCommonHistLYZ2PROD"));
565 pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
566 (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ2PROD"));
48385911 567
0c380f17 568
b9ac26e9 569 pHistR0theta = dynamic_cast<TH1D*>
570 (outputListHistos->FindObject("First_Flow_r0theta_LYZPROD"));
0c380f17 571
572 pHistQsumforChi = dynamic_cast<TH1F*>
573 (outputListHistos->FindObject("Flow_QsumforChi_LYZPROD"));
574
575 //Get the histograms from the output list
576 for(Int_t theta = 0;theta<iNtheta;theta++){
577 TString nameRP = "AliFlowLYZHist2RP_";
578 nameRP += theta;
579 pLYZHist2RP[theta] = dynamic_cast<AliFlowLYZHist2*>
580 (outputListHistos->FindObject(nameRP));
581 TString namePOI = "AliFlowLYZHist2POI_";
582 namePOI += theta;
583 pLYZHist2POI[theta] = dynamic_cast<AliFlowLYZHist2*>
584 (outputListHistos->FindObject(namePOI));
585 }
586 pHistProReDenom = dynamic_cast<TProfile*>
587 (outputListHistos->FindObject("Second_FlowPro_ReDenom_LYZPROD"));
588 pHistProImDenom = dynamic_cast<TProfile*>
589 (outputListHistos->FindObject("Second_FlowPro_ImDenom_LYZPROD"));
590
591 pHistProVetaRP = dynamic_cast<TProfile*>
592 (outputListHistos->FindObject("Second_FlowPro_VetaRP_LYZPROD"));
593 pHistProVetaPOI = dynamic_cast<TProfile*>
594 (outputListHistos->FindObject("Second_FlowPro_VetaPOI_LYZPROD"));
595 pHistProVPtRP = dynamic_cast<TProfile*>
596 (outputListHistos->FindObject("Second_FlowPro_VPtRP_LYZPROD"));
597 pHistProVPtPOI = dynamic_cast<TProfile*>
598 (outputListHistos->FindObject("Second_FlowPro_VPtPOI_LYZPROD"));
599
600 //Set the histogram pointers and call Finish()
601 if (pCommonHist && pCommonHistResults && pLYZHist2RP[0] && pLYZHist2POI[0] &&
b9ac26e9 602 pHistR0theta && pHistProReDenom && pHistProImDenom && pHistProVetaRP &&
0c380f17 603 pHistProVetaPOI && pHistProVPtRP && pHistProVPtPOI && pHistQsumforChi) {
604 this->SetCommonHists(pCommonHist);
605 this->SetCommonHistsRes(pCommonHistResults);
606 this->SetHist2RP(pLYZHist2RP);
607 this->SetHist2POI(pLYZHist2POI);
b9ac26e9 608 this->SetHistR0theta(pHistR0theta);
0c380f17 609 this->SetHistProReDenom(pHistProReDenom);
610 this->SetHistProImDenom(pHistProImDenom);
611 this->SetHistProVetaRP(pHistProVetaRP);
612 this->SetHistProVetaPOI(pHistProVetaPOI);
613 this->SetHistProVPtRP(pHistProVPtRP);
614 this->SetHistProVPtPOI(pHistProVPtPOI);
615 this->SetHistQsumforChi(pHistQsumforChi);
616 }
617 else {
618 cout<<"WARNING: Histograms needed to run Finish() secondrun (PROD) are not accessable!"<<endl;
619 }
48385911 620 }
0c380f17 621 } //secondrun
622 //outputListHistos->Print();
bee2efdc 623 delete [] pLYZHist1;
624 delete [] pLYZHist2RP;
625 delete [] pLYZHist2POI;
0c380f17 626 } //listhistos
627 else {
628 cout << "histogram list pointer is empty in method AliFlowAnalysisWithLeeYangZeros::GetOutputHistograms() " << endl;}
48385911 629}
0c380f17 630
f1d945a1 631 //-----------------------------------------------------------------------
632 Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
633{
634 //finish method
635 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
e8c3ff94 636
f1d945a1 637 //define variables for both runs
882ffd6a 638 Double_t dJ01 = 2.405;
bee2efdc 639 Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
b9ac26e9 640
9d062fe3 641 //set the event number
a93de0f0 642 SetEventNumber((int)fCommonHists->GetHistQ()->GetEntries());
b9ac26e9 643
5b40431d 644 //Get multiplicity for RP selection
645 Double_t dMultRP = fCommonHists->GetHistMultRP()->GetMean();
646 if (fDebug) cout<<"The average multiplicity is "<<dMultRP<<endl;
647
9d062fe3 648 //set the sum of Q vectors
649 fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
650 SetQ2sum(fHistQsumforChi->GetBinContent(3));
651
f1d945a1 652 if (fFirstRun){
b9ac26e9 653
654 //define variables for the first run
882ffd6a 655 Double_t dR0 = 0;
656 Double_t dVtheta = 0;
657 Double_t dv = 0;
658 Double_t dV = 0;
b9ac26e9 659
660 //reset histograms in case of merged output files
661 fHistR0theta->Reset();
662 fHistVtheta->Reset();
663
882ffd6a 664 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 665 {
666 //get the first minimum r0
667 fHist1[theta]->FillGtheta();
882ffd6a 668 dR0 = fHist1[theta]->GetR0();
4930f9f3 669
f1d945a1 670 //calculate integrated flow
e77f0120 671 if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) { dVtheta = dJ01/dR0; }
4930f9f3 672 else { cout<<"r0 is not found! Leaving LYZ analysis."<<endl; return kFALSE; }
673
882ffd6a 674 //for estimating systematic error resulting from d0
aa5e7c3f 675 Double_t dBinsize =0.;
bee2efdc 676 if (fUseSum){ dBinsize = (AliFlowLYZConstants::GetMaster()->GetMaxSUM())/(AliFlowLYZConstants::GetMaster()->GetNbins());}
677 else { dBinsize = (AliFlowLYZConstants::GetMaster()->GetMaxPROD())/(AliFlowLYZConstants::GetMaster()->GetNbins());}
3201e2d2 678 Double_t dVplus = -1.;
679 Double_t dVmin = -1.;
e77f0120 680 if (!TMath::AreEqualAbs(dR0+dBinsize, 0., 1e-100)) {dVplus = dJ01/(dR0+dBinsize);}
681 if (!TMath::AreEqualAbs(dR0-dBinsize, 0., 1e-100)) {dVmin = dJ01/(dR0-dBinsize);}
5b40431d 682 //convert V to v
683 Double_t dvplus = -1.;
684 Double_t dvmin= -1.;
685 if (fUseSum){
686 //for SUM: V=v because the Q-vector is scaled by 1/M
687 dv = dVtheta;
688 dvplus = dVplus;
689 dvmin = dVmin; }
690 else {
691 //for PRODUCT: v=V/M
e77f0120 692 if (!TMath::AreEqualAbs(dMultRP, 0., 1e-100)){
5b40431d 693 dv = dVtheta/dMultRP;
694 dvplus = dVplus/dMultRP;
695 dvmin = dVmin/dMultRP; }}
696
e8c3ff94 697 if (fDebug) cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
f1d945a1 698
699 //fill the histograms
b9ac26e9 700 fHistR0theta->SetBinContent(theta+1,dR0);
701 fHistR0theta->SetBinError(theta+1,0.0);
702 fHistVtheta->SetBinContent(theta+1,dVtheta);
703 fHistVtheta->SetBinError(theta+1,0.0);
f1d945a1 704 //get average value of fVtheta = fV
882ffd6a 705 dV += dVtheta;
f1d945a1 706 } //end of loop over theta
707
708 //get average value of fVtheta = fV
882ffd6a 709 dV /=iNtheta;
5b40431d 710 if (!fUseSum) { if (dMultRP!=0.){dV /=dMultRP;}} //scale with multiplicity for PRODUCT
aa5e7c3f 711
f1d945a1 712 //sigma2 and chi
882ffd6a 713 Double_t dSigma2 = 0;
714 Double_t dChi= 0;
f1d945a1 715 if (fEventNumber!=0) {
0092f3c2 716 *fQsum /= fEventNumber;
aa5e7c3f 717 //cout<<"fQsum is "<<fQsum->X()<<" "<<fQsum->Y()<<endl;
f1d945a1 718 fQ2sum /= fEventNumber;
aa5e7c3f 719 //cout<<"fQ2sum is "<<fQ2sum<<endl;
882ffd6a 720 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
aa5e7c3f 721 //cout<<"dSigma2 is "<<dSigma2<<endl;
882ffd6a 722 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
723 else dChi = -1.;
62e36168 724 fCommonHistsRes->FillChi(dChi);
e8c3ff94 725
726 cout<<"*************************************"<<endl;
727 cout<<"*************************************"<<endl;
728 cout<<" Integrated flow from "<<endl;
c741f5d0 729 if (fUseSum) {
730 cout<<" Lee-Yang Zeroes SUM "<<endl;}
731 else {
732 cout<<" Lee-Yang Zeroes PRODUCT "<<endl;}
e8c3ff94 733 cout<<endl;
81bbfdbc 734 cout<<"Chi = "<<dChi<<endl;
5b40431d 735 //cout<<endl;
f1d945a1 736 }
737
738 // recalculate statistical errors on integrated flow
739 //combining 5 theta angles to 1 relative error BP eq. 89
882ffd6a 740 Double_t dRelErr2comb = 0.;
9d062fe3 741 Int_t iEvts = fEventNumber;
1eda9b5e 742 if (iEvts!=0) {
743 for (Int_t theta=0;theta<iNtheta;theta++){
744 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
745 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 746 TMath::Cos(dTheta));
1eda9b5e 747 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 748 TMath::Cos(dTheta));
1eda9b5e 749 dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
882ffd6a 750 TMath::BesselJ1(dJ01)))*
1eda9b5e 751 (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
752 dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
753 }
754 dRelErr2comb /= iNtheta;
f1d945a1 755 }
882ffd6a 756 Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
f1d945a1 757
758 //copy content of profile into TH1D and add error
882ffd6a 759 Double_t dv2pro = dV; //in the case that fv is equal to fV
760 Double_t dv2Err = dv2pro*dRelErrcomb ;
81bbfdbc 761 cout<<"dV = "<<dv2pro<<" +- "<<dv2Err<<endl;
e8c3ff94 762 cout<<endl;
763 cout<<"*************************************"<<endl;
764 cout<<"*************************************"<<endl;
81bbfdbc 765 fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
f1d945a1 766
767
768 if (fDebug) cout<<"****histograms filled****"<<endl;
769
f1d945a1 770 return kTRUE;
f1d945a1 771 } //firstrun
772
773
b9ac26e9 774 else { //second run to calculate differential flow
f1d945a1 775
b9ac26e9 776 //declare variables for the second run
e8c3ff94 777 TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
f1d945a1 778 Int_t m = 1;
779 TComplex i = TComplex::I();
882ffd6a 780 Double_t dBesselRatio[3] = {1., 1.202, 2.69};
bee2efdc 781 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
782 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
9d062fe3 783
882ffd6a 784 Double_t dR0 = 0.;
785 Double_t dVtheta = 0.;
786 Double_t dV = 0.;
787 Double_t dReDenom = 0.;
788 Double_t dImDenom = 0.;
b9ac26e9 789
790 //reset histograms in case of merged output files
791 if (fUseSum) {
792 fHistReDtheta->Reset();
793 fHistImDtheta->Reset();
794 }
795 fHistProVetaRP ->Reset();
796 fHistProVetaPOI->Reset();
797 fHistProVPtRP ->Reset();
798 fHistProVPtPOI ->Reset();
799
800 //scale fHistR0theta by the number of merged files to undo the merging
801 if (!fHistR0theta) { cout<<"Hist pointer R0theta in file does not exist"<<endl; }
802 else {
803 Int_t iEntries = (int)fHistR0theta->GetEntries();
804 if (iEntries > iNtheta){
805 //for each individual file fHistR0theta has iNtheta entries
806 Int_t iFiles = iEntries/iNtheta;
807 cout<<iFiles<<" files were merged!"<<endl;
808 fHistR0theta->Scale(1./iFiles);
809 }
810
811 //loop over theta
812 for (Int_t theta=0;theta<iNtheta;theta++) {
813 dR0 = fHistR0theta->GetBinContent(theta+1); //histogram starts at bin 1
882ffd6a 814 if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
e77f0120 815 if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) dVtheta = dJ01/dR0;
882ffd6a 816 dV += dVtheta;
f1d945a1 817 // BP Eq. 9 -> Vn^theta = j01/r0^theta
24696303 818 if (!fHistProReDenom || !fHistProImDenom) {
e8c3ff94 819 cout << "Hist pointer fDenom in file does not exist" <<endl;
24696303 820 cout<< "Leaving LYZ second pass analysis!" <<endl;
821 return kFALSE;
e8c3ff94 822 } else {
882ffd6a 823 dReDenom = fHistProReDenom->GetBinContent(theta+1);
824 dImDenom = fHistProImDenom->GetBinContent(theta+1);
e8c3ff94 825 cDenom(dReDenom,dImDenom);
b9ac26e9 826
e8c3ff94 827 //for new method and use by others (only with the sum generating function):
828 if (fUseSum) {
829 cDtheta = dR0*cDenom/dJ01;
b9ac26e9 830 fHistReDtheta->SetBinContent(theta+1,cDtheta.Re());
831 fHistReDtheta->SetBinError(theta+1,0.0);
832 fHistImDtheta->SetBinContent(theta+1,cDtheta.Im());
833 fHistImDtheta->SetBinError(theta+1,0.0);
e8c3ff94 834 }
f1d945a1 835
e8c3ff94 836 cDenom *= TComplex::Power(i, m-1);
b9ac26e9 837
e8c3ff94 838 //v as a function of eta for RP selection
839 for (Int_t be=1;be<=iNbinsEta;be++) {
b9ac26e9 840 Double_t dReRatioRP = 0.0;
841 Double_t dEtaRP = fHist2RP[theta]->GetBinCenter(be);
e8c3ff94 842 cNumerRP = fHist2RP[theta]->GetNumerEta(be);
843 if (cNumerRP.Rho()==0) {
844 if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
e8c3ff94 845 }
e77f0120 846 else if (TMath::AreEqualAbs(cDenom.Rho(), 0, 1e-100)) {
e8c3ff94 847 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
e8c3ff94 848 }
849 else {
850 dReRatioRP = (cNumerRP/cDenom).Re();
851 }
b9ac26e9 852 Double_t dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
e8c3ff94 853 fHistProVetaRP->Fill(dEtaRP,dVetaRP);
854 } //loop over bins be
f1d945a1 855
e8c3ff94 856
857 //v as a function of eta for POI selection
858 for (Int_t be=1;be<=iNbinsEta;be++) {
b9ac26e9 859 Double_t dReRatioPOI = 0.0;
860 Double_t dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
e8c3ff94 861 cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
862 if (cNumerPOI.Rho()==0) {
863 if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
e8c3ff94 864 }
865 else if (cDenom.Rho()==0) {
866 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
e8c3ff94 867 }
868 else {
869 dReRatioPOI = (cNumerPOI/cDenom).Re();
870 }
b9ac26e9 871 Double_t dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
e8c3ff94 872 fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
873 } //loop over bins be
874
875
876 //v as a function of Pt for RP selection
877 for (Int_t bp=1;bp<=iNbinsPt;bp++) {
b9ac26e9 878 Double_t dReRatioRP = 0.0;
879 Double_t dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
e8c3ff94 880 cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
881 if (cNumerRP.Rho()==0) {
882 if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
e8c3ff94 883 }
884 else if (cDenom.Rho()==0) {
885 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
e8c3ff94 886 }
887 else {
888 dReRatioRP = (cNumerRP/cDenom).Re();
889 }
b9ac26e9 890 Double_t dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
e8c3ff94 891 fHistProVPtRP->Fill(dPtRP,dVPtRP);
892 } //loop over bins bp
893
894
e8c3ff94 895 //v as a function of Pt for POI selection
896 for (Int_t bp=1;bp<=iNbinsPt;bp++) {
b9ac26e9 897 Double_t dReRatioPOI = 0.0;
898 Double_t dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
e8c3ff94 899 cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
900 if (cNumerPOI.Rho()==0) {
901 if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
e8c3ff94 902 }
903 else if (cDenom.Rho()==0) {
904 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
e8c3ff94 905 }
906 else {
907 dReRatioPOI = (cNumerPOI/cDenom).Re();
908 }
b9ac26e9 909 Double_t dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
e8c3ff94 910 fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
911 } //loop over bins bp
912
f1d945a1 913 }
e8c3ff94 914 }
f1d945a1 915
916 }//end of loop over theta
917
918 //calculate the average of fVtheta = fV
882ffd6a 919 dV /= iNtheta;
5b40431d 920 if (!fUseSum) { if (dMultRP!=0.) { dV /=dMultRP; }} //scale by the multiplicity for PRODUCT
e77f0120 921 if (TMath::AreEqualAbs(dV, 0., 1e-100)) { cout<<"dV = 0! Leaving LYZ analysis."<<endl; return kFALSE; }
f1d945a1 922
923 //sigma2 and chi (for statistical error calculations)
882ffd6a 924 Double_t dSigma2 = 0;
925 Double_t dChi= 0;
f1d945a1 926 if (fEventNumber!=0) {
0092f3c2 927 *fQsum /= fEventNumber;
f1d945a1 928 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
929 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
930 fQ2sum /= fEventNumber;
931 //cout<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a 932 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
933 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
934 else dChi = -1.;
62e36168 935 fCommonHistsRes->FillChi(dChi);
e8c3ff94 936
937 // recalculate statistical errors on integrated flow
938 //combining 5 theta angles to 1 relative error BP eq. 89
939 Double_t dRelErr2comb = 0.;
940 Int_t iEvts = fEventNumber;
941 if (iEvts!=0) {
942 for (Int_t theta=0;theta<iNtheta;theta++){
943 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
944 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
945 TMath::Cos(dTheta));
946 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
947 TMath::Cos(dTheta));
948 dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
949 TMath::BesselJ1(dJ01)))*
950 (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
951 dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
952 }
953 dRelErr2comb /= iNtheta;
954 }
955 Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
956 Double_t dVErr = dV*dRelErrcomb ;
81bbfdbc 957 fCommonHistsRes->FillIntegratedFlow(dV,dVErr);
e8c3ff94 958
959 cout<<"*************************************"<<endl;
960 cout<<"*************************************"<<endl;
961 cout<<" Integrated flow from "<<endl;
c741f5d0 962 if (fUseSum) {
963 cout<<" Lee-Yang Zeroes SUM "<<endl;}
964 else {
965 cout<<" Lee-Yang Zeroes PRODUCT "<<endl;}
e8c3ff94 966 cout<<endl;
81bbfdbc 967 cout<<"Chi = "<<dChi<<endl;
968 cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
5b40431d 969 //cout<<endl;
f1d945a1 970 }
971
e8c3ff94 972 //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
973
974 //v as a function of eta for RP selection
975 for(Int_t b=0;b<iNbinsEta;b++) {
976 Double_t dv2pro = fHistProVetaRP->GetBinContent(b);
977 Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
978 Double_t dErrdifcomb = 0.; //set error to zero
979 Double_t dErr2difcomb = 0.; //set error to zero
980 //calculate error
e77f0120 981 if (!TMath::AreEqualAbs(dNprime, 0., 1e-100)) {
e8c3ff94 982 for (Int_t theta=0;theta<iNtheta;theta++) {
983 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
984 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
985 TMath::Cos(dTheta));
986 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
987 TMath::Cos(dTheta));
988 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
989 TMath::BesselJ1(dJ01)))*
990 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
991 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
992 } //loop over theta
993 }
994
e77f0120 995 if (!TMath::AreEqualAbs(dErr2difcomb, 0., 1e-100)) {
e8c3ff94 996 dErr2difcomb /= iNtheta;
997 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
998 }
999 else {dErrdifcomb = 0.;}
1000 //fill TH1D
1001 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
1002 } //loop over bins b
1003
1004
1005 //v as a function of eta for POI selection
1006 for(Int_t b=0;b<iNbinsEta;b++) {
1007 Double_t dv2pro = fHistProVetaPOI->GetBinContent(b);
1008 Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
1009 Double_t dErrdifcomb = 0.; //set error to zero
1010 Double_t dErr2difcomb = 0.; //set error to zero
1011 //calculate error
1012 if (dNprime!=0.) {
1013 for (Int_t theta=0;theta<iNtheta;theta++) {
1014 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1015 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1016 TMath::Cos(dTheta));
1017 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1018 TMath::Cos(dTheta));
1019 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1020 TMath::BesselJ1(dJ01)))*
1021 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1022 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1023 } //loop over theta
1024 }
1025
1026 if (dErr2difcomb!=0.) {
1027 dErr2difcomb /= iNtheta;
1028 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1029 }
1030 else {dErrdifcomb = 0.;}
1031 //fill TH1D
1032 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
1033 } //loop over bins b
1034
1035
1036
1037 //v as a function of Pt for RP selection
b7cb54d5 1038 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
81bbfdbc 1039 Double_t dVRP = 0.;
1040 Double_t dSum = 0.;
1041 Double_t dErrV =0.;
882ffd6a 1042 for(Int_t b=0;b<iNbinsPt;b++) {
e8c3ff94 1043 Double_t dv2pro = fHistProVPtRP->GetBinContent(b);
1044 Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
1045 Double_t dErrdifcomb = 0.; //set error to zero
1046 Double_t dErr2difcomb = 0.; //set error to zero
1047 //calculate error
1048 if (dNprime!=0.) {
1049 for (Int_t theta=0;theta<iNtheta;theta++) {
1050 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1051 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1052 TMath::Cos(dTheta));
1053 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1054 TMath::Cos(dTheta));
1055 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1056 TMath::BesselJ1(dJ01)))*
1057 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1058 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1059 } //loop over theta
1060 }
1061
1062 if (dErr2difcomb!=0.) {
1063 dErr2difcomb /= iNtheta;
1064 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1065 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1066 }
1067 else {dErrdifcomb = 0.;}
1068
1069 //fill TH1D
1070 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
81bbfdbc 1071 //calculate integrated flow for RP selection
b7cb54d5 1072 if (fHistPtRP){
1073 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
81bbfdbc 1074 dVRP += dv2pro*dYieldPt;
1075 dSum +=dYieldPt;
1076 dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
b7cb54d5 1077 } else { cout<<"fHistPtRP is NULL"<<endl; }
e8c3ff94 1078 } //loop over bins b
81bbfdbc 1079
e77f0120 1080 if (!TMath::AreEqualAbs(dSum, 0., 1e-100)) {
81bbfdbc 1081 dVRP /= dSum; //the pt distribution should be normalised
1082 dErrV /= (dSum*dSum);
1083 dErrV = TMath::Sqrt(dErrV);
1084 }
1085 cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
1086 //cout<<endl;
1087 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
1088
1089
e8c3ff94 1090 //v as a function of Pt for POI selection
b7cb54d5 1091 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
e8c3ff94 1092 Double_t dVPOI = 0.;
81bbfdbc 1093 dSum = 0.;
1094 dErrV =0.;
e8c3ff94 1095
1096 for(Int_t b=0;b<iNbinsPt;b++) {
1097 Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
1098 Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
882ffd6a 1099 Double_t dErrdifcomb = 0.; //set error to zero
1100 Double_t dErr2difcomb = 0.; //set error to zero
f1d945a1 1101 //calculate error
882ffd6a 1102 if (dNprime!=0.) {
1103 for (Int_t theta=0;theta<iNtheta;theta++) {
1104 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1105 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1106 TMath::Cos(dTheta));
1107 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1108 TMath::Cos(dTheta));
1109 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1110 TMath::BesselJ1(dJ01)))*
1111 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1112 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
f1d945a1 1113 } //loop over theta
1114 }
1115
882ffd6a 1116 if (dErr2difcomb!=0.) {
1117 dErr2difcomb /= iNtheta;
0828acf6 1118 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
882ffd6a 1119 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
f1d945a1 1120 }
882ffd6a 1121 else {dErrdifcomb = 0.;}
f1d945a1 1122
1123 //fill TH1D
e8c3ff94 1124 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
1125 //calculate integrated flow for POI selection
b7cb54d5 1126 if (fHistPtPOI){
1127 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
e8c3ff94 1128 dVPOI += dv2pro*dYieldPt;
1129 dSum +=dYieldPt;
1130 dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
b7cb54d5 1131 } else { cout<<"fHistPtPOI is NULL"<<endl; }
f1d945a1 1132 } //loop over bins b
e8c3ff94 1133
1134 if (dSum != 0.) {
1135 dVPOI /= dSum; //the pt distribution should be normalised
1136 dErrV /= (dSum*dSum);
1137 dErrV = TMath::Sqrt(dErrV);
1138 }
81bbfdbc 1139 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
e8c3ff94 1140 cout<<endl;
1141 cout<<"*************************************"<<endl;
1142 cout<<"*************************************"<<endl;
1143 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
1144
f1d945a1 1145 } //secondrun
1146
951817d5 1147 //cout<<"----LYZ analysis finished....----"<<endl<<endl;
f1d945a1 1148
1149 return kTRUE;
1150}
1151
1152
1153//-----------------------------------------------------------------------
1154
0092f3c2 1155 Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
f1d945a1 1156{
1157 // Get event quantities from AliFlowEvent for all particles
1158
1159 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
1160
0092f3c2 1161 if (!anEvent){
f1d945a1 1162 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1163 return kFALSE;
1164 }
1165
1166 //define variables
882ffd6a 1167 TComplex cExpo, cGtheta, cGthetaNew, cZ;
bee2efdc 1168 Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1169 Int_t iNbins = AliFlowLYZConstants::GetMaster()->GetNbins();
b9ac26e9 1170
f1d945a1 1171 //calculate flow
882ffd6a 1172 Double_t dOrder = 2.;
f1d945a1 1173
1174 //get the Q vector
882ffd6a 1175 AliFlowVector vQ = anEvent->GetQ();
f1d945a1 1176 //weight by the multiplicity
882ffd6a 1177 Double_t dQX = 0;
1178 Double_t dQY = 0;
e77f0120 1179 if (!TMath::AreEqualAbs(vQ.GetMult(), 0., 1e-100)) {
882ffd6a 1180 dQX = vQ.X()/vQ.GetMult();
1181 dQY = vQ.Y()/vQ.GetMult();
f2ec07bf 1182 }
882ffd6a 1183 vQ.Set(dQX,dQY);
9d062fe3 1184
f1d945a1 1185 //for chi calculation:
882ffd6a 1186 *fQsum += vQ;
9d062fe3 1187 fHistQsumforChi->SetBinContent(1,fQsum->X());
1188 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 1189 fQ2sum += vQ.Mod2();
9d062fe3 1190 fHistQsumforChi->SetBinContent(3,fQ2sum);
b9ac26e9 1191
1192 for (Int_t theta=0;theta<iNtheta;theta++) {
1193 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 1194
b9ac26e9 1195 //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1196 Double_t dQtheta = GetQtheta(vQ, dTheta);
f1d945a1 1197
b9ac26e9 1198 for (Int_t bin=1;bin<=iNbins;bin++) {
1199 Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED???
1200 if (fUseSum) {
1201 //calculate the sum generating function
1202 cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
1203 cGtheta = TComplex::Exp(cExpo);
1204 }
1205 else {
1206 //calculate the product generating function
1207 cGtheta = GetGrtheta(anEvent, dR, dTheta);
1208 if (cGtheta.Rho2() > 100.) break;
1209 }
1210 //fill real and imaginary part of cGtheta
1211 fHist1[theta]->Fill(dR,cGtheta);
1212 } //loop over bins
1213 } //loop over theta
1214
f1d945a1 1215 return kTRUE;
f1d945a1 1216
1217}
1218
1219 //-----------------------------------------------------------------------
0092f3c2 1220 Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
f1d945a1 1221{
1222 //for differential flow
1223
1224 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
1225
0092f3c2 1226 if (!anEvent){
f1d945a1 1227 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1228 return kFALSE;
1229 }
1230
1231 //define variables
e8c3ff94 1232 TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
882ffd6a 1233 Double_t dR0 = 0.;
e8c3ff94 1234 Double_t dCosTermRP = 0.;
1235 Double_t dCosTermPOI = 0.;
f1d945a1 1236 Double_t m = 1.;
882ffd6a 1237 Double_t dOrder = 2.;
bee2efdc 1238 Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
b9ac26e9 1239
f1d945a1 1240 //get the Q vector
882ffd6a 1241 AliFlowVector vQ = anEvent->GetQ();
f1d945a1 1242 //weight by the multiplicity
882ffd6a 1243 Double_t dQX = 0.;
1244 Double_t dQY = 0.;
1245 if (vQ.GetMult() != 0) {
1246 dQX = vQ.X()/vQ.GetMult();
1247 dQY = vQ.Y()/vQ.GetMult();
f2ec07bf 1248 }
9d062fe3 1249 vQ.Set(dQX,dQY);
1250
f1d945a1 1251 //for chi calculation:
882ffd6a 1252 *fQsum += vQ;
9d062fe3 1253 fHistQsumforChi->SetBinContent(1,fQsum->X());
1254 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 1255 fQ2sum += vQ.Mod2();
9d062fe3 1256 fHistQsumforChi->SetBinContent(3,fQ2sum);
f1d945a1 1257
882ffd6a 1258 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 1259 {
882ffd6a 1260 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 1261
882ffd6a 1262 //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1263 Double_t dQtheta = GetQtheta(vQ, dTheta);
b9ac26e9 1264
f1d945a1 1265 //denominator for differential v
b9ac26e9 1266 if (fHistR0theta) {
1267 dR0 = fHistR0theta->GetBinContent(theta+1);
f1d945a1 1268 }
b9ac26e9 1269 else { cout <<"pointer fHistR0theta does not exist" << endl;
f1d945a1 1270 }
b9ac26e9 1271
e8c3ff94 1272 if (fUseSum) //sum generating function
f1d945a1 1273 {
882ffd6a 1274 cExpo(0.,dR0*dQtheta);
1275 cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
f1d945a1 1276 //loop over tracks in event
882ffd6a 1277 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1278 for (Int_t i=0;i<iNumberOfTracks;i++) {
1279 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i);
1280 if (pTrack) {
b9ac26e9 1281 Double_t dEta = pTrack->Eta();
1282 Double_t dPt = pTrack->Pt();
1283 Double_t dPhi = pTrack->Phi();
b7cb54d5 1284 if (pTrack->InRPSelection()) { // RP selection
e8c3ff94 1285 dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
1286 cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
b9ac26e9 1287 if (cNumerRP.Rho()==0) { cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
1288 if (fDebug) { cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;}
e8c3ff94 1289 if (fHist2RP[theta]) {
b9ac26e9 1290 fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1291 }
f1d945a1 1292 }
b7cb54d5 1293 if (pTrack->InPOISelection()) { //POI selection
e8c3ff94 1294 dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
1295 cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
b9ac26e9 1296 if (cNumerPOI.Rho()==0) { cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
1297 if (fDebug) { cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;}
e8c3ff94 1298 if (fHist2POI[theta]) {
b9ac26e9 1299 fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1300 }
e8c3ff94 1301 }
e8c3ff94 1302 } //if track
f1d945a1 1303 else {cerr << "no particle!!!"<<endl;}
1304 } //loop over tracks
f1d945a1 1305 } //sum
b9ac26e9 1306 else { //product generating function
1307 cDenom = GetDiffFlow(anEvent, dR0, theta);
1308 }//product
1309
f1d945a1 1310 if (fHistProReDenom && fHistProImDenom) {
882ffd6a 1311 fHistProReDenom->Fill(theta,cDenom.Re()); //fill the real part of fDenom
1312 fHistProImDenom->Fill(theta,cDenom.Im()); //fill the imaginary part of fDenom
f1d945a1 1313 }
b9ac26e9 1314 else { cout << "Pointers to cDenom mising" << endl;}
1315
f1d945a1 1316 }//end of loop over theta
1317
1318 return kTRUE;
b9ac26e9 1319
f1d945a1 1320}
1321 //-----------------------------------------------------------------------
882ffd6a 1322 Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
f1d945a1 1323{
882ffd6a 1324 //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
f1d945a1 1325 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
1326
882ffd6a 1327 Double_t dQtheta = 0.;
1328 Double_t dOrder = 2.;
f1d945a1 1329
882ffd6a 1330 dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
f1d945a1 1331
882ffd6a 1332 return dQtheta;
f1d945a1 1333
1334}
1335
1336
1337//-----------------------------------------------------------------------
80f72ed6 1338TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* const anEvent, Double_t aR, Double_t aTheta)
f1d945a1 1339{
1340 // Product Generating Function for LeeYangZeros method
1341 // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1342
1343 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1344
1345
882ffd6a 1346 TComplex cG = TComplex::One();
1347 Double_t dOrder = 2.;
5b40431d 1348 Double_t dWgt = 1.;
1349 //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
8fbe1a98 1350
882ffd6a 1351 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1 1352
882ffd6a 1353 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
f1d945a1 1354 {
882ffd6a 1355 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1356 if (pTrack){
b7cb54d5 1357 if (pTrack->InRPSelection()) {
882ffd6a 1358 Double_t dPhi = pTrack->Phi();
1359 Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
1360 TComplex cGi(1., dGIm);
1361 cG *= cGi; //product over all tracks
f1d945a1 1362 }
1363 }
1364 else {cerr << "no particle pointer !!!"<<endl;}
1365 }//loop over tracks
1366
882ffd6a 1367 return cG;
f1d945a1 1368
1369}
1370
1371
1372//-----------------------------------------------------------------------
80f72ed6 1373TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* const anEvent, Double_t aR0, Int_t theta)
f1d945a1 1374{
1375 // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1376 // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1377 // Also for v1 mixed harmonics: DF Eq. 5
1378 // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1379
1380 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1381
882ffd6a 1382 TComplex cG = TComplex::One();
1383 TComplex cdGr0(0.,0.);
1384 Double_t dOrder = 2.;
1385 Double_t dWgt = 1.;
5b40431d 1386 //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1387
882ffd6a 1388 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1 1389
bee2efdc 1390 Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
882ffd6a 1391 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 1392
e8c3ff94 1393 //for the denominator (use all RP selected particles)
882ffd6a 1394 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
f1d945a1 1395 {
882ffd6a 1396 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1397 if (pTrack){
b7cb54d5 1398 if (pTrack->InRPSelection()) {
882ffd6a 1399 Double_t dPhi = pTrack->Phi();
1400 Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
f1d945a1 1401 //GetGr0theta
882ffd6a 1402 Double_t dGIm = aR0 * dCosTerm;
1403 TComplex cGi(1., dGIm);
e8c3ff94 1404 TComplex cCosTermComplex(1., aR0*dCosTerm);
882ffd6a 1405 cG *= cGi; //product over all tracks
f1d945a1 1406 //GetdGr0theta
882ffd6a 1407 cdGr0 +=(dCosTerm / cCosTermComplex); //sum over all tracks
f1d945a1 1408 }
1409 } //if particle
1410 else {cerr << "no particle!!!"<<endl;}
1411 }//loop over tracks
1412
e8c3ff94 1413 //for the numerator
882ffd6a 1414 for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1 1415 {
882ffd6a 1416 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1417 if (pTrack){
e8c3ff94 1418 Double_t dEta = pTrack->Eta();
1419 Double_t dPt = pTrack->Pt();
1420 Double_t dPhi = pTrack->Phi();
1421 Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1422 TComplex cCosTermComplex(1.,aR0*dCosTerm);
1423 //RP selection
b7cb54d5 1424 if (pTrack->InRPSelection()) {
e8c3ff94 1425 TComplex cNumerRP = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
1426 fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1427 }
1428 //POI selection
b7cb54d5 1429 if (pTrack->InPOISelection()) {
e8c3ff94 1430 TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
1431 fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
f1d945a1 1432 }
1433 } //if particle
1434 else {cerr << "no particle pointer!!!"<<endl;}
1435 }//loop over tracks
1436
e8c3ff94 1437 TComplex cDenom = cG*cdGr0;
882ffd6a 1438 return cDenom;
f1d945a1 1439
1440}
1441
1442//-----------------------------------------------------------------------
1443