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