]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix product method
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / 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
27#include "Riostream.h" //needed as include
e8c3ff94 28#include "TObject.h" //needed as include
f1d945a1 29#include "AliFlowCommonConstants.h" //needed as include
30#include "AliFlowLYZConstants.h" //needed as include
31#include "AliFlowAnalysisWithLeeYangZeros.h"
32#include "AliFlowLYZHist1.h" //needed as include
33#include "AliFlowLYZHist2.h" //needed as include
34#include "AliFlowCommonHist.h" //needed as include
35#include "AliFlowCommonHistResults.h" //needed as include
36#include "AliFlowEventSimple.h" //needed as include
37#include "AliFlowTrackSimple.h" //needed as include
38
39class AliFlowVector;
40
41#include "TMath.h" //needed as include
42#include "TFile.h" //needed as include
88e00a8a 43#include "TList.h"
f1d945a1 44
45class TComplex;
46class TProfile;
47class TH1F;
48class TH1D;
49
50
f1d945a1 51ClassImp(AliFlowAnalysisWithLeeYangZeros)
52
53 //-----------------------------------------------------------------------
54
55 AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
882ffd6a 56 fQsum(NULL),
57 fQ2sum(0),
58 fEventNumber(0),
59 fFirstRun(kTRUE),
60 fUseSum(kTRUE),
61 fDoubleLoop(kFALSE),
62 fDebug(kFALSE),
88e00a8a 63 fHistList(NULL),
9d062fe3 64 fFirstRunList(NULL),
882ffd6a 65 fHistProVtheta(NULL),
e8c3ff94 66 fHistProVetaRP(NULL),
67 fHistProVetaPOI(NULL),
68 fHistProVPtRP(NULL),
69 fHistProVPtPOI(NULL),
882ffd6a 70 fHistProR0theta(NULL),
71 fHistProReDenom(NULL),
72 fHistProImDenom(NULL),
73 fHistProReDtheta(NULL),
74 fHistProImDtheta(NULL),
9d062fe3 75 fHistQsumforChi(NULL),
882ffd6a 76 fCommonHists(NULL),
77 fCommonHistsRes(NULL)
f1d945a1 78
79{
80 //default constructor
81 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
82
88e00a8a 83 fHistList = new TList();
9d062fe3 84 fFirstRunList = new TList();
f1d945a1 85
86 for(Int_t i = 0;i<5;i++)
87 {
88 fHist1[i]=0;
e8c3ff94 89 fHist2RP[i]=0;
90 fHist2POI[i]=0;
f1d945a1 91 }
92
0092f3c2 93 fQsum = new TVector2();
f1d945a1 94
95}
96
f1d945a1 97//-----------------------------------------------------------------------
98
99
100 AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
101 {
102 //default destructor
103 if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
0092f3c2 104 delete fQsum;
88e00a8a 105 delete fHistList;
9d062fe3 106 delete fFirstRunList;
107
f1d945a1 108 }
109
7ebf0358 110//-----------------------------------------------------------------------
111
112void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString* outputFileName)
113{
114 //store the final results in output .root file
ebaded83 115
116 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
117 if (GetFirstRun()) {
118 output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
119 }
6a209f0c 120 else {
121 output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
122 }
123 delete output;
124}
125
126//-----------------------------------------------------------------------
127
128void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString outputFileName)
129{
130 //store the final results in output .root file
131
132 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
133 if (GetFirstRun()) {
134 output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
135 }
ebaded83 136 else {
137 output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
138 }
139 delete output;
7ebf0358 140}
141
142//-----------------------------------------------------------------------
f1d945a1 143
144Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
145{
146 //init method
147 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
148
f1d945a1 149 // Book histograms
882ffd6a 150 Int_t iNtheta = AliFlowLYZConstants::kTheta;
151 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
152 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
153
154 Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
155 Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
156 Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();
157 Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();
9d062fe3 158
159 //for control histograms
7ebf0358 160 if (fFirstRun){ fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1");
9d062fe3 161 fHistList->Add(fCommonHists);
7ebf0358 162 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1");
163 fHistList->Add(fCommonHistsRes);
164 }
165 else { fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2");
166 fHistList->Add(fCommonHists);
167 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2");
9d062fe3 168 fHistList->Add(fCommonHistsRes);
7ebf0358 169 }
170
171 fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZ","Flow_QsumforChi_LYZ",3,-1.,2.);
9d062fe3 172 fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
173 fHistQsumforChi->SetYTitle("value");
174 fHistList->Add(fHistQsumforChi);
f1d945a1 175
176 //for first loop over events
177 if (fFirstRun){
882ffd6a 178 fHistProR0theta = new TProfile("First_FlowPro_r0theta_LYZ","First_FlowPro_r0theta_LYZ",iNtheta,-0.5,iNtheta-0.5);
f1d945a1 179 fHistProR0theta->SetXTitle("#theta");
180 fHistProR0theta->SetYTitle("r_{0}^{#theta}");
88e00a8a 181 fHistList->Add(fHistProR0theta);
f1d945a1 182
882ffd6a 183 fHistProVtheta = new TProfile("First_FlowPro_Vtheta_LYZ","First_FlowPro_Vtheta_LYZ",iNtheta,-0.5,iNtheta-0.5);
f1d945a1 184 fHistProVtheta->SetXTitle("#theta");
185 fHistProVtheta->SetYTitle("V_{n}^{#theta}");
88e00a8a 186 fHistList->Add(fHistProVtheta);
f1d945a1 187
188 //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
882ffd6a 189 for (Int_t theta=0;theta<iNtheta;theta++) {
9d062fe3 190 TString name = "AliFlowLYZHist1_";
191 name += theta;
192 fHist1[theta]=new AliFlowLYZHist1(theta, name);
193 fHistList->Add(fHist1[theta]);
f1d945a1 194 }
9d062fe3 195
f1d945a1 196 }
197 //for second loop over events
198 else {
882ffd6a 199 fHistProReDenom = new TProfile("Second_FlowPro_ReDenom_LYZ","Second_FlowPro_ReDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5);
f1d945a1 200 fHistProReDenom->SetXTitle("#theta");
201 fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
88e00a8a 202 fHistList->Add(fHistProReDenom);
f1d945a1 203
882ffd6a 204 fHistProImDenom = new TProfile("Second_FlowPro_ImDenom_LYZ","Second_FlowPro_ImDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5);
f1d945a1 205 fHistProImDenom->SetXTitle("#theta");
206 fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
88e00a8a 207 fHistList->Add(fHistProImDenom);
f1d945a1 208
e8c3ff94 209 fHistProVetaRP = new TProfile("Second_FlowPro_VetaRP_LYZ","Second_FlowPro_VetaRP_LYZ",iNbinsEta,dEtaMin,dEtaMax);
210 fHistProVetaRP->SetXTitle("rapidity");
211 fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
212 fHistList->Add(fHistProVetaRP);
213
214 fHistProVetaPOI = new TProfile("Second_FlowPro_VetaPOI_LYZ","Second_FlowPro_VetaPOI_LYZ",iNbinsEta,dEtaMin,dEtaMax);
215 fHistProVetaPOI->SetXTitle("rapidity");
216 fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
217 fHistList->Add(fHistProVetaPOI);
218
219 fHistProVPtRP = new TProfile("Second_FlowPro_VPtRP_LYZ","Second_FlowPro_VPtRP_LYZ",iNbinsPt,dPtMin,dPtMax);
220 fHistProVPtRP->SetXTitle("Pt");
221 fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
222 fHistList->Add(fHistProVPtRP);
f1d945a1 223
e8c3ff94 224 fHistProVPtPOI = new TProfile("Second_FlowPro_VPtPOI_LYZ","Second_FlowPro_VPtPOI_LYZ",iNbinsPt,dPtMin,dPtMax);
225 fHistProVPtPOI->SetXTitle("p_{T}");
226 fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
227 fHistList->Add(fHistProVPtPOI);
f1d945a1 228
882ffd6a 229 fHistProReDtheta = new TProfile("Second_FlowPro_ReDtheta_LYZ","Second_FlowPro_ReDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5);
f1d945a1 230 fHistProReDtheta->SetXTitle("#theta");
231 fHistProReDtheta->SetYTitle("Re(D^{#theta})");
88e00a8a 232 fHistList->Add(fHistProReDtheta);
f1d945a1 233
882ffd6a 234 fHistProImDtheta = new TProfile("Second_FlowPro_ImDtheta_LYZ","Second_FlowPro_ImDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5);
f1d945a1 235 fHistProImDtheta->SetXTitle("#theta");
236 fHistProImDtheta->SetYTitle("Im(D^{#theta})");
88e00a8a 237 fHistList->Add(fHistProImDtheta);
f1d945a1 238
239 //class AliFlowLYZHist2 defines the histograms:
882ffd6a 240 for (Int_t theta=0;theta<iNtheta;theta++) {
e8c3ff94 241 TString nameRP = "AliFlowLYZHist2RP_";
242 nameRP += theta;
81bbfdbc 243 fHist2RP[theta]=new AliFlowLYZHist2(theta, "RP", nameRP);
e8c3ff94 244 fHistList->Add(fHist2RP[theta]);
245
246 TString namePOI = "AliFlowLYZHist2POI_";
247 namePOI += theta;
81bbfdbc 248 fHist2POI[theta]=new AliFlowLYZHist2(theta, "POI", namePOI);
e8c3ff94 249 fHistList->Add(fHist2POI[theta]);
f1d945a1 250 }
9d062fe3 251
252 //read histogram fHistProR0theta from the first run list
253 if (fFirstRunList) {
254 fHistProR0theta = (TProfile*)fFirstRunList->FindObject("First_FlowPro_r0theta_LYZ");
88e00a8a 255 if (!fHistProR0theta) {cout<<"fHistProR0theta has a NULL pointer!"<<endl;}
256 fHistList->Add(fHistProR0theta);
9d062fe3 257 } else { cout<<"list is NULL pointer!"<<endl; }
258
f1d945a1 259 }
260
261
262 if (fDebug) cout<<"****Histograms initialised****"<<endl;
263
264 fEventNumber = 0; //set event counter to zero
265
266 return kTRUE;
267}
268
269 //-----------------------------------------------------------------------
270
0092f3c2 271Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* anEvent)
f1d945a1 272{
273 //make method
274 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
275
276 //get tracks from event
0092f3c2 277 if (anEvent) {
f1d945a1 278 if (fFirstRun){
0092f3c2 279 fCommonHists->FillControlHistograms(anEvent);
280 FillFromFlowEvent(anEvent);
f1d945a1 281 }
282 else {
0092f3c2 283 fCommonHists->FillControlHistograms(anEvent);
284 SecondFillFromFlowEvent(anEvent);
f1d945a1 285 }
286 }
287
288 else {
289 cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
290 return kFALSE;
291 }
951817d5 292 // cout<<"^^^^read event "<<fEventNumber<<endl;
f1d945a1 293 fEventNumber++;
294
295
296 return kTRUE;
297}
298
299
300 //-----------------------------------------------------------------------
301 Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
302{
303 //finish method
304 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
e8c3ff94 305
f1d945a1 306 //define variables for both runs
882ffd6a 307 Double_t dJ01 = 2.405;
308 Int_t iNtheta = AliFlowLYZConstants::kTheta;
9d062fe3 309 //set the event number
310 SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
1eda9b5e 311 //cout<<"number of events processed is "<<fEventNumber<<endl;
9d062fe3 312
313 //set the sum of Q vectors
314 fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
315 SetQ2sum(fHistQsumforChi->GetBinContent(3));
316
f1d945a1 317 if (fFirstRun){
882ffd6a 318 Double_t dR0 = 0;
319 Double_t dVtheta = 0;
320 Double_t dv = 0;
321 Double_t dV = 0;
322 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 323 {
324 //get the first minimum r0
325 fHist1[theta]->FillGtheta();
882ffd6a 326 dR0 = fHist1[theta]->GetR0();
4930f9f3 327
f1d945a1 328 //calculate integrated flow
4930f9f3 329 if (dR0!=0) { dVtheta = dJ01/dR0; }
330 else { cout<<"r0 is not found! Leaving LYZ analysis."<<endl; return kFALSE; }
331
882ffd6a 332 //for estimating systematic error resulting from d0
333 Double_t dBinsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins);
334 Double_t dVplus = dJ01/(dR0+dBinsize);
335 Double_t dVmin = dJ01/(dR0-dBinsize);
336 dv = dVtheta;
337 Double_t dvplus = dVplus;
338 Double_t dvmin = dVmin;
e8c3ff94 339 if (fDebug) cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
f1d945a1 340
341 //fill the histograms
882ffd6a 342 fHistProR0theta->Fill(theta,dR0);
343 fHistProVtheta->Fill(theta,dVtheta);
f1d945a1 344 //get average value of fVtheta = fV
882ffd6a 345 dV += dVtheta;
f1d945a1 346 } //end of loop over theta
347
348 //get average value of fVtheta = fV
882ffd6a 349 dV /=iNtheta;
4930f9f3 350
f1d945a1 351 //sigma2 and chi
882ffd6a 352 Double_t dSigma2 = 0;
353 Double_t dChi= 0;
f1d945a1 354 if (fEventNumber!=0) {
0092f3c2 355 *fQsum /= fEventNumber;
f1d945a1 356 fQ2sum /= fEventNumber;
882ffd6a 357 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
358 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
359 else dChi = -1.;
e8c3ff94 360 fCommonHistsRes->FillChiRP(dChi);
361
362 cout<<"*************************************"<<endl;
363 cout<<"*************************************"<<endl;
364 cout<<" Integrated flow from "<<endl;
365 cout<<" Lee-Yang Zeroes "<<endl;
366 cout<<endl;
81bbfdbc 367 cout<<"Chi = "<<dChi<<endl;
e8c3ff94 368 cout<<endl;
f1d945a1 369 }
370
371 // recalculate statistical errors on integrated flow
372 //combining 5 theta angles to 1 relative error BP eq. 89
882ffd6a 373 Double_t dRelErr2comb = 0.;
9d062fe3 374 Int_t iEvts = fEventNumber;
1eda9b5e 375 if (iEvts!=0) {
376 for (Int_t theta=0;theta<iNtheta;theta++){
377 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
378 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 379 TMath::Cos(dTheta));
1eda9b5e 380 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 381 TMath::Cos(dTheta));
1eda9b5e 382 dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
882ffd6a 383 TMath::BesselJ1(dJ01)))*
1eda9b5e 384 (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
385 dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
386 }
387 dRelErr2comb /= iNtheta;
f1d945a1 388 }
882ffd6a 389 Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
f1d945a1 390
391 //copy content of profile into TH1D and add error
882ffd6a 392 Double_t dv2pro = dV; //in the case that fv is equal to fV
393 Double_t dv2Err = dv2pro*dRelErrcomb ;
81bbfdbc 394 cout<<"dV = "<<dv2pro<<" +- "<<dv2Err<<endl;
e8c3ff94 395 cout<<endl;
396 cout<<"*************************************"<<endl;
397 cout<<"*************************************"<<endl;
81bbfdbc 398 fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
f1d945a1 399
400
401 if (fDebug) cout<<"****histograms filled****"<<endl;
402
f1d945a1 403 return kTRUE;
404 fEventNumber =0; //set to zero for second round over events
405 } //firstrun
406
407
408 else { //second run
409
410 //calculate differential flow
411 //declare variables
e8c3ff94 412 TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
f1d945a1 413 Int_t m = 1;
414 TComplex i = TComplex::I();
882ffd6a 415 Double_t dBesselRatio[3] = {1., 1.202, 2.69};
416 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
417 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
f1d945a1 418
e8c3ff94 419 Double_t dEtaRP, dPtRP, dReRatioRP, dVetaRP, dVPtRP, dEtaPOI, dPtPOI, dReRatioPOI, dVetaPOI, dVPtPOI;
9d062fe3 420
882ffd6a 421 Double_t dR0 = 0.;
422 Double_t dVtheta = 0.;
423 Double_t dV = 0.;
424 Double_t dReDenom = 0.;
425 Double_t dImDenom = 0.;
426 for (Int_t theta=0;theta<iNtheta;theta++) { //loop over theta
e8c3ff94 427 if (!fHistProR0theta) {
428 cout << "Hist pointer R0theta in file does not exist" <<endl;
429 } else {
430 dR0 = fHistProR0theta->GetBinContent(theta+1); //histogram starts at bin 1
882ffd6a 431 if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
432 if (dR0!=0) dVtheta = dJ01/dR0;
433 dV += dVtheta;
f1d945a1 434 // BP Eq. 9 -> Vn^theta = j01/r0^theta
e8c3ff94 435 if (!fHistProReDenom && !fHistProImDenom) {
436 cout << "Hist pointer fDenom in file does not exist" <<endl;
437 } else {
882ffd6a 438 dReDenom = fHistProReDenom->GetBinContent(theta+1);
439 dImDenom = fHistProImDenom->GetBinContent(theta+1);
e8c3ff94 440 cDenom(dReDenom,dImDenom);
f1d945a1 441
e8c3ff94 442 //for new method and use by others (only with the sum generating function):
443 if (fUseSum) {
444 cDtheta = dR0*cDenom/dJ01;
445 fHistProReDtheta->Fill(theta,cDtheta.Re());
446 fHistProImDtheta->Fill(theta,cDtheta.Im());
447 }
f1d945a1 448
e8c3ff94 449 cDenom *= TComplex::Power(i, m-1);
450 //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
f1d945a1 451
e8c3ff94 452 //v as a function of eta for RP selection
453 for (Int_t be=1;be<=iNbinsEta;be++) {
454 dEtaRP = fHist2RP[theta]->GetBinCenter(be);
455 cNumerRP = fHist2RP[theta]->GetNumerEta(be);
456 if (cNumerRP.Rho()==0) {
457 if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
458 dReRatioRP = 0;
459 }
460 else if (cDenom.Rho()==0) {
461 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
462 dReRatioRP = 0;
463 }
464 else {
465 dReRatioRP = (cNumerRP/cDenom).Re();
466 }
467 dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
468 fHistProVetaRP->Fill(dEtaRP,dVetaRP);
469 } //loop over bins be
f1d945a1 470
e8c3ff94 471
472 //v as a function of eta for POI selection
473 for (Int_t be=1;be<=iNbinsEta;be++) {
474 dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
475 cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
476 if (cNumerPOI.Rho()==0) {
477 if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
478 dReRatioPOI = 0;
479 }
480 else if (cDenom.Rho()==0) {
481 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
482 dReRatioPOI = 0;
483 }
484 else {
485 dReRatioPOI = (cNumerPOI/cDenom).Re();
486 }
487 dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
488 fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
489 } //loop over bins be
490
491
492 //v as a function of Pt for RP selection
493 for (Int_t bp=1;bp<=iNbinsPt;bp++) {
494 dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
495 cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
496 if (cNumerRP.Rho()==0) {
497 if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
498 dReRatioRP = 0;
499 }
500 else if (cDenom.Rho()==0) {
501 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
502 dReRatioRP = 0;
503 }
504 else {
505 dReRatioRP = (cNumerRP/cDenom).Re();
506 }
507 dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
508 fHistProVPtRP->Fill(dPtRP,dVPtRP);
509 } //loop over bins bp
510
511
512
513 //v as a function of Pt for POI selection
514 for (Int_t bp=1;bp<=iNbinsPt;bp++) {
515 dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
516 cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
517 if (cNumerPOI.Rho()==0) {
518 if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
519 dReRatioPOI = 0;
520 }
521 else if (cDenom.Rho()==0) {
522 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
523 dReRatioPOI = 0;
524 }
525 else {
526 dReRatioPOI = (cNumerPOI/cDenom).Re();
527 }
528 dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
529 fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
530 } //loop over bins bp
531
f1d945a1 532 }
e8c3ff94 533 }
f1d945a1 534
535 }//end of loop over theta
536
537 //calculate the average of fVtheta = fV
882ffd6a 538 dV /= iNtheta;
4930f9f3 539 if (dV==0.) { cout<<"dV = 0! Leaving LYZ analysis."<<endl; return kFALSE; }
f1d945a1 540
541 //sigma2 and chi (for statistical error calculations)
882ffd6a 542 Double_t dSigma2 = 0;
543 Double_t dChi= 0;
f1d945a1 544 if (fEventNumber!=0) {
0092f3c2 545 *fQsum /= fEventNumber;
f1d945a1 546 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
547 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
548 fQ2sum /= fEventNumber;
549 //cout<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a 550 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
551 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
552 else dChi = -1.;
e8c3ff94 553 fCommonHistsRes->FillChiRP(dChi);
554
555 // recalculate statistical errors on integrated flow
556 //combining 5 theta angles to 1 relative error BP eq. 89
557 Double_t dRelErr2comb = 0.;
558 Int_t iEvts = fEventNumber;
559 if (iEvts!=0) {
560 for (Int_t theta=0;theta<iNtheta;theta++){
561 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
562 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
563 TMath::Cos(dTheta));
564 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
565 TMath::Cos(dTheta));
566 dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
567 TMath::BesselJ1(dJ01)))*
568 (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
569 dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
570 }
571 dRelErr2comb /= iNtheta;
572 }
573 Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
574 Double_t dVErr = dV*dRelErrcomb ;
81bbfdbc 575 fCommonHistsRes->FillIntegratedFlow(dV,dVErr);
e8c3ff94 576
577 cout<<"*************************************"<<endl;
578 cout<<"*************************************"<<endl;
579 cout<<" Integrated flow from "<<endl;
580 cout<<" Lee-Yang Zeroes "<<endl;
581 cout<<endl;
81bbfdbc 582 cout<<"Chi = "<<dChi<<endl;
583 cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
584 cout<<endl;
f1d945a1 585 }
586
e8c3ff94 587 //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
588
589 //v as a function of eta for RP selection
590 for(Int_t b=0;b<iNbinsEta;b++) {
591 Double_t dv2pro = fHistProVetaRP->GetBinContent(b);
592 Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
593 Double_t dErrdifcomb = 0.; //set error to zero
594 Double_t dErr2difcomb = 0.; //set error to zero
595 //calculate error
596 if (dNprime!=0.) {
597 for (Int_t theta=0;theta<iNtheta;theta++) {
598 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
599 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
600 TMath::Cos(dTheta));
601 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
602 TMath::Cos(dTheta));
603 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
604 TMath::BesselJ1(dJ01)))*
605 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
606 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
607 } //loop over theta
608 }
609
610 if (dErr2difcomb!=0.) {
611 dErr2difcomb /= iNtheta;
612 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
613 }
614 else {dErrdifcomb = 0.;}
615 //fill TH1D
616 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
617 } //loop over bins b
618
619
620 //v as a function of eta for POI selection
621 for(Int_t b=0;b<iNbinsEta;b++) {
622 Double_t dv2pro = fHistProVetaPOI->GetBinContent(b);
623 Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
624 Double_t dErrdifcomb = 0.; //set error to zero
625 Double_t dErr2difcomb = 0.; //set error to zero
626 //calculate error
627 if (dNprime!=0.) {
628 for (Int_t theta=0;theta<iNtheta;theta++) {
629 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
630 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
631 TMath::Cos(dTheta));
632 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
633 TMath::Cos(dTheta));
634 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
635 TMath::BesselJ1(dJ01)))*
636 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
637 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
638 } //loop over theta
639 }
640
641 if (dErr2difcomb!=0.) {
642 dErr2difcomb /= iNtheta;
643 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
644 }
645 else {dErrdifcomb = 0.;}
646 //fill TH1D
647 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
648 } //loop over bins b
649
650
651
652 //v as a function of Pt for RP selection
b7cb54d5 653 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
81bbfdbc 654 Double_t dVRP = 0.;
655 Double_t dSum = 0.;
656 Double_t dErrV =0.;
882ffd6a 657 for(Int_t b=0;b<iNbinsPt;b++) {
e8c3ff94 658 Double_t dv2pro = fHistProVPtRP->GetBinContent(b);
659 Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
660 Double_t dErrdifcomb = 0.; //set error to zero
661 Double_t dErr2difcomb = 0.; //set error to zero
662 //calculate error
663 if (dNprime!=0.) {
664 for (Int_t theta=0;theta<iNtheta;theta++) {
665 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
666 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
667 TMath::Cos(dTheta));
668 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
669 TMath::Cos(dTheta));
670 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
671 TMath::BesselJ1(dJ01)))*
672 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
673 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
674 } //loop over theta
675 }
676
677 if (dErr2difcomb!=0.) {
678 dErr2difcomb /= iNtheta;
679 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
680 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
681 }
682 else {dErrdifcomb = 0.;}
683
684 //fill TH1D
685 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
81bbfdbc 686 //calculate integrated flow for RP selection
b7cb54d5 687 if (fHistPtRP){
688 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
81bbfdbc 689 dVRP += dv2pro*dYieldPt;
690 dSum +=dYieldPt;
691 dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
b7cb54d5 692 } else { cout<<"fHistPtRP is NULL"<<endl; }
e8c3ff94 693 } //loop over bins b
81bbfdbc 694
695 if (dSum != 0.) {
696 dVRP /= dSum; //the pt distribution should be normalised
697 dErrV /= (dSum*dSum);
698 dErrV = TMath::Sqrt(dErrV);
699 }
700 cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
701 //cout<<endl;
702 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
703
704
e8c3ff94 705 //v as a function of Pt for POI selection
b7cb54d5 706 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
e8c3ff94 707 Double_t dVPOI = 0.;
81bbfdbc 708 dSum = 0.;
709 dErrV =0.;
e8c3ff94 710
711 for(Int_t b=0;b<iNbinsPt;b++) {
712 Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
713 Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
882ffd6a 714 Double_t dErrdifcomb = 0.; //set error to zero
715 Double_t dErr2difcomb = 0.; //set error to zero
f1d945a1 716 //calculate error
882ffd6a 717 if (dNprime!=0.) {
718 for (Int_t theta=0;theta<iNtheta;theta++) {
719 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
720 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
721 TMath::Cos(dTheta));
722 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
723 TMath::Cos(dTheta));
724 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
725 TMath::BesselJ1(dJ01)))*
726 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
727 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
f1d945a1 728 } //loop over theta
729 }
730
882ffd6a 731 if (dErr2difcomb!=0.) {
732 dErr2difcomb /= iNtheta;
0828acf6 733 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
882ffd6a 734 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
f1d945a1 735 }
882ffd6a 736 else {dErrdifcomb = 0.;}
f1d945a1 737
738 //fill TH1D
e8c3ff94 739 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
740 //calculate integrated flow for POI selection
b7cb54d5 741 if (fHistPtPOI){
742 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
e8c3ff94 743 dVPOI += dv2pro*dYieldPt;
744 dSum +=dYieldPt;
745 dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
b7cb54d5 746 } else { cout<<"fHistPtPOI is NULL"<<endl; }
f1d945a1 747 } //loop over bins b
e8c3ff94 748
749 if (dSum != 0.) {
750 dVPOI /= dSum; //the pt distribution should be normalised
751 dErrV /= (dSum*dSum);
752 dErrV = TMath::Sqrt(dErrV);
753 }
81bbfdbc 754 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
e8c3ff94 755 cout<<endl;
756 cout<<"*************************************"<<endl;
757 cout<<"*************************************"<<endl;
758 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
759
f1d945a1 760 } //secondrun
761
951817d5 762 //cout<<"----LYZ analysis finished....----"<<endl<<endl;
f1d945a1 763
764 return kTRUE;
765}
766
767
768//-----------------------------------------------------------------------
769
0092f3c2 770 Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
f1d945a1 771{
772 // Get event quantities from AliFlowEvent for all particles
773
774 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
775
0092f3c2 776 if (!anEvent){
f1d945a1 777 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
778 return kFALSE;
779 }
780
781 //define variables
882ffd6a 782 TComplex cExpo, cGtheta, cGthetaNew, cZ;
783 Int_t iNtheta = AliFlowLYZConstants::kTheta;
784 Int_t iNbins = AliFlowLYZConstants::kNbins;
f1d945a1 785
786
787 //calculate flow
882ffd6a 788 Double_t dOrder = 2.;
f1d945a1 789
790 //get the Q vector
882ffd6a 791 AliFlowVector vQ = anEvent->GetQ();
f1d945a1 792 //weight by the multiplicity
882ffd6a 793 Double_t dQX = 0;
794 Double_t dQY = 0;
795 if (vQ.GetMult() != 0) {
796 dQX = vQ.X()/vQ.GetMult();
797 dQY = vQ.Y()/vQ.GetMult();
f2ec07bf 798 }
882ffd6a 799 vQ.Set(dQX,dQY);
9d062fe3 800
f1d945a1 801 //for chi calculation:
882ffd6a 802 *fQsum += vQ;
9d062fe3 803 fHistQsumforChi->SetBinContent(1,fQsum->X());
804 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 805 fQ2sum += vQ.Mod2();
9d062fe3 806 fHistQsumforChi->SetBinContent(3,fQ2sum);
f1d945a1 807 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
808
882ffd6a 809 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 810 {
882ffd6a 811 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 812
882ffd6a 813 //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
814 Double_t dQtheta = GetQtheta(vQ, dTheta);
f1d945a1 815
882ffd6a 816 for (Int_t bin=1;bin<=iNbins;bin++)
f1d945a1 817 {
882ffd6a 818 Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED???
819 //if (theta == 0) cerr<<"cerr::dR = "<<dR<<endl;
f1d945a1 820 if (fUseSum)
821 {
822 //calculate the sum generating function
e8c3ff94 823 cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
882ffd6a 824 cGtheta = TComplex::Exp(cExpo);
f1d945a1 825 }
826 else
827 {
828 //calculate the product generating function
e8c3ff94 829 cGtheta = GetGrtheta(anEvent, dR, dTheta);
882ffd6a 830 if (cGtheta.Rho2() > 100.) break;
f1d945a1 831 }
832
882ffd6a 833 fHist1[theta]->Fill(dR,cGtheta); //fill real and imaginary part of cGtheta
f1d945a1 834
835 } //loop over bins
836 } //loop over theta
837
838
839 return kTRUE;
840
841
842}
843
844 //-----------------------------------------------------------------------
0092f3c2 845 Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
f1d945a1 846{
847 //for differential flow
848
849 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
850
0092f3c2 851 if (!anEvent){
f1d945a1 852 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
853 return kFALSE;
854 }
855
856 //define variables
e8c3ff94 857 TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
882ffd6a 858 Double_t dPhi, dEta, dPt;
859 Double_t dR0 = 0.;
e8c3ff94 860 Double_t dCosTermRP = 0.;
861 Double_t dCosTermPOI = 0.;
f1d945a1 862 Double_t m = 1.;
882ffd6a 863 Double_t dOrder = 2.;
864 Int_t iNtheta = AliFlowLYZConstants::kTheta;
f1d945a1 865
882ffd6a 866
f1d945a1 867 //get the Q vector
882ffd6a 868 AliFlowVector vQ = anEvent->GetQ();
f1d945a1 869 //weight by the multiplicity
882ffd6a 870 Double_t dQX = 0.;
871 Double_t dQY = 0.;
872 if (vQ.GetMult() != 0) {
873 dQX = vQ.X()/vQ.GetMult();
874 dQY = vQ.Y()/vQ.GetMult();
f2ec07bf 875 }
9d062fe3 876 vQ.Set(dQX,dQY);
877
f1d945a1 878 //for chi calculation:
882ffd6a 879 *fQsum += vQ;
9d062fe3 880 fHistQsumforChi->SetBinContent(1,fQsum->X());
881 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 882 fQ2sum += vQ.Mod2();
9d062fe3 883 fHistQsumforChi->SetBinContent(3,fQ2sum);
f1d945a1 884
882ffd6a 885 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 886 {
882ffd6a 887 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 888
882ffd6a 889 //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
890 Double_t dQtheta = GetQtheta(vQ, dTheta);
891 //cerr<<"dQtheta for fdenom = "<<dQtheta<<endl;
f1d945a1 892
893 //denominator for differential v
894 if (fHistProR0theta) {
882ffd6a 895 dR0 = fHistProR0theta->GetBinContent(theta+1);
f1d945a1 896 }
897 else {
88e00a8a 898 cout <<"pointer fHistProR0theta does not exist" << endl;
f1d945a1 899 }
882ffd6a 900 //cerr<<"dR0 = "<<dR0 <<endl;
f1d945a1 901
e8c3ff94 902 if (fUseSum) //sum generating function
f1d945a1 903 {
882ffd6a 904 cExpo(0.,dR0*dQtheta);
905 cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
f1d945a1 906 //loop over tracks in event
882ffd6a 907 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
908 for (Int_t i=0;i<iNumberOfTracks;i++) {
909 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i);
910 if (pTrack) {
e8c3ff94 911 dEta = pTrack->Eta();
912 dPt = pTrack->Pt();
913 dPhi = pTrack->Phi();
b7cb54d5 914 if (pTrack->InRPSelection()) { // RP selection
e8c3ff94 915 dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
916 cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
917 if (cNumerRP.Rho()==0) {cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
918 if (fDebug) cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;
919 if (fHist2RP[theta]) {
920 fHist2RP[theta]->Fill(dEta,dPt,cNumerRP); }
f1d945a1 921 }
b7cb54d5 922 if (pTrack->InPOISelection()) { //POI selection
e8c3ff94 923 dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
924 cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
925 if (cNumerPOI.Rho()==0) {cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
926 if (fDebug) cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;
927 if (fHist2POI[theta]) {
928 fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI); }
929 }
3abe32c8 930 // else {
931 // cout << "fHist2 pointer mising" <<endl;
932 // }
e8c3ff94 933 } //if track
f1d945a1 934 else {cerr << "no particle!!!"<<endl;}
935 } //loop over tracks
936
937 } //sum
e8c3ff94 938 else //product generating function
f1d945a1 939 {
882ffd6a 940 cDenom = GetDiffFlow(anEvent, dR0, theta);
f1d945a1 941
942 }//product
943 if (fHistProReDenom && fHistProImDenom) {
882ffd6a 944 fHistProReDenom->Fill(theta,cDenom.Re()); //fill the real part of fDenom
945 fHistProImDenom->Fill(theta,cDenom.Im()); //fill the imaginary part of fDenom
f1d945a1 946 }
947 else {
882ffd6a 948 cout << "Pointers to cDenom mising" << endl;
f1d945a1 949 }
950
951 }//end of loop over theta
952
953 return kTRUE;
954
955
956}
957 //-----------------------------------------------------------------------
882ffd6a 958 Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
f1d945a1 959{
882ffd6a 960 //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
f1d945a1 961 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
962
882ffd6a 963 Double_t dQtheta = 0.;
964 Double_t dOrder = 2.;
f1d945a1 965
882ffd6a 966 dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
f1d945a1 967
882ffd6a 968 return dQtheta;
f1d945a1 969
970}
971
972
973//-----------------------------------------------------------------------
80f72ed6 974TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* const anEvent, Double_t aR, Double_t aTheta)
f1d945a1 975{
976 // Product Generating Function for LeeYangZeros method
977 // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
978
979 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
980
981
882ffd6a 982 TComplex cG = TComplex::One();
983 Double_t dOrder = 2.;
8fbe1a98 984 Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
985
882ffd6a 986 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1 987
882ffd6a 988 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
f1d945a1 989 {
882ffd6a 990 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
991 if (pTrack){
b7cb54d5 992 if (pTrack->InRPSelection()) {
882ffd6a 993 Double_t dPhi = pTrack->Phi();
994 Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
995 TComplex cGi(1., dGIm);
996 cG *= cGi; //product over all tracks
f1d945a1 997 }
998 }
999 else {cerr << "no particle pointer !!!"<<endl;}
1000 }//loop over tracks
1001
882ffd6a 1002 return cG;
f1d945a1 1003
1004}
1005
1006
1007//-----------------------------------------------------------------------
80f72ed6 1008TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* const anEvent, Double_t aR0, Int_t theta)
f1d945a1 1009{
1010 // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1011 // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1012 // Also for v1 mixed harmonics: DF Eq. 5
1013 // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1014
1015 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1016
882ffd6a 1017 TComplex cG = TComplex::One();
1018 TComplex cdGr0(0.,0.);
1019 Double_t dOrder = 2.;
1020 Double_t dWgt = 1.;
f1d945a1 1021
882ffd6a 1022 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1 1023
882ffd6a 1024 Int_t iNtheta = AliFlowLYZConstants::kTheta;
1025 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 1026
e8c3ff94 1027 //for the denominator (use all RP selected particles)
882ffd6a 1028 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
f1d945a1 1029 {
882ffd6a 1030 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1031 if (pTrack){
b7cb54d5 1032 if (pTrack->InRPSelection()) {
882ffd6a 1033 Double_t dPhi = pTrack->Phi();
1034 Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
f1d945a1 1035 //GetGr0theta
882ffd6a 1036 Double_t dGIm = aR0 * dCosTerm;
1037 TComplex cGi(1., dGIm);
e8c3ff94 1038 TComplex cCosTermComplex(1., aR0*dCosTerm);
882ffd6a 1039 cG *= cGi; //product over all tracks
f1d945a1 1040 //GetdGr0theta
882ffd6a 1041 cdGr0 +=(dCosTerm / cCosTermComplex); //sum over all tracks
f1d945a1 1042 }
1043 } //if particle
1044 else {cerr << "no particle!!!"<<endl;}
1045 }//loop over tracks
1046
e8c3ff94 1047 //for the numerator
882ffd6a 1048 for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1 1049 {
882ffd6a 1050 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1051 if (pTrack){
e8c3ff94 1052 Double_t dEta = pTrack->Eta();
1053 Double_t dPt = pTrack->Pt();
1054 Double_t dPhi = pTrack->Phi();
1055 Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1056 TComplex cCosTermComplex(1.,aR0*dCosTerm);
1057 //RP selection
b7cb54d5 1058 if (pTrack->InRPSelection()) {
e8c3ff94 1059 TComplex cNumerRP = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
1060 fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1061 }
1062 //POI selection
b7cb54d5 1063 if (pTrack->InPOISelection()) {
e8c3ff94 1064 TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
1065 fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
f1d945a1 1066 }
1067 } //if particle
1068 else {cerr << "no particle pointer!!!"<<endl;}
1069 }//loop over tracks
1070
e8c3ff94 1071 TComplex cDenom = cG*cdGr0;
882ffd6a 1072 return cDenom;
f1d945a1 1073
1074}
1075
1076//-----------------------------------------------------------------------
1077