]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowLeeYangZerosMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
gain set to 1 for all ch
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowLeeYangZerosMaker.cxx
CommitLineData
f456b167 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"
21#include "AliFlowLeeYangZerosMaker.h"
22#include "AliFlowEvent.h"
23#include "AliFlowSelection.h"
24#include "AliFlowConstants.h" //??
25#include "AliFlowLYZHist1.h"
26#include "AliFlowLYZHist2.h"
27#include "AliFlowLYZConstants.h" //??
f456b167 28
29#include "TMath.h"
30#include "TComplex.h"
31#include "TProfile.h"
32#include "TObjArray.h"
33#include "TFile.h"
34#include "TVector.h"
35#include "TVector2.h"
36#include "TGraphErrors.h"
37#include "TCanvas.h"
38
39class TTree;
40class TH1F;
41class TH1D;
42class TVector3;
43class TProfile2D;
44class TObject;
45
46//class Riostream; //does not compile
47//class TMath; //does not compile
48//class TVector; //does not compile
49
50//Description: Maker to analyze Flow using the LeeYangZeros method
51// Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
52// Practical Guide (PG): J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
53// Adapted from StFlowLeeYangZerosMaker.cxx
54// by Markus Oldenberg and Art Poskanzer, LBNL
55// with advice from Jean-Yves Ollitrault and Nicolas Borghini
56//
57//Author: Naomi van der Kolk (kolk@nikhef.nl)
58
59
60
61ClassImp(AliFlowLeeYangZerosMaker)
62
63 //-----------------------------------------------------------------------
64
65 AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker():
66 fFirstRun(kTRUE),
67 fUseSum(kTRUE),
68 fDebug(kFALSE),
69 fHistFileName(0),
70 fHistFile(0),
71 fSummaryFile(0),
72 firstRunFile(0)
73
74{
75 //default constructor
76 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker default constructor****"<<endl;
77
78 fFlowSelect = new AliFlowSelection();
79 if (fDebug) { cerr<<"****fFlowSelect in constructor AliFlowLeeYangZerosMaker ("<<fFlowSelect<<")****"<<endl;}
80 // output file (histograms)
81 TString fHistFileName = "flowLYZAnalysPlot.root" ;
82}
83
84
85AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker(const AliFlowSelection* flowSelect):
86 fFirstRun(kTRUE),
87 fUseSum(kTRUE),
88 fDebug(kFALSE),
89 fHistFileName(0),
90 fHistFile(0),
91 fSummaryFile(0),
92 firstRunFile(0)
93{
94 //custum constructor
95 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker custum constructor****"<<endl;
96
97 if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect); }
98 else {
99 fFlowSelect = new AliFlowSelection() ;
100 if (fDebug) cerr<<"****fFlowSelect in constructor AliFlowLeeYangZerosMaker ("<<fFlowSelect<<")****"<<endl;
101 }
102 // output file (histograms)
103 TString fHistFileName = "flowLYZAnalysPlot.root" ;
104}
105
106 //-----------------------------------------------------------------------
107
108
109 AliFlowLeeYangZerosMaker::~AliFlowLeeYangZerosMaker()
110 {
111 //default destructor
112 if (fDebug) cout<<"****~AliFlowLeeYangZerosMaker****"<<endl;
113 delete fHistFile;
114 }
115
116 //-----------------------------------------------------------------------
117
118Bool_t AliFlowLeeYangZerosMaker::Init()
119{
120 //init method
121 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Init()****"<<endl;
122
123 // Open output files (->plots)
124 fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
125 //fHistFile->cd() ; //all histograms will be saved in this file
126
127 //for each harmonic ???
128 fQsum.Set(0.,0.);
129 fQ2sum = 0.;
130
131 // Book histograms
132 Int_t fNtheta = AliFlowLYZConstants::kTheta;
133
134
135 //for control histograms
136 fHistOrigMult = new TH1F("Control_FlowLYZ_OrigMult", "Control_FlowLYZ_OrigMult",1000, 0., 10000.);
137 fHistOrigMult->SetXTitle("Original Multiplicity");
138 fHistOrigMult->SetYTitle("Counts");
139
140 fHistMult = new TH1F("Control_FlowLYZ_Mult", "Control_FlowLYZ_Mult",1000, 0., 10000.);
141 fHistMult->SetXTitle("Multiplicity from selection");
142 fHistMult->SetYTitle("Counts");
143
144 fHistQ = new TH1F("Control_FlowLYZ_Q","Control_FlowLYZ_Q",500, 0., 10.);
145 fHistQ->SetXTitle("Qvector");
146 fHistQ->SetYTitle("Counts");
147
148 fHistPt = new TH1F("Control_FlowLYZ_Pt","Control_FlowLYZ_Pt",200, 0., 10.);
149 fHistPt->SetXTitle("Pt (GeV/c)");
150 fHistPt->SetYTitle("Counts");
151
152 fHistPhi = new TH1F("Control_FlowLYZ_Phi","Control_FlowLYZ_Phi",70, 0., 7.);
153 fHistPhi->SetXTitle("Phi");
154 fHistPhi->SetYTitle("Counts");
155
156 fHistEta = new TH1F("Control_FlowLYZ_Eta","Control_FlowLYZ_Eta",40, 0., 2.);
157 fHistEta->SetXTitle("Eta");
158 fHistEta->SetYTitle("Counts");
159
160 fHistQtheta = new TH1F("Control_FlowLYZ_Qtheta","Control_FlowLYZ_Qtheta",50,-1000.,1000.);
161 fHistQtheta->SetXTitle("Qtheta");
162 fHistQtheta->SetYTitle("Counts");
163
164 if (fFirstRun){
165 //for first loop over events
166 fHistProR0thetaHar1 = new TProfile("First_FlowProLYZ_r0theta_Har1","First_FlowProLYZ_r0theta_Har1",fNtheta,-0.5,fNtheta-0.5);
167 fHistProR0thetaHar1->SetXTitle("#theta");
168 fHistProR0thetaHar1->SetYTitle("r_{0}^{#theta}");
169
170 fHistProR0thetaHar2 = new TProfile("First_FlowProLYZ_r0theta_Har2","First_FlowProLYZ_r0theta_Har2",fNtheta,-0.5,fNtheta-0.5);
171 fHistProR0thetaHar2->SetXTitle("#theta");
172 fHistProR0thetaHar2->SetYTitle("r_{0}^{#theta}");
173
174 fHistProVthetaHar1 = new TProfile("First_FlowProLYZ_Vtheta_Har1","First_FlowProLYZ_Vtheta_Har1",fNtheta,-0.5,fNtheta-0.5);
175 fHistProVthetaHar1->SetXTitle("#theta");
176 fHistProVthetaHar1->SetYTitle("V_{n}^{#theta}");
177
178 fHistProVthetaHar2 = new TProfile("First_FlowProLYZ_Vtheta_Har2","First_FlowProLYZ_Vtheta_Har2",fNtheta,-0.5,fNtheta-0.5);
179 fHistProVthetaHar2->SetXTitle("#theta");
180 fHistProVthetaHar2->SetYTitle("V_{n}^{#theta}");
181
182 fHistProVR0 = new TProfile("First_FlowProLYZ_vR0","First_FlowProLYZ_vR0",2,0.5,2.5,-100.,100.);
183 fHistProVR0->SetXTitle("Harmonic");
184 fHistProVR0->SetYTitle("v integrated from r0 (%)");
185
186 fHistVR0 = new TH1D("First_FlowLYZ_vR0","First_FlowLYZ_vR0",2,0.5,2.5);
187 fHistVR0->SetXTitle("Harmonic");
188 fHistVR0->SetYTitle("v integrated from r0 (%)");
189
190 fHistProV = new TProfile("First_FlowProLYZ_V","First_FlowProLYZ_V",2,0.5,2.5,-1000.,1000.);
191 fHistProV->SetXTitle("Harmonic");
192 fHistProV->SetYTitle("v integrated");
193
194 //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
195 for (Int_t j=0;j<AliFlowConstants::kHars;j++)
196 {
197 for (Int_t theta=0;theta<fNtheta;theta++)
198 {
199 fHist1[j][theta]=new AliFlowLYZHist1(theta,j+1);
200 }
201 }
202 }
203 else {
204 //for second loop over events
205 fHistProReDenomHar1 = new TProfile("Second_FlowProLYZ_ReDenom_Har1","Second_FlowProLYZ_ReDenom_Har1" , fNtheta, -0.5, fNtheta-0.5);
206 fHistProReDenomHar1->SetXTitle("#theta");
207 fHistProReDenomHar1->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
208
209 fHistProReDenomHar2 = new TProfile("Second_FlowProLYZ_ReDenom_Har2","Second_FlowProLYZ_ReDenom_Har2" , fNtheta, -0.5, fNtheta-0.5);
210 fHistProReDenomHar2->SetXTitle("#theta");
211 fHistProReDenomHar2->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
212
213 fHistProImDenomHar1 = new TProfile("Second_FlowProLYZ_ImDenom_Har1","Second_FlowProLYZ_ImDenom_Har1" , fNtheta, -0.5, fNtheta-0.5);
214 fHistProImDenomHar1->SetXTitle("#theta");
215 fHistProImDenomHar1->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
216
217 fHistProImDenomHar2 = new TProfile("Second_FlowProLYZ_ImDenom_Har2","Second_FlowProLYZ_ImDenom_Har2" , fNtheta, -0.5, fNtheta-0.5);
218 fHistProImDenomHar2->SetXTitle("#theta");
219 fHistProImDenomHar2->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
220
221 fHistProVetaHar1 = new TProfile("Second_FlowProLYZ_Veta_Har1","Second_FlowProLYZ_Veta_Har1",40,-2.,2.);
222 fHistProVetaHar1->SetXTitle("rapidity");
223 fHistProVetaHar1->SetYTitle("v (%)");
224
225 fHistProVetaHar2 = new TProfile("Second_FlowProLYZ_Veta_Har2","Second_FlowProLYZ_Veta_Har2",40,-2.,2.);
226 fHistProVetaHar2->SetXTitle("rapidity");
227 fHistProVetaHar2->SetYTitle("v (%)");
228
229 fHistProVPtHar1 = new TProfile("Second_FlowProLYZ_VPt_Har1","Second_FlowProLYZ_VPt_Har1",100,0.,10.);
230 fHistProVPtHar1->SetXTitle("Pt");
231 fHistProVPtHar1->SetYTitle("v (%)");
232
233 fHistProVPtHar2 = new TProfile("Second_FlowProLYZ_VPt_Har2","Second_FlowProLYZ_VPt_Har2",100,0.,10.);
234 fHistProVPtHar2->SetXTitle("Pt");
235 fHistProVPtHar2->SetYTitle("v (%)");
236
237 fHistVPtHar2 = new TH1D("Second_FlowLYZ_VPt_Har2","Second_FlowLYZ_VPt_Har2",100,0.,10.);
238 fHistVPtHar2->SetXTitle("Pt");
239 fHistVPtHar2->SetYTitle("v (%)");
240
241 fHistProReDtheta = new TProfile("Second_FlowProLYZ_ReDtheta_Har2","Second_FlowProLYZ_ReDtheta_Har2",fNtheta, -0.5, fNtheta-0.5);
242 fHistProReDtheta->SetXTitle("#theta");
243 fHistProReDtheta->SetYTitle("Re(D^{#theta})");
244
245 fHistProImDtheta = new TProfile("Second_FlowProLYZ_ImDtheta_Har2","Second_FlowProLYZ_ImDtheta_Har2",fNtheta, -0.5, fNtheta-0.5);
246 fHistProImDtheta->SetXTitle("#theta");
247 fHistProImDtheta->SetYTitle("Im(D^{#theta})");
248
249
250 //class AliFlowLYZHist2 defines the histograms:
251 for (Int_t j=0;j<AliFlowConstants::kHars;j++)
252 {
253 for (Int_t theta=0;theta<fNtheta;theta++)
254 {
255 fHist2[j][theta]=new AliFlowLYZHist2(theta,j+1);
256 }
257 }
258
259 //read hists from first run file
260 //firstRunFile = new TFile("fof_flowLYZAnal_firstrun.root","READ"); //default is read
261 if (firstRunFile->IsZombie()){ //check if file exists
262 cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
263 exit(-1);
264 } else if (firstRunFile->IsOpen()){
265 cout<<"----firstRunFile is open----"<<endl<<endl;
266 fHistProVthetaHar1 = (TProfile*)firstRunFile->Get("First_FlowProLYZ_Vtheta_Har1");
267 fHistProVthetaHar2 = (TProfile*)firstRunFile->Get("First_FlowProLYZ_Vtheta_Har2");
268 fHistProR0thetaHar2 = (TProfile*)firstRunFile->Get("First_FlowProLYZ_r0theta_Har2");
269 fHistProV = (TProfile*)firstRunFile->Get("First_FlowProLYZ_V");
270 }
271 }
272
273
274 if (fDebug) cout<<"****Histograms initialised****"<<endl;
275 if (fDebug) cout<<"****fFlowSelect in Init() "<<fFlowSelect<<"****"<<endl;
276
277 fEventNumber = 0; //set event counter to zero
278 /*
279 if (fUseSum)
280 {
281 //initialize LYZ summary class
282 fLYZSummary = new AliFlowLYZSummary();
283 fSummaryFile = new TFile("fFlowSummary.root","RECREATE","Flow LYZ summary file");
284 fSummaryFile->SetCompressionLevel(1);
285 fFlowTree = new TTree("FlowTree", "Flow Summary Tree");
286 fFlowTree->SetAutoSave(1000000); // autosave when 1 Mbyte written
287 fFlowTree->Branch("fLYZSummary","AliFlowLYZSummary",&fLYZSummary,25000,99);
288 }
289 */
290 return kTRUE;
291}
292
293 //-----------------------------------------------------------------------
294
295Bool_t AliFlowLeeYangZerosMaker::Make(AliFlowEvent* fFlowEvent)
296{
297 //make method
298 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Make()****"<<endl;
299
300 //get tracks from event
301 if (fFlowEvent) {
302 fFlowTracks = fFlowEvent->TrackCollection();
303 if (fDebug) cout<<"****fFlowSelect in Make() "<<fFlowSelect<<"****"<<endl;
304 if (fDebug) fFlowSelect->PrintSelectionList() ;
305 if (fDebug) fFlowSelect->PrintList() ;
306
307 if (fFlowSelect && fFlowSelect->Select(fFlowEvent)) // check if event is selected
308 {
309 fFlowTracks = fFlowEvent->TrackCollection(); //get tracks from event
310 fFlowEvent->SetSelections(fFlowSelect) ; // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
311 if (fFirstRun){
312 MakeControlHistograms(fFlowEvent);
313 FillFromFlowEvent(fFlowEvent);
314 }
315 else {
316 MakeControlHistograms(fFlowEvent);
317 SecondFillFromFlowEvent(fFlowEvent);
318 }
319 }
320 }
321 else {
322 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
323 return kFALSE;
324 }
325
326 fEventNumber++;
327
328 return kTRUE;
329}
330
331
332 //-----------------------------------------------------------------------
333 Bool_t AliFlowLeeYangZerosMaker::Finish()
334{
335 //finish method
336 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Finish()****"<<endl;
337
338 //define variables for both runs
339 Float_t fR0 = 0;
340 Float_t fv = 0;
341 Float_t fVtheta = 0;
342 Float_t fSigma2 = 0;
343 Float_t fChi= 0;
344 Float_t fJ01 = 2.405;
345 Int_t fNtheta = AliFlowLYZConstants::kTheta;
346 Float_t fV = 0;
347
348 if (fFirstRun){
349 for (Int_t j=0;j<AliFlowConstants::kHars;j++)
350 //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
351 {
352 Float_t fMeanMult = fHistMult->GetMean();
353 //cerr<<"fMeanMult = "<<fMeanMult<<endl;
354
355 for (Int_t theta=0;theta<fNtheta;theta++)
356 {
357 //get the first minimum r0
358 fHist1[j][theta]->FillGtheta();
359 fR0 = fHist1[j][theta]->GetR0();
360 //cerr<<"fR0 = "<<fR0<<endl;
361
362 //calculate integrated flow
363 if (fR0!=0) fVtheta = fJ01/fR0;
364 if (fMeanMult!=0.)
365 {
366 fv = fVtheta/fMeanMult;
367 cerr<<"fv = "<<fv<<endl;
368 }
369
370
371 //fill the histograms
372 if (j==0)
373 {
374 fHistProR0thetaHar1->Fill(theta,fR0);
375 fHistProVthetaHar1->Fill(theta,fVtheta);
376 fHistProV->Fill(j+1,fVtheta); //profile takes care of the averaging over theta.
377 }
378 if (j==1)
379 {
380 fHistProR0thetaHar2->Fill(theta,fR0);
381 fHistProVthetaHar2->Fill(theta,fVtheta);
382 fHistProV->Fill(j+1,fVtheta); //profile takes care of the averaging over theta.
383 }
384
385 fHistProVR0->Fill(j+1,fv*100); //*100 to get % //profile takes care of the averaging over theta.
386
387 } //end of loop over theta
388
389 //sigma2 and chi
390 if (j==1) //second harmonic only temporarily
391 {
392 if (fEventNumber!=0) {
393 fQsum /= fEventNumber;
394 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
395 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
396 fQ2sum /= fEventNumber;
397 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
398 fV = fHistProV->GetBinContent(j+1);
399 fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.); //BP eq. 62
400 if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
401 else fChi = -1.;
402 cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
403 }
404
405 } //j==1
406
407 } //end of loop over harmonics
408
409 // recalculate statistical errors on integrated flow
410 //combining 5 theta angles to 1 relative error BP eq. 89
411 Double_t fRelErr2comb = 0.;
412 Int_t nEvts = fEventNumber;
413 for (Int_t theta=0;theta<fNtheta;theta++){
414 Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
415 Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
416 TMath::Cos(fTheta));
417 Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
418 TMath::Cos(fTheta));
419 fRelErr2comb += (1/(2*nEvts*(fJ01*fJ01)*TMath::BesselJ1(fJ01)*
420 TMath::BesselJ1(fJ01)))*
421 (fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2)) +
422 fAmincomb*TMath::BesselJ0(2*fJ01)*TMath::Cos(fTheta/2));
423 }
424 fRelErr2comb /= fNtheta;
425 Double_t fRelErrcomb = TMath::Sqrt(fRelErr2comb);
426 cerr<<"fRelErrcomb = "<<fRelErrcomb<<endl;
427
428 //copy content of profile into TH1D and add error
429 for(Int_t b=0;b<2;b++){
430 Double_t fv2pro = fHistProVR0->GetBinContent(b+1);
431 fHistVR0->SetBinContent(b+1,fv2pro);
432 Double_t fv2Err = fv2pro*fRelErrcomb ;
433 cerr<<"fv2pro +- fv2Err = "<<fv2pro<<" +- "<<fv2Err<<" for bin "<<b+1<<endl;
434 fHistVR0->SetBinError(b+1,fv2Err);
435 }
436
437
438 if (fDebug) cout<<"****histograms filled****"<<endl;
439 /*
440 if (fUseSum)
441 {
442 if (fSummaryFile->IsOpen())
443 {
444 fSummaryFile->Write(0,TObject::kOverwrite);
445 fSummaryFile->Close();
446 }
447 }
448 */
449
450 //save histograms in file //temp for testing selector
451 fHistFile->cd();
452 fHistFile->Write();
453
454 return kTRUE;
455 fEventNumber =0; //set to zero for second round over events
456 } //firstrun
457
458
459 else { //second run
460
461 //calculate differential flow
462 //declare variables
463 Float_t fEta, fPt, fReRatio, fVeta, fVPt;
464 Float_t fReDenom = 0;
465 Float_t fImDenom = 0;
466 Double_t fR0 = 0;
467 TComplex fDenom, fNumer, fDtheta;
468 Int_t m = 1;
469 TComplex i = TComplex::I();
470 Float_t fBesselRatio[3] = {1., 1.202, 2.69};
471 Double_t fErrdifcomb = 0.; //set error to zero
472 Double_t fErr2difcomb = 0.; //set error to zero
473
474 /*
475 //v1 integrated
476 for (Int_t j=0;j<AliFlowConstants::kHars;j++) //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
477 {
478 for (Int_t theta=0;theta<fNtheta;theta++)
479 {
480 //get the first minimum r0
481 //fR0 = GetR0(fHist1[j][theta]->FillGtheta());
482
483 fReDenom = fHistProReDenomHar2->GetBinContent(theta+1);
484 fImDenom = fHistProImDenomHar2->GetBinContent(theta+1);
485 TComplex fDenom(fReDenom,fImDenom);
486
487 //complete!!
488
489 }//end of loop over theta
490 }//end of loop over harmonics
491 */
492
493 //differential flow
494 //for (Int_t j=0;j<AliFlowConstants::kHars;j++) //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
495 //{
496 Int_t j=1; //temp only harm 2
497 //fFlowSelect->SetHarmonic(j); //not needed here?
498
499 for (Int_t theta=0;theta<fNtheta;theta++) { //loop over theta
500 if (j==0) {
501 fR0 = fHistProR0thetaHar1->GetBinContent(theta+1);
502 fVtheta = 2.405/fR0;
503 //fVtheta = fHistProVthetaHar1->GetBinContent(theta+1); // BP Eq. 9 -> Vn^theta = j01/r0^theta
504 fReDenom = fHistProReDenomHar1->GetBinContent(theta+1);
505 fImDenom = fHistProImDenomHar1->GetBinContent(theta+1);
506 }
507 if (j==1) {
508 fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
509 if (fDebug) cerr<<"fR0 = "<<fR0<<endl;
510 fVtheta = 2.405/fR0; // BP Eq. 9 -> Vn^theta = j01/r0^theta
511 fReDenom = fHistProReDenomHar2->GetBinContent(theta+1);
512 fImDenom = fHistProImDenomHar2->GetBinContent(theta+1);
513 }
514
515 fDenom(fReDenom,fImDenom);
516
517
518 //for new method and use by others (only with the sum generating function):
519 if (fUseSum) {
520 fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
521 fDtheta = fR0*fDenom;
522 fHistProReDtheta->Fill(theta,fDtheta.Re());
523 fHistProImDtheta->Fill(theta,fDtheta.Im());
524 }
525
526 fDenom *= TComplex::Power(i, m-1);
527 //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
528
529 //v as a function of eta
530 Int_t fEtaBins = AliFlowLYZConstants::kEtaBins;
531 for (Int_t be=1;be<=fEtaBins;be++) {
532 fEta = fHist2[j][theta]->GetBinCenter(be);
533 fNumer = fHist2[j][theta]->GetfNumer(be);
534 if (fNumer.Rho()==0) {
535 if (fDebug) cerr<<"WARNING: modulus of fNumer is zero in Finish()"<<endl;
536 fReRatio = 0;
537 }
538 else if (fDenom.Rho()==0) {
539 if (fDebug) cerr<<"WARNING: modulus of fDenom is zero"<<endl;
540 fReRatio = 0;
541 }
542 else {
543 //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
544 fReRatio = (fNumer/fDenom).Re();
545 }
546
547 fVeta = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
548 //if ( j==1 && theta==0) cerr<<"eta = "<<fEta<<" cerr::fReRatio for eta = "<<fReRatio<<" cerr::fVeta for eta = "<<fVeta<<endl;
549
550 if (j==0) fHistProVetaHar1->Fill(fEta,fVeta);
551 if (j==1) fHistProVetaHar2->Fill(fEta,fVeta);
552 } //loop over bins be
553
554 //v as a function of Pt
555 Int_t fPtBins = AliFlowLYZConstants::kPtBins;
556 for (Int_t bp=1;bp<=fPtBins;bp++) {
557 fPt = fHist2[j][theta]->GetBinCenterPt(bp);
558 fNumer = fHist2[j][theta]->GetfNumerPt(bp);
559 if (fNumer.Rho()==0) {
560 if (fDebug) cerr<<"modulus of fNumer is zero"<<endl;
561 fReRatio = 0;
562 }
563 else if (fDenom.Rho()==0) {
564 if (fDebug) cerr<<"modulus of fDenom is zero"<<endl;
565 fReRatio = 0;
566 }
567 else {
568 //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
569 fReRatio = (fNumer/fDenom).Re();
570 }
571
572 fVPt = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
573 //cerr<<"fBesselRatio[m-1] = "<<fBesselRatio[m-1]<<endl; //checked ok
574 //if ( j==1 && theta==0) cerr<<"pt = "<<fPt<<" cerr::fReRatio for pt = "<<fReRatio<<" cerr::fVPt for pt = "<<fVeta<<endl;
575
576 if (j==0) fHistProVPtHar1->Fill(fPt,fVPt);
577 if (j==1) fHistProVPtHar2->Fill(fPt,fVPt);
578 } //loop over bins bp
579
580 }//end of loop over theta
581
582
583 //sigma2 and chi (for statistical error calculations)
584 if (j==1) { //second harmonic only temporarily
585
586 if (fEventNumber!=0) {
587 fQsum /= fEventNumber;
588 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
589 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
590 fQ2sum /= fEventNumber;
591 cerr<<"fQ2sum = "<<fQ2sum<<endl;
592 fV = fHistProV->GetBinContent(j+1);
593 fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.); //BP eq. 62
594 if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
595 else fChi = -1.;
596 cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
597 }
598
599 } //j==1
600
601
602 //copy content of profile into TH1D and add error
603 for(Int_t b=0;b<100;b++){
604 Double_t fv2pro = fHistProVPtHar2->GetBinContent(b);
605 //calculate error
606 for (Int_t theta=0;theta<fNtheta;theta++) {
607 Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
608 Int_t Nprime = fHist2[j][theta]->GetNprimePt(b);
609 //cerr<<"Nprime = "<<Nprime<<endl;
610 if (Nprime!=0.) {
611 Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
612 TMath::Cos(fTheta));
613 Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
614 TMath::Cos(fTheta));
615 fErr2difcomb += (TMath::Cos(fTheta)/(4*Nprime*TMath::BesselJ1(fJ01)*
616 TMath::BesselJ1(fJ01)))*
617 ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) -
618 (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
619 }
620 } //loop over theta
621
622 if (fErr2difcomb!=0.) {
623 fErr2difcomb /= fNtheta;
624 fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
625 //cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
626 }
627 else {fErrdifcomb = 0.;}
628
629 //fill TH1D
630 if (j==1) {
631 fHistVPtHar2->SetBinContent(b,fv2pro);
632 fHistVPtHar2->SetBinError(b,fErrdifcomb);
633 }
634 //check that error is set
635 //if (j==1) { cout<<"difference between calculated error and error in hostogram: "<< fErrdifcomb - fHistVPtHar2->GetBinError(b)<<endl; }
636
637 } //loop over bins b
638
639
640 //} //end of loop over harmonics //temporarily out
641
642 //save histograms in file
643 fHistFile->cd();
644 fHistFile->Write();
645 fHistFile->Close();
646 //Note that when the file is closed, all histograms and Trees in memory associated with this file are deleted
647 if (fDebug) cout<<"****Histograms saved and fHistFile closed, all histograms deleted****"<<endl;
648
649 //close the first run file
650 firstRunFile->Close();
651
652
653 } //secondrun
654
655 delete fFlowSelect;
656 cout<<"----LYZ analysis finished....----"<<endl<<endl;
657
658 return kTRUE;
659}
660
661
662//-----------------------------------------------------------------------
663 Bool_t AliFlowLeeYangZerosMaker::MakeControlHistograms(AliFlowEvent* fFlowEvent)
664{
665 //contol histograms of pt, eta, phi, Q, mult
666 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::MakeControlHistograms()****"<<endl;
667
668 //define variables
669 TVector2 fQ;
670 Float_t fPt, fPhi, fEta, fQmult;
671
672 if (!fFlowEvent){
673 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
674 return kFALSE;
675 }
676
677 if (!fFlowSelect){
678 cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
679 return kFALSE;
680 }
681
682 //set selection and harmonic
683 fFlowSelect->SetSelection(1);
684 fFlowSelect->SetHarmonic(1); //second harmonic
685
686 //cerr<<"selection in MakeControlHistograms()"<<endl;
687 //fFlowSelect->PrintSelectionList() ;
688 //fFlowSelect->PrintList() ;
689
690 fFlowTracks = fFlowEvent->TrackCollection();
691 Int_t fNumberOfTracks = fFlowTracks->GetEntries();
692 fHistOrigMult->Fill(fNumberOfTracks);
693 //cerr<<"fNumberOfTracks = "<<fNumberOfTracks<<endl;
694 Int_t fMult = fFlowEvent->Mult(fFlowSelect); // Multiplicity of tracks in the specified Selection
695 fHistMult->Fill(fMult);
696 //cerr<<"Mult = "<<fMult<<endl;
697
698 fQ = fFlowEvent ->Q(fFlowSelect);
699 fQmult = fQ.Mod()/TMath::Sqrt(fNumberOfTracks);
700 fHistQ->Fill(fQmult);
701
702 Int_t tempmult = 0; //for testing
703 for (Int_t i=0;i<fNumberOfTracks;i++)
704 {
705 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ; //get object at array position i
706 //if (fFlowSelect->SelectPart(fFlowTrack)) //if track is selected
707 if (fFlowSelect->Select(fFlowTrack))
708 {
709 fPt = fFlowTrack->Pt();
710 fHistPt->Fill(fPt);
711 tempmult++;
712 fPhi = fFlowTrack->Phi();
713 if (fPhi<0.) fPhi+=2*TMath::Pi();
714 fHistPhi->Fill(fPhi);
715 fEta = fFlowTrack->Eta();
716 fHistEta->Fill(fEta);
717 }
718 }
719 if (fMult!=tempmult){cerr<<"ERROR: Mult() is not tempmult! "<<fMult<<" :: "<<tempmult<<endl<<endl;}
720 //else {cerr<<"Mult()= tempmult "<<fMult<<" :: "<<tempmult<<endl<<endl;}
721
722 return kTRUE;
723
724
725}
726
727
728
729//-----------------------------------------------------------------------
730
731 Bool_t AliFlowLeeYangZerosMaker::FillFromFlowEvent(AliFlowEvent* fFlowEvent)
732{
733 // Get event quantities from AliFlowEvent for all particles
734
735 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::FillFromFlowEvent()****"<<endl;
736
737 if (!fFlowEvent){
738 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
739 return kFALSE;
740 }
741
742 if (!fFlowSelect){
743 cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
744 return kFALSE;
745 }
746
747
748 //define variables
749 TComplex fExpo, fGtheta, fGthetaNew, fZ;
750 //Int_t m;
751 Int_t fNtheta = AliFlowLYZConstants::kTheta;
752 Int_t fNbins = AliFlowLYZConstants::kNbins;
753
754
755 //calculate flow
756 fFlowSelect->SetSelection(1);
757
758 for (Int_t j=0;j<AliFlowConstants::kHars;j++) //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
759 {
760 Int_t m=1;
761 fFlowSelect->SetHarmonic(j);
762 Float_t fOrder = (double)((j+1)/m);
763 //cerr<<"fOrder = "<<fOrder<<endl;
764
765 //get the Q vector from the FlowEvent
766 TVector2 fQ = fFlowEvent->Q(fFlowSelect);
767 //for chi calculation:
768 if (j==1) //second harmonic only temporarily
769 {
770 fQsum += fQ;
771 fQ2sum += fQ.Mod2();
772 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
773 }
774
775 fFlowTracks = fFlowEvent->TrackCollection();
776
777 for (Int_t theta=0;theta<fNtheta;theta++)
778 {
779 fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;
780 //cerr<<"fTheta = "<<fTheta<<endl;
781
782 //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
783 fQtheta = GetQtheta(fFlowSelect,fFlowTracks,fTheta);
784 //save fQtheta in AliFlowLYZSummary class
785
786 //something
787 //AliFlowLYZSummary::SetQtheta(theta,fQtheta);
788
789 if (j==1 && theta==0) fHistQtheta->Fill(fQtheta);
790 //cerr<<"fQtheta = "<<fQtheta<<endl;
791
792 for (Int_t bin=1;bin<=fNbins;bin++)
793 {
794 Float_t fR = fHist1[j][theta]->GetBinCenter(bin); //bincentre of bins in histogram
795 //if (theta == 0) cerr<<"cerr::fR = "<<fR<<endl;
796 if (fUseSum)
797 {
798 //calculate the sum generating function
799 fExpo(0.,fR*fQtheta); //Re=0 ; Im=fR*fQtheta
800 fGtheta = TComplex::Exp(fExpo);
801 }
802 else
803 {
804 //calculate the product generating function
805 fGtheta = GetGrtheta(fFlowSelect,fR,fTheta); //make this function
806 if (fGtheta.Rho2() > 100.) break;
807 }
808
809 fHist1[j][theta]->Fill(fR,fGtheta); //fill real and imaginary part of fGtheta
810
811 } //loop over bins
812 } //loop over theta
813 } //loop over harmonics
814
815 return kTRUE;
816
817
818}
819
820 //-----------------------------------------------------------------------
821 Bool_t AliFlowLeeYangZerosMaker::SecondFillFromFlowEvent(AliFlowEvent* fFlowEvent)
822{
823 //for differential flow
824
825 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::SecondFillFromFlowEvent()****"<<endl;
826
827 if (!fFlowEvent){
828 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
829 return kFALSE;
830 }
831
832 if (!fFlowSelect){
833 cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
834 return kFALSE;
835 }
836
837 //define variables
838 //TVector2 fQ;
839 TComplex fExpo, fDenom, fNumer,fCosTermComplex;
840 Float_t fOrder, fR0, fPhi, fEta, fPt;
841 Double_t fCosTerm;
842 Double_t m = 1.;
843 Int_t fNtheta = AliFlowLYZConstants::kTheta;
844
845 //calculate flow
846 fFlowSelect->SetSelection(1);
847
848 //cerr<<"selection in SecondFillFromFlowEvent()"<<endl;
849 //fFlowSelect->PrintSelectionList() ;
850 //fFlowSelect->PrintList() ;
851
852
853 for (Int_t j=0;j<AliFlowConstants::kHars;j++) //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
854 {
855 m=1;
856 fFlowSelect->SetHarmonic(j);
857 fOrder = (double)((j+1)/m);
858
859 //get the Q vector from the FlowEvent
860 TVector2 fQ = fFlowEvent->Q(fFlowSelect);
861 //for chi calculation:
862 if (j==1) //second harmonic only temporarily
863 {
864 fQsum += fQ;
865 fQ2sum += fQ.Mod2();
866 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
867 }
868
869 fFlowTracks = fFlowEvent->TrackCollection();
870
871 for (Int_t theta=0;theta<fNtheta;theta++)
872 {
873 fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;
874
875 //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
876 fQtheta = GetQtheta(fFlowSelect,fFlowTracks,fTheta);
877 //cerr<<"fQtheta for fdenom = "<<fQtheta<<endl;
878
879 /*
880 if (j==0) //first harmonic
881 {
882 //denominator for differential v
883 fR0 = fHistProR0thetaHar1->GetBinContent(theta+1);
884 fExpo(0.,fR0*fQtheta);
885 fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
886 //cerr<<"fDenom.Re() = "<<fDenom.Re()<<endl;
887 //cerr<<"fDenom.Im() = "<<fDenom.Im()<<endl;
888
889 //denominator for differential v
890 // ****put in product generating function!!
891
892
893 fHistProReDenomHar1->Fill(theta,fDenom.Re()); //fill the real part of fDenom
894 fHistProImDenomHar1->Fill(theta,fDenom.Im()); //fill the imaginary part of fDenom
895 }
896 */
897
898 if (j==1) //second harmonic
899 {
900 //denominator for differential v
901 fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
902 //cerr<<"fR0 = "<<fR0 <<endl;
903
904 if (fUseSum) //sum generating function
905 {
906 fExpo(0.,fR0*fQtheta);
907 fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
908 //loop over tracks in event
909 fFlowTracks = fFlowEvent->TrackCollection();
910 Int_t fNumberOfTracks = fFlowTracks->GetEntries();
911 for (Int_t i=0;i<fNumberOfTracks;i++)
912 {
913 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ; //get object at array position i
914 if (fFlowSelect->SelectPart(fFlowTrack)) //if track is selected
915 {
916 fPhi = fFlowTrack->Phi();
917 fCosTerm = cos(m*fOrder*(fPhi-fTheta));
918 //cerr<<"fCosTerm = "<<fCosTerm <<endl;
919 fNumer = fCosTerm*(TComplex::Exp(fExpo));
920 if (fNumer.Rho()==0) {cerr<<"WARNING: modulus of fNumer is zero in SecondFillFromFlowEvent"<<endl;}
921 if (fDebug) cerr<<"modulus of fNumer is "<<fNumer.Rho()<<endl;
922 fEta = fFlowTrack->Eta(); //rapidity
923 fPt = fFlowTrack->Pt();
924 fHist2[j][theta]->Fill(fEta,fPt,fNumer);
925 } //if
926 } //loop over tracks
927
928 } //sum
929 else //product generating function
930 {
931 fDenom = GetDiffFlow(fFlowSelect,fR0,theta);
932
933 }//product
934
935 fHistProReDenomHar2->Fill(theta,fDenom.Re()); //fill the real part of fDenom
936 fHistProImDenomHar2->Fill(theta,fDenom.Im()); //fill the imaginary part of fDenom
937 }//j==1
938
939
940 }//end of loop over theta
941
942 }//loop over harmonics
943
944
945 return kTRUE;
946
947
948
949}
950 //-----------------------------------------------------------------------
951 Double_t AliFlowLeeYangZerosMaker::GetQtheta(AliFlowSelection* fFlowSelect, TObjArray* fFlowTracks, Float_t fTheta)
952{
953 //calculate Qtheta. Qtheta is the sum over all particles of cos(fOrder*(fPhi-fTheta)) BP eq. 3
954 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetQtheta()****"<<endl;
955
956 if (!fFlowSelect){
957 cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
958 return kFALSE;
959 }
960
961
962 Double_t fQtheta = 0.;
963 Int_t fHarmonic = fFlowSelect->Har();
964 Double_t fOrder = (double)(fHarmonic+1);
965
966 Int_t fNumberOfTracks = fFlowTracks->GetEntries();
967 //cerr<<"GetQtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
968
969 for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
970 {
971 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ; //get object at array position i
972 if(fFlowSelect->Select(fFlowTrack)) //if track is selected //gives the same number of particles as Mult(fFlowSelect) method
973 {
974 Float_t fPhi = fFlowTrack->Phi();
975 fQtheta += cos(fOrder*(fPhi-fTheta));
976 }
977
978 }//loop over tracks
979
980 return fQtheta;
981
982}
983
984//-----------------------------------------------------------------------
985TComplex AliFlowLeeYangZerosMaker::GetGrtheta(AliFlowSelection* fFlowSelect, Float_t fR, Float_t fTheta)
986{
987 // Product Generating Function for LeeYangZeros method
988 // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
989
990 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetGrtheta()****"<<endl;
991
992 if (!fFlowSelect){
993 cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
994 return kFALSE;
995 }
996
997 TComplex fG = TComplex::One();
998 Int_t fHarmonic = fFlowSelect->Har();
999 Double_t fOrder = (double)(fHarmonic+1);
1000 Double_t fWgt = 1.;
1001
1002 Int_t fNumberOfTracks = fFlowTracks->GetEntries();
1003 //cerr<<"GetGrtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
1004
1005 for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
1006 {
1007 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ; //get object at array position i
1008 if (fFlowSelect->SelectPart(fFlowTrack)) //if track is selected
1009 {
1010 Float_t fPhi = fFlowTrack->Phi();
1011 Double_t fGIm = fR * fWgt*cos(fOrder*(fPhi - fTheta));
1012 TComplex fGi(1., fGIm);
1013 fG *= fGi; //product over all tracks
1014 }//if
1015 }//loop over tracks
1016
1017 return fG;
1018
1019 }
1020
1021
1022//-----------------------------------------------------------------------
1023TComplex AliFlowLeeYangZerosMaker::GetDiffFlow(AliFlowSelection* fFlowSelect, Float_t fR0, Int_t theta)
1024{
1025 // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1026 // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1027 // Also for v1 mixed harmonics: DF Eq. 5
1028 // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1029
1030 if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetGrtheta()****"<<endl;
1031
1032 if (!fFlowSelect){
1033 cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
1034 return kFALSE;
1035 }
1036
1037
1038 TComplex fG = TComplex::One();
1039 TComplex fdGr0(0.,0.);
1040 Int_t fHarmonic = fFlowSelect->Har();
1041 Double_t fOrder = (double)(fHarmonic+1);
1042 Double_t fWgt = 1.;
1043
1044 Int_t fNumberOfTracks = fFlowTracks->GetEntries();
1045 //cerr<<"GetGrtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
1046
1047 Int_t fNtheta = AliFlowLYZConstants::kTheta;
1048 Float_t fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;
1049
1050 for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
1051 {
1052 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ; //get object at array position i
1053 if (fFlowSelect->SelectPart(fFlowTrack)) //if track is selected
1054 {
1055 Float_t fPhi = fFlowTrack->Phi();
1056 Double_t fCosTerm = fWgt*cos(fOrder*(fPhi - fTheta));
1057 //GetGr0theta
1058 Double_t fGIm = fR0 * fCosTerm;
1059 TComplex fGi(1., fGIm);
1060 fG *= fGi; //product over all tracks
1061 //GetdGr0theta
1062 TComplex fCosTermComplex(1., fR0*fCosTerm);
1063 fdGr0 +=(fCosTerm / fCosTermComplex); //sum over all tracks
1064 }//if
1065 }//loop over tracks
1066
1067 for (Int_t i=0;i<fNumberOfTracks;i++)
1068 {
1069 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ; //get object at array position i
1070 if (fFlowSelect->SelectPart(fFlowTrack)) //if track is selected
1071 {
1072 Float_t fPhi = fFlowTrack->Phi();
1073 Double_t fCosTerm = cos(fOrder*(fPhi-fTheta));
1074 TComplex fCosTermComplex(1.,fR0*fCosTerm);
1075
1076 TComplex fNumer = fG*fCosTerm/fCosTermComplex; //PG Eq. 9
1077 Float_t fEta = fFlowTrack->Eta(); //rapidity
1078 Float_t fPt = fFlowTrack->Pt();
1079 fHist2[1][theta]->Fill(fEta,fPt,fNumer);
1080 }//if
1081 }//loop over tracks
1082
1083 TComplex fDenom = fG*fdGr0;
1084 return fDenom;
1085
1086 }
1087
1088//-------------------------------------------------