]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowAnalysisWithLeeYangZeros.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moved old reacton plane code to /oldcode
[u/mrichter/AliRoot.git] / PWG2 / FLOW / 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
16/*
17$Log$
18*/
19
20#include "Riostream.h" //needed as include
21#include "AliFlowCommonConstants.h" //needed as include
22#include "AliFlowLYZConstants.h" //needed as include
23#include "AliFlowAnalysisWithLeeYangZeros.h"
24#include "AliFlowLYZHist1.h" //needed as include
25#include "AliFlowLYZHist2.h" //needed as include
26#include "AliFlowCommonHist.h" //needed as include
27#include "AliFlowCommonHistResults.h" //needed as include
28#include "AliFlowEventSimple.h" //needed as include
29#include "AliFlowTrackSimple.h" //needed as include
30
31class AliFlowVector;
32
33#include "TMath.h" //needed as include
34#include "TFile.h" //needed as include
35
36class TComplex;
37class TProfile;
38class TH1F;
39class TH1D;
40
41
42//Description: Maker to analyze Flow using the LeeYangZeros method
43// Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
44// Practical Guide (PG): J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
45// Adapted from StFlowLeeYangZerosMaker.cxx
46// by Markus Oldenberg and Art Poskanzer, LBNL
47// with advice from Jean-Yves Ollitrault and Nicolas Borghini
48//
49//Author: Naomi van der Kolk (kolk@nikhef.nl)
50
51
52
53ClassImp(AliFlowAnalysisWithLeeYangZeros)
54
55 //-----------------------------------------------------------------------
56
57 AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
58 fQ2sum(0),
59 fQtheta(0),
60 fEventNumber(0),
61 fMult(0),
62 fNbins(0),
63 fTheta(0),
64 fEvent(0),
65 fTrack(0),
66 fFirstRun(kTRUE),
67 fUseSum(kTRUE),
68 fDoubleLoop(kFALSE),
69 fDebug(kFALSE),
70 fHistFileName(0),
71 fHistFile(0),
72 fSummaryFile(0),
73 firstRunFileName(0),
74 firstRunFile(0),
75 fHistProVtheta(0),
76 fHistProVeta(0),
77 fHistProVPt(0),
78 fHistProR0theta(0),
79 fHistProReDenom(0),
80 fHistProImDenom(0),
81 fHistProReDtheta(0),
82 fHistProImDtheta(0),
83 fCommonHists(0),
84 fCommonHistsRes(0)
85
86{
87 //default constructor
88 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
89
90 // output file (histograms)
91 TString fHistFileName = "outputFromLYZAnalysis.root" ;
92
93 for(Int_t i = 0;i<5;i++)
94 {
95 fHist1[i]=0;
96 fHist2[i]=0;
97 }
98
99 fQ.Set(0.,0.);
100 fQsum.Set(0.,0.);
101
102}
103
104
105//-----------------------------------------------------------------------
106
107
108 AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
109 {
110 //default destructor
111 if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
112 delete fHistFile;
113 }
114
115 //-----------------------------------------------------------------------
116
117Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
118{
119 //init method
120 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
121
122 // Open output files (->plots)
123 fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
124
125
126 // Book histograms
127 Int_t fNtheta = AliFlowLYZConstants::kTheta;
128 Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
129 Int_t fNbinsEta = AliFlowCommonConstants::GetNbinsEta();
130
131 Double_t fPtMin = AliFlowCommonConstants::GetPtMin();
132 Double_t fPtMax = AliFlowCommonConstants::GetPtMax();
133 Double_t fEtaMin = AliFlowCommonConstants::GetEtaMin();
134 Double_t fEtaMax = AliFlowCommonConstants::GetEtaMax();
135
136
137 //for control histograms
138 fCommonHists = new AliFlowCommonHist("LYZ");
139 fCommonHistsRes = new AliFlowCommonHistResults("LYZ");
140
141 //for first loop over events
142 if (fFirstRun){
143 fHistProR0theta = new TProfile("First_FlowPro_r0theta_LYZ","First_FlowPro_r0theta_LYZ",fNtheta,-0.5,fNtheta-0.5);
144 fHistProR0theta->SetXTitle("#theta");
145 fHistProR0theta->SetYTitle("r_{0}^{#theta}");
146
147 fHistProVtheta = new TProfile("First_FlowPro_Vtheta_LYZ","First_FlowPro_Vtheta_LYZ",fNtheta,-0.5,fNtheta-0.5);
148 fHistProVtheta->SetXTitle("#theta");
149 fHistProVtheta->SetYTitle("V_{n}^{#theta}");
150
151 //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
152 for (Int_t theta=0;theta<fNtheta;theta++) {
153 fHist1[theta]=new AliFlowLYZHist1(theta);
154 }
155
156 }
157 //for second loop over events
158 else {
159 fHistProReDenom = new TProfile("Second_FlowPro_ReDenom_LYZ","Second_FlowPro_ReDenom_LYZ" , fNtheta, -0.5, fNtheta-0.5);
160 fHistProReDenom->SetXTitle("#theta");
161 fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
162
163 fHistProImDenom = new TProfile("Second_FlowPro_ImDenom_LYZ","Second_FlowPro_ImDenom_LYZ" , fNtheta, -0.5, fNtheta-0.5);
164 fHistProImDenom->SetXTitle("#theta");
165 fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
166
167 fHistProVeta = new TProfile("Second_FlowPro_Veta_LYZ","Second_FlowPro_Veta_LYZ",fNbinsEta,fEtaMin,fEtaMax);
168 fHistProVeta->SetXTitle("rapidity");
169 fHistProVeta->SetYTitle("v (%)");
170
171 fHistProVPt = new TProfile("Second_FlowPro_VPt_LYZ","Second_FlowPro_VPt_LYZ",fNbinsPt,fPtMin,fPtMax);
172 fHistProVPt->SetXTitle("Pt");
173 fHistProVPt->SetYTitle("v (%)");
174
175 fHistProReDtheta = new TProfile("Second_FlowPro_ReDtheta_LYZ","Second_FlowPro_ReDtheta_LYZ",fNtheta, -0.5, fNtheta-0.5);
176 fHistProReDtheta->SetXTitle("#theta");
177 fHistProReDtheta->SetYTitle("Re(D^{#theta})");
178
179 fHistProImDtheta = new TProfile("Second_FlowPro_ImDtheta_LYZ","Second_FlowPro_ImDtheta_LYZ",fNtheta, -0.5, fNtheta-0.5);
180 fHistProImDtheta->SetXTitle("#theta");
181 fHistProImDtheta->SetYTitle("Im(D^{#theta})");
182
183 //class AliFlowLYZHist2 defines the histograms:
184 for (Int_t theta=0;theta<fNtheta;theta++) {
185 fHist2[theta]=new AliFlowLYZHist2(theta);
186 }
187
188 //read hists from first run file
189 //firstRunFile = new TFile("fof_flowLYZAnal_firstrun.root","READ"); //default is read
190 if (firstRunFile->IsZombie()){ //check if file exists
191 cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
192 exit(-1);
193 } else if (firstRunFile->IsOpen()){
194 cout<<"----firstRunFile is open----"<<endl<<endl;
195 fHistProR0theta = (TProfile*)firstRunFile->Get("First_FlowPro_r0theta_LYZ");
196 }
197 }
198
199
200 if (fDebug) cout<<"****Histograms initialised****"<<endl;
201
202 fEventNumber = 0; //set event counter to zero
203
204 return kTRUE;
205}
206
207 //-----------------------------------------------------------------------
208
209Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* fEvent)
210{
211 //make method
212 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
213
214 //get tracks from event
215 if (fEvent) {
216 if (fFirstRun){
217 fCommonHists->FillControlHistograms(fEvent);
218 FillFromFlowEvent(fEvent);
219 }
220 else {
221 fCommonHists->FillControlHistograms(fEvent);
222 SecondFillFromFlowEvent(fEvent);
223 }
224 }
225
226 else {
227 cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
228 return kFALSE;
229 }
230 cout<<"^^^^read event "<<fEventNumber<<endl;
231 fEventNumber++;
232
233
234 return kTRUE;
235}
236
237
238 //-----------------------------------------------------------------------
239 Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
240{
241 //finish method
242 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
243
244 //define variables for both runs
245 Double_t fJ01 = 2.405;
246 Int_t fNtheta = AliFlowLYZConstants::kTheta;
247
248 if (fFirstRun){
249 Double_t fR0 = 0;
250 Double_t fVtheta = 0;
251 Double_t fv = 0;
252 Double_t fV = 0;
253 for (Int_t theta=0;theta<fNtheta;theta++)
254 {
255 //get the first minimum r0
256 fHist1[theta]->FillGtheta();
257 fR0 = fHist1[theta]->GetR0();
258
259 //calculate integrated flow
260 if (fR0!=0) fVtheta = fJ01/fR0;
261 //for estimating systematic error resulting from r0
262 Double_t binsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins);
263 Double_t fVplus = fJ01/(fR0+binsize);
264 Double_t fVmin = fJ01/(fR0-binsize);
265 fv = fVtheta;
266 Double_t fvplus = fVplus;
267 Double_t fvmin = fVmin;
268 cout<<"fv = "<<fv<<" and fvplus = "<<fvplus<< " and fvmin = "<<fvmin<<endl;
269
270 //fill the histograms
271 fHistProR0theta->Fill(theta,fR0);
272 fHistProVtheta->Fill(theta,fVtheta);
273 //get average value of fVtheta = fV
274 fV += fVtheta;
275 } //end of loop over theta
276
277 //get average value of fVtheta = fV
278 fV /=fNtheta;
279
280 //sigma2 and chi
281 Double_t fSigma2 = 0;
282 Double_t fChi= 0;
283 if (fEventNumber!=0) {
284 fQsum /= fEventNumber;
285 fQ2sum /= fEventNumber;
286 fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.); //BP eq. 62
287 if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
288 else fChi = -1.;
289 fCommonHistsRes->FillChi(fChi);
290 cout<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
291 }
292
293 // recalculate statistical errors on integrated flow
294 //combining 5 theta angles to 1 relative error BP eq. 89
295 Double_t fRelErr2comb = 0.;
296 Int_t nEvts = fEventNumber;
297 for (Int_t theta=0;theta<fNtheta;theta++){
298 Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
299 Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
300 TMath::Cos(fTheta));
301 Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
302 TMath::Cos(fTheta));
303 fRelErr2comb += (1/(2*nEvts*(fJ01*fJ01)*TMath::BesselJ1(fJ01)*
304 TMath::BesselJ1(fJ01)))*
305 (fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2)) +
306 fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2)));
307 }
308 fRelErr2comb /= fNtheta;
309 Double_t fRelErrcomb = TMath::Sqrt(fRelErr2comb);
310 cout<<"fRelErrcomb = "<<fRelErrcomb<<endl;
311
312 //copy content of profile into TH1D and add error
313 Double_t fv2pro = fV; //in the case that fv is equal to fV
314 Double_t fv2Err = fv2pro*fRelErrcomb ;
315 cout<<"fv2pro +- fv2Err = "<<fv2pro<<" +- "<<fv2Err<<endl;
316 fCommonHistsRes->FillIntegratedFlow(fv2pro, fv2Err);
317
318
319 if (fDebug) cout<<"****histograms filled****"<<endl;
320
321 //save histograms in file //temp for testing selector
322 fHistFile->cd();
323 fHistFile->Write();
324
325 return kTRUE;
326 fEventNumber =0; //set to zero for second round over events
327 } //firstrun
328
329
330 else { //second run
331
332 //calculate differential flow
333 //declare variables
334 TComplex fDenom, fNumer, fDtheta;
335 Int_t m = 1;
336 TComplex i = TComplex::I();
337 Double_t fBesselRatio[3] = {1., 1.202, 2.69};
338 Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
339 Int_t fNbinsEta = AliFlowCommonConstants::GetNbinsEta();
340
341 Double_t fEta, fPt, fReRatio, fVeta, fVPt;
342
343
344 Double_t fR0 = 0.;
345 Double_t fVtheta = 0.;
346 Double_t fV = 0.;
347 Double_t fReDenom = 0.;
348 Double_t fImDenom = 0.;
349 for (Int_t theta=0;theta<fNtheta;theta++) { //loop over theta
350 if (fHistProR0theta) {
351 fR0 = fHistProR0theta->GetBinContent(theta+1);
352 if (fDebug) cerr<<"fR0 = "<<fR0<<endl;
353 if (fR0!=0) fVtheta = fJ01/fR0;
354 fV += fVtheta;
355 // BP Eq. 9 -> Vn^theta = j01/r0^theta
356 if (fHistProReDenom && fHistProImDenom) {
357 fReDenom = fHistProReDenom->GetBinContent(theta+1);
358 fImDenom = fHistProImDenom->GetBinContent(theta+1);
359 }
360 else {
361 cout << "Hist pointer fDenom in file does not exist" <<endl;
362 }
363
364 }
365 else {
366 cout << "Hist pointer R0theta in file does not exist" <<endl;
367 }
368 //} //loop over theta
369
370 fDenom(fReDenom,fImDenom);
371
372
373 //for new method and use by others (only with the sum generating function):
374 if (fUseSum) {
375 fR0 = fHistProR0theta->GetBinContent(theta+1);
376 fDtheta = fR0*fDenom/fJ01;
377 fHistProReDtheta->Fill(theta,fDtheta.Re());
378 fHistProImDtheta->Fill(theta,fDtheta.Im());
379 }
380
381 fDenom *= TComplex::Power(i, m-1);
382 //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
383
384 //v as a function of eta
385 for (Int_t be=1;be<=fNbinsEta;be++) {
386 fEta = fHist2[theta]->GetBinCenter(be);
387 fNumer = fHist2[theta]->GetfNumer(be);
388 if (fNumer.Rho()==0) {
389 if (fDebug) cerr<<"WARNING: modulus of fNumer is zero in Finish()"<<endl;
390 fReRatio = 0;
391 }
392 else if (fDenom.Rho()==0) {
393 if (fDebug) cerr<<"WARNING: modulus of fDenom is zero"<<endl;
394 fReRatio = 0;
395 }
396 else {
397 //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
398 fReRatio = (fNumer/fDenom).Re();
399 }
400
401 fVeta = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
402 //if ( j==1 && theta==0) cerr<<"eta = "<<fEta<<" cerr::fReRatio for eta = "<<fReRatio<<" cerr::fVeta for eta = "<<fVeta<<endl;
403
404 fHistProVeta->Fill(fEta,fVeta);
405 } //loop over bins be
406
407 //v as a function of Pt
408 Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
409 for (Int_t bp=1;bp<=fNbinsPt;bp++) {
410 fPt = fHist2[theta]->GetBinCenterPt(bp);
411 fNumer = fHist2[theta]->GetfNumerPt(bp);
412 if (fNumer.Rho()==0) {
413 if (fDebug) cerr<<"modulus of fNumer is zero"<<endl;
414 fReRatio = 0;
415 }
416 else if (fDenom.Rho()==0) {
417 if (fDebug) cerr<<"modulus of fDenom is zero"<<endl;
418 fReRatio = 0;
419 }
420 else {
421 //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
422 fReRatio = (fNumer/fDenom).Re();
423 }
424
425 fVPt = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
426
427 fHistProVPt->Fill(fPt,fVPt);
428 } //loop over bins bp
429
430 }//end of loop over theta
431
432 //calculate the average of fVtheta = fV
433 fV /= fNtheta;
434
435 //sigma2 and chi (for statistical error calculations)
436 Double_t fSigma2 = 0;
437 Double_t fChi= 0;
438 if (fEventNumber!=0) {
439 fQsum /= fEventNumber;
440 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
441 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
442 fQ2sum /= fEventNumber;
443 //cout<<"fQ2sum = "<<fQ2sum<<endl;
444 fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.); //BP eq. 62
445 if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
446 else fChi = -1.;
447 fCommonHistsRes->FillChi(fChi);
448 cout<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
449 }
450
451 //copy content of profile into TH1D and add error
452 for(Int_t b=0;b<fNbinsPt;b++) {
453 Double_t fv2pro = fHistProVPt->GetBinContent(b);
454 Double_t fNprime = fCommonHists->GetEntriesInPtBin(b);
455 Double_t fErrdifcomb = 0.; //set error to zero
456 Double_t fErr2difcomb = 0.; //set error to zero
457 //calculate error
458 if (fNprime!=0.) {
459 for (Int_t theta=0;theta<fNtheta;theta++) {
460 Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
461 Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
462 TMath::Cos(fTheta));
463 Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
464 TMath::Cos(fTheta));
465 fErr2difcomb += (TMath::Cos(fTheta)/(4*fNprime*TMath::BesselJ1(fJ01)*
466 TMath::BesselJ1(fJ01)))*
467 ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) -
468 (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
469 } //loop over theta
470 }
471
472 if (fErr2difcomb!=0.) {
473 fErr2difcomb /= fNtheta;
474 fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
475 //cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
476 }
477 else {fErrdifcomb = 0.;}
478
479 //fill TH1D
480 fCommonHistsRes->FillDifferentialFlow(b, fv2pro, fErrdifcomb);
481 } //loop over bins b
482
483
484 //save histograms in file
485 fHistFile->cd();
486 fHistFile->Write();
487 fHistFile->Close();
488 //Note that when the file is closed, all histograms and Trees in memory associated with this file are deleted
489 if (fDebug) cout<<"****Histograms saved and fHistFile closed, all histograms deleted****"<<endl;
490
491 //close the first run file
492 firstRunFile->Close();
493
494
495 } //secondrun
496
497 cout<<"----LYZ analysis finished....----"<<endl<<endl;
498
499 return kTRUE;
500}
501
502
503//-----------------------------------------------------------------------
504
505 Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* fEvent)
506{
507 // Get event quantities from AliFlowEvent for all particles
508
509 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
510
511 if (!fEvent){
512 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
513 return kFALSE;
514 }
515
516 //define variables
517 TComplex fExpo, fGtheta, fGthetaNew, fZ;
518 Int_t fNtheta = AliFlowLYZConstants::kTheta;
519 Int_t fNbins = AliFlowLYZConstants::kNbins;
520
521
522 //calculate flow
523 Double_t fOrder = 2.;
524
525 //get the Q vector
526 fQ = fEvent->GetQ();
527 //weight by the multiplicity
528 Double_t fQX = fQ.X()/fQ.GetMult();
529 Double_t fQY = fQ.Y()/fQ.GetMult();
530 fQ.Set(fQX,fQY);
531 //for chi calculation:
532 fQsum += fQ;
533 fQ2sum += fQ.Mod2();
534 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
535
536 for (Int_t theta=0;theta<fNtheta;theta++)
537 {
538 fTheta = ((double)theta/fNtheta)*TMath::Pi()/fOrder;
539
540 //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
541 fQtheta = GetQtheta(fQ, fTheta);
542
543 for (Int_t bin=1;bin<=fNbins;bin++)
544 {
545 Double_t fR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED???
546 //if (theta == 0) cerr<<"cerr::fR = "<<fR<<endl;
547 if (fUseSum)
548 {
549 //calculate the sum generating function
550 fExpo(0.,fR*fQtheta); //Re=0 ; Im=fR*fQtheta
551 fGtheta = TComplex::Exp(fExpo);
552 }
553 else
554 {
555 //calculate the product generating function
556 fGtheta = GetGrtheta(fEvent, fR, fTheta); //make this function
557 if (fGtheta.Rho2() > 100.) break;
558 }
559
560 fHist1[theta]->Fill(fR,fGtheta); //fill real and imaginary part of fGtheta
561
562 } //loop over bins
563 } //loop over theta
564
565
566 return kTRUE;
567
568
569}
570
571 //-----------------------------------------------------------------------
572 Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* fEvent)
573{
574 //for differential flow
575
576 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
577
578 if (!fEvent){
579 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
580 return kFALSE;
581 }
582
583 //define variables
584 //TVector2 fQ;
585 TComplex fExpo, fDenom, fNumer,fCosTermComplex;
586 Double_t fOrder, fPhi, fEta, fPt;
587 Double_t fR0 = 0.;
588 Double_t fCosTerm;
589 Double_t m = 1.;
590 Int_t fNtheta = AliFlowLYZConstants::kTheta;
591
592 //calculate flow
593 fOrder = 2.;
594
595 //get the Q vector
596 fQ = fEvent->GetQ();
597 //weight by the multiplicity
598 Double_t fQX = fQ.X()/fQ.GetMult();
599 Double_t fQY = fQ.Y()/fQ.GetMult();
600 fQ.Set(fQX,fQY);
601 //for chi calculation:
602 fQsum += fQ;
603 fQ2sum += fQ.Mod2();
604
605 for (Int_t theta=0;theta<fNtheta;theta++)
606 {
607 fTheta = ((double)theta/fNtheta)*TMath::Pi()/fOrder;
608
609 //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
610 fQtheta = GetQtheta(fQ, fTheta);
611 //cerr<<"fQtheta for fdenom = "<<fQtheta<<endl;
612
613 //denominator for differential v
614 if (fHistProR0theta) {
615 fR0 = fHistProR0theta->GetBinContent(theta+1);
616 }
617 else {
618 cout <<"pointer fHistProR0Theta does not exist" << endl;
619 }
620 //cerr<<"fR0 = "<<fR0 <<endl;
621
622 if (fUseSum) //sum generating function
623 {
624 fExpo(0.,fR0*fQtheta);
625 fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
626 //loop over tracks in event
627 Int_t fNumberOfTracks = fEvent->NumberOfTracks();
628 for (Int_t i=0;i<fNumberOfTracks;i++) {
629 fTrack = fEvent->GetTrack(i);
630 if (fTrack) {
631 if (fTrack->UseForDifferentialFlow()) {
632 fEta = fTrack->Eta();
633 fPt = fTrack->Pt();
634 fPhi = fTrack->Phi();
635 fCosTerm = cos(m*fOrder*(fPhi-fTheta));
636 //cerr<<"fCosTerm = "<<fCosTerm <<endl;
637 fNumer = fCosTerm*(TComplex::Exp(fExpo));
638 if (fNumer.Rho()==0) {cerr<<"WARNING: modulus of fNumer is zero in SecondFillFromFlowEvent"<<endl;}
639 if (fDebug) cerr<<"modulus of fNumer is "<<fNumer.Rho()<<endl;
640 if (fHist2[theta]) {
641 fHist2[theta]->Fill(fEta,fPt,fNumer);
642 }
643 else {
644 cout << "fHist2 pointer mising" <<endl;
645 }
646 }
647 } //if particle
648 else {cerr << "no particle!!!"<<endl;}
649 } //loop over tracks
650
651 } //sum
652 else //product generating function
653 {
654 fDenom = GetDiffFlow(fEvent, fR0, theta);
655
656 }//product
657 if (fHistProReDenom && fHistProImDenom) {
658 fHistProReDenom->Fill(theta,fDenom.Re()); //fill the real part of fDenom
659 fHistProImDenom->Fill(theta,fDenom.Im()); //fill the imaginary part of fDenom
660 }
661 else {
662 cout << "Pointers to fDenom mising" << endl;
663 }
664
665 }//end of loop over theta
666
667 return kTRUE;
668
669
670}
671 //-----------------------------------------------------------------------
672 Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(TVector2 fQ, Double_t fTheta)
673{
674 //calculate Qtheta. Qtheta is the sum over all particles of cos(fOrder*(fPhi-fTheta)) BP eq. 3
675 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
676
677 Double_t fQtheta = 0.;
678 Double_t fOrder = 2.;
679
680 fQtheta = fQ.X()*cos(fOrder*fTheta)+fQ.Y()*sin(fOrder*fTheta);
681
682 return fQtheta;
683
684}
685
686
687//-----------------------------------------------------------------------
688TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* fEvent, Double_t fR, Double_t fTheta)
689{
690 // Product Generating Function for LeeYangZeros method
691 // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
692
693 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
694
695
696 TComplex fG = TComplex::One();
697 Double_t fOrder = 2.;
698 Double_t fWgt = 1.;
699
700 Int_t fNumberOfTracks = fEvent->NumberOfTracks();
701
702 for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
703 {
704 fTrack = fEvent->GetTrack(i) ;
705 if (fTrack){
706 if (fTrack->UseForIntegratedFlow()) {
707 Double_t fPhi = fTrack->Phi();
708 Double_t fGIm = fR * fWgt*cos(fOrder*(fPhi - fTheta));
709 TComplex fGi(1., fGIm);
710 fG *= fGi; //product over all tracks
711 }
712 }
713 else {cerr << "no particle pointer !!!"<<endl;}
714 }//loop over tracks
715
716 return fG;
717
718}
719
720
721//-----------------------------------------------------------------------
722TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* fEvent, Double_t fR0, Int_t theta)
723{
724 // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
725 // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
726 // Also for v1 mixed harmonics: DF Eq. 5
727 // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
728
729 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
730
731 TComplex fG = TComplex::One();
732 TComplex fdGr0(0.,0.);
733 Double_t fOrder = 2.;
734 Double_t fWgt = 1.;
735
736 Int_t fNumberOfTracks = fEvent->NumberOfTracks();
737
738 Int_t fNtheta = AliFlowLYZConstants::kTheta;
739 Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi()/fOrder;
740
741 for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
742 {
743 fTrack = fEvent->GetTrack(i) ;
744 if (fTrack){
745 if (fTrack->UseForDifferentialFlow()) {
746 Double_t fPhi = fTrack->Phi();
747 Double_t fCosTerm = fWgt*cos(fOrder*(fPhi - fTheta));
748 //GetGr0theta
749 Double_t fGIm = fR0 * fCosTerm;
750 TComplex fGi(1., fGIm);
751 fG *= fGi; //product over all tracks
752 //GetdGr0theta
753 TComplex fCosTermComplex(1., fR0*fCosTerm);
754 fdGr0 +=(fCosTerm / fCosTermComplex); //sum over all tracks
755 }
756 } //if particle
757 else {cerr << "no particle!!!"<<endl;}
758 }//loop over tracks
759
760
761 for (Int_t i=0;i<fNumberOfTracks;i++)
762 {
763 fTrack = fEvent->GetTrack(i) ;
764 if (fTrack){
765 if (fTrack->UseForDifferentialFlow()) {
766 Double_t fEta = fTrack->Eta();
767 Double_t fPt = fTrack->Pt();
768 Double_t fPhi = fTrack->Phi();
769 Double_t fCosTerm = cos(fOrder*(fPhi-fTheta));
770 TComplex fCosTermComplex(1.,fR0*fCosTerm);
771 TComplex fNumer = fG*fCosTerm/fCosTermComplex; //PG Eq. 9
772 fHist2[theta]->Fill(fEta,fPt,fNumer);
773 }
774 } //if particle
775 else {cerr << "no particle pointer!!!"<<endl;}
776 }//loop over tracks
777
778 TComplex fDenom = fG*fdGr0;
779 return fDenom;
780
781}
782
783//-----------------------------------------------------------------------
784