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