AliPhysics  b555aef (b555aef)
AliFlowAnalysisWithLeeYangZeros.cxx
Go to the documentation of this file.
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 
17 //Description: Maker to analyze Flow using the LeeYangZeros method
18 // Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
19 // Practical Guide (PG): J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
20 // Adapted from StFlowLeeYangZerosMaker.cxx
21 // by Markus Oldenberg and Art Poskanzer, LBNL
22 // with advice from Jean-Yves Ollitrault and Nicolas Borghini
23 //
24 //Author: Naomi van der Kolk (kolk@nikhef.nl)
26 
27 #include "Riostream.h"
28 #include "TObject.h"
29 #include "TMath.h"
30 #include "TFile.h"
31 #include "TList.h"
32 #include "TVector2.h"
33 #include "TH1F.h"
34 #include "TComplex.h"
35 #include "TProfile.h"
36 #include "TDirectoryFile.h"
37 
38 #include "AliFlowCommonConstants.h"
39 #include "AliFlowLYZConstants.h"
41 #include "AliFlowLYZHist1.h"
42 #include "AliFlowLYZHist2.h"
43 #include "AliFlowCommonHist.h"
45 #include "AliFlowEventSimple.h"
46 #include "AliFlowTrackSimple.h"
47 #include "AliFlowVector.h"
48 
49 using std::endl;
50 using std::cout;
51 using std::cerr;
53 
54  //-----------------------------------------------------------------------
55 
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  fHistVtheta(NULL),
67  fHistProVetaRP(NULL),
68  fHistProVetaPOI(NULL),
69  fHistProVPtRP(NULL),
70  fHistProVPtPOI(NULL),
71  fHistR0theta(NULL),
72  fHistProReDenom(NULL),
73  fHistProImDenom(NULL),
74  fHistReDtheta(NULL),
75  fHistImDtheta(NULL),
76  fHistQsumforChi(NULL),
77  fCommonHists(NULL),
78  fCommonHistsRes(NULL)
79 
80 {
81  //default constructor
82  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
83 
84  fHistList = new TList();
85  fFirstRunList = new TList();
86 
87  for(Int_t i = 0;i<5;i++)
88  {
89  fHist1[i]=0;
90  fHist2RP[i]=0;
91  fHist2POI[i]=0;
92  }
93 
94  fQsum = new TVector2();
95 
96 }
97 
98 //-----------------------------------------------------------------------
99 
100 
102  {
103  //default destructor
104  if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
105  delete fQsum;
106  delete fHistList;
107  delete fFirstRunList;
108 
109  }
110 
111 //-----------------------------------------------------------------------
112 
114 {
115  //store the final results in output .root file
116 
117  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
118  if (GetFirstRun()) {
119  //output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
120  if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
121  else {fHistList->SetName("cobjLYZ1PROD");}
122  fHistList->SetOwner(kTRUE);
123  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
124  }
125  else {
126  //output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
127  if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
128  else { fHistList->SetName("cobjLYZ2PROD"); }
129  fHistList->SetOwner(kTRUE);
130  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
131  }
132  delete output;
133 }
134 
135 //-----------------------------------------------------------------------
136 
138 {
139  //store the final results in output .root file
140 
141  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
142  if (GetFirstRun()) {
143  //output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
144  if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
145  else {fHistList->SetName("cobjLYZ1PROD");}
146  fHistList->SetOwner(kTRUE);
147  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
148  }
149  else {
150  //output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
151  if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
152  else { fHistList->SetName("cobjLYZ2PROD"); }
153  fHistList->SetOwner(kTRUE);
154  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
155  }
156  delete output;
157 }
158 
159 //-----------------------------------------------------------------------
160 
161 void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TDirectoryFile *outputFileName)
162 {
163  //store the final results in output .root file
164  if (GetFirstRun()) {
165  if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
166  else {fHistList->SetName("cobjLYZ1PROD");}
167  fHistList->SetOwner(kTRUE);
168  outputFileName->Add(fHistList);
169  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
170  }
171  else {
172  if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
173  else { fHistList->SetName("cobjLYZ2PROD"); }
174  fHistList->SetOwner(kTRUE);
175  outputFileName->Add(fHistList);
176  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
177  }
178 }
179 
180 //-----------------------------------------------------------------------
181 
183 {
184  //init method
185  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
186 
187  //save old value and prevent histograms from being added to directory
188  //to avoid name clashes in case multiple analaysis objects are used
189  //in an analysis
190  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
191  TH1::AddDirectory(kFALSE);
192 
193  // Book histograms
197 
202 
203  //for control histograms
204  if (fFirstRun){
205  if (fUseSum) {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1SUM");}
206  else {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1PROD");}
207  fHistList->Add(fCommonHists);
208  if (fUseSum) {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1SUM");}
209  else {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1PROD");}
211  }
212  else {
213  if (fUseSum) {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2SUM");}
214  else {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2PROD");}
215  fHistList->Add(fCommonHists);
216  if (fUseSum) {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2SUM");}
217  else {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2PROD");}
218  fHistList->Add(fCommonHistsRes);
219  }
220 
221  TString nameChiHist;
222  if (fUseSum) {nameChiHist = "Flow_QsumforChi_LYZSUM";}
223  else {nameChiHist = "Flow_QsumforChi_LYZPROD";}
224  fHistQsumforChi = new TH1F(nameChiHist.Data(),nameChiHist.Data(),3,-1.,2.);
225  fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
226  fHistQsumforChi->SetYTitle("value");
228 
229  //for first loop over events
230  if (fFirstRun){
231  TString nameR0Hist;
232  if (fUseSum) {nameR0Hist = "First_Flow_r0theta_LYZSUM";}
233  else {nameR0Hist = "First_Flow_r0theta_LYZPROD";}
234  fHistR0theta = new TH1D(nameR0Hist.Data(),nameR0Hist.Data(),iNtheta,-0.5,iNtheta-0.5);
235  fHistR0theta->SetXTitle("#theta");
236  fHistR0theta->SetYTitle("r_{0}^{#theta}");
237  fHistList->Add(fHistR0theta);
238 
239  TString nameVHist;
240  if (fUseSum) {nameVHist = "First_Flow_Vtheta_LYZSUM";}
241  else {nameVHist = "First_Flow_Vtheta_LYZPROD";}
242  fHistVtheta = new TH1D(nameVHist.Data(),nameVHist.Data(),iNtheta,-0.5,iNtheta-0.5);
243  fHistVtheta->SetXTitle("#theta");
244  fHistVtheta->SetYTitle("V_{n}^{#theta}");
245  fHistList->Add(fHistVtheta);
246 
247  //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta
248  for (Int_t theta=0;theta<iNtheta;theta++) {
249  TString nameHist1;
250  if (fUseSum) {nameHist1 = "AliFlowLYZHist1_";}
251  else {nameHist1 = "AliFlowLYZHist1_";}
252  nameHist1 += theta;
253  fHist1[theta]=new AliFlowLYZHist1(theta, nameHist1, fUseSum);
254  fHistList->Add(fHist1[theta]);
255  }
256 
257  }
258  //for second loop over events
259  else {
260  TString nameReDenomHist;
261  if (fUseSum) {nameReDenomHist = "Second_FlowPro_ReDenom_LYZSUM";}
262  else {nameReDenomHist = "Second_FlowPro_ReDenom_LYZPROD";}
263  fHistProReDenom = new TProfile(nameReDenomHist.Data(),nameReDenomHist.Data(), iNtheta, -0.5, iNtheta-0.5);
264  fHistProReDenom->SetXTitle("#theta");
265  fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
267 
268  TString nameImDenomHist;
269  if (fUseSum) {nameImDenomHist = "Second_FlowPro_ImDenom_LYZSUM";}
270  else {nameImDenomHist = "Second_FlowPro_ImDenom_LYZPROD";}
271  fHistProImDenom = new TProfile(nameImDenomHist.Data(),nameImDenomHist.Data(), iNtheta, -0.5, iNtheta-0.5);
272  fHistProImDenom->SetXTitle("#theta");
273  fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
275 
276  TString nameVetaRPHist;
277  if (fUseSum) {nameVetaRPHist = "Second_FlowPro_VetaRP_LYZSUM";}
278  else {nameVetaRPHist = "Second_FlowPro_VetaRP_LYZPROD";}
279  fHistProVetaRP = new TProfile(nameVetaRPHist.Data(),nameVetaRPHist.Data(),iNbinsEta,dEtaMin,dEtaMax);
280  fHistProVetaRP->SetXTitle("rapidity");
281  fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
283 
284  TString nameVetaPOIHist;
285  if (fUseSum) {nameVetaPOIHist = "Second_FlowPro_VetaPOI_LYZSUM";}
286  else {nameVetaPOIHist = "Second_FlowPro_VetaPOI_LYZPROD";}
287  fHistProVetaPOI = new TProfile(nameVetaPOIHist.Data(),nameVetaPOIHist.Data(),iNbinsEta,dEtaMin,dEtaMax);
288  fHistProVetaPOI->SetXTitle("rapidity");
289  fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
291 
292  TString nameVPtRPHist;
293  if (fUseSum) {nameVPtRPHist = "Second_FlowPro_VPtRP_LYZSUM";}
294  else {nameVPtRPHist = "Second_FlowPro_VPtRP_LYZPROD";}
295  fHistProVPtRP = new TProfile(nameVPtRPHist.Data(),nameVPtRPHist.Data(),iNbinsPt,dPtMin,dPtMax);
296  fHistProVPtRP->SetXTitle("Pt");
297  fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
298  fHistList->Add(fHistProVPtRP);
299 
300  TString nameVPtPOIHist;
301  if (fUseSum) {nameVPtPOIHist = "Second_FlowPro_VPtPOI_LYZSUM";}
302  else {nameVPtPOIHist = "Second_FlowPro_VPtPOI_LYZPROD";}
303  fHistProVPtPOI = new TProfile(nameVPtPOIHist.Data(),nameVPtPOIHist.Data(),iNbinsPt,dPtMin,dPtMax);
304  fHistProVPtPOI->SetXTitle("p_{T}");
305  fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
307 
308  if (fUseSum){
309  fHistReDtheta = new TH1D("Second_Flow_ReDtheta_LYZSUM","Second_Flow_ReDtheta_LYZSUM",iNtheta, -0.5, (double)iNtheta-0.5);
310  fHistReDtheta->SetXTitle("#theta");
311  fHistReDtheta->SetYTitle("Re(D^{#theta})");
312  fHistList->Add(fHistReDtheta);
313 
314  fHistImDtheta = new TH1D("Second_Flow_ImDtheta_LYZSUM","Second_Flow_ImDtheta_LYZSUM",iNtheta, -0.5, (double)iNtheta-0.5);
315  fHistImDtheta->SetXTitle("#theta");
316  fHistImDtheta->SetYTitle("Im(D^{#theta})");
317  fHistList->Add(fHistImDtheta);
318  }
319 
320  //class AliFlowLYZHist2 defines the histograms:
321  for (Int_t theta=0;theta<iNtheta;theta++) {
322  TString nameRP = "AliFlowLYZHist2RP_";
323  nameRP += theta;
324  fHist2RP[theta]=new AliFlowLYZHist2(theta, "RP", nameRP, fUseSum);
325  fHistList->Add(fHist2RP[theta]);
326 
327  TString namePOI = "AliFlowLYZHist2POI_";
328  namePOI += theta;
329  fHist2POI[theta]=new AliFlowLYZHist2(theta, "POI", namePOI, fUseSum);
330  fHistList->Add(fHist2POI[theta]);
331  }
332 
333  //read histogram fHistR0theta from the first run list
334  if (fFirstRunList) {
335  if (fUseSum) { fHistR0theta = (TH1D*)fFirstRunList->FindObject("First_Flow_r0theta_LYZSUM");}
336  else{ fHistR0theta = (TH1D*)fFirstRunList->FindObject("First_Flow_r0theta_LYZPROD");}
337  if (!fHistR0theta) {cout<<"fHistR0theta has a NULL pointer!"<<endl;}
338  fHistList->Add(fHistR0theta);
339  } else { cout<<"list is NULL pointer!"<<endl; }
340 
341  }
342 
343 
344  if (fDebug) cout<<"****Histograms initialised****"<<endl;
345 
346  fEventNumber = 0; //set event counter to zero
347 
348  //resore old stuff
349  TH1::AddDirectory(oldHistAddStatus);
350 
351  return kTRUE;
352 }
353 
354  //-----------------------------------------------------------------------
355 
357 {
358  //make method
359  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
360 
361  //get tracks from event
362  if (anEvent) {
363  if (fFirstRun){
365  FillFromFlowEvent(anEvent);
366  }
367  else {
369  SecondFillFromFlowEvent(anEvent);
370  }
371  }
372 
373  else {
374  cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
375  return kFALSE;
376  }
377  // cout<<"^^^^read event "<<fEventNumber<<endl;
378  fEventNumber++;
379 
380 
381  return kTRUE;
382 }
383 
384  //-----------------------------------------------------------------------
386  // get the pointers to all output histograms before calling Finish()
387 
388  const Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
389 
390  if (outputListHistos) {
391 
392  //define histograms for first and second run
393  AliFlowCommonHist *pCommonHist = NULL;
394  AliFlowCommonHistResults *pCommonHistResults = NULL;
395  TH1D* pHistR0theta = NULL;
396  TH1D* pHistVtheta = NULL;
397  TProfile* pHistProReDenom = NULL;
398  TProfile* pHistProImDenom = NULL;
399  TProfile* pHistProVetaRP = NULL;
400  TProfile* pHistProVetaPOI = NULL;
401  TProfile* pHistProVPtRP = NULL;
402  TProfile* pHistProVPtPOI = NULL;
403  TH1F* pHistQsumforChi = NULL;
404  //AliFlowLYZHist1 *pLYZHist1[iNtheta] = {NULL}; //array of pointers to AliFlowLYZHist1
405  //AliFlowLYZHist2 *pLYZHist2RP[iNtheta] = {NULL}; //array of pointers to AliFlowLYZHist2
406  //AliFlowLYZHist2 *pLYZHist2POI[iNtheta] = {NULL}; //array of pointers to AliFlowLYZHist2
407  AliFlowLYZHist1 **pLYZHist1 = new AliFlowLYZHist1*[iNtheta]; //array of pointers to AliFlowLYZHist1
408  AliFlowLYZHist2 **pLYZHist2RP = new AliFlowLYZHist2*[iNtheta]; //array of pointers to AliFlowLYZHist2
409  AliFlowLYZHist2 **pLYZHist2POI = new AliFlowLYZHist2*[iNtheta]; //array of pointers to AliFlowLYZHist2
410  for (Int_t i=0; i<iNtheta; i++)
411  {
412  pLYZHist1[i] = NULL;
413  pLYZHist2RP[i] = NULL;
414  pLYZHist2POI[i] = NULL;
415  }
416 
417  if (GetFirstRun()) { //first run
418  //Get the common histograms from the output list
419  if (GetUseSum()){
420  pCommonHist = dynamic_cast<AliFlowCommonHist*>
421  (outputListHistos->FindObject("AliFlowCommonHistLYZ1SUM"));
422  pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
423  (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ1SUM"));
424  //Get the histograms from the output list
425  for(Int_t theta = 0;theta<iNtheta;theta++){
426  TString name = "AliFlowLYZHist1_";
427  name += theta;
428  pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
429  (outputListHistos->FindObject(name));
430  }
431  pHistVtheta = dynamic_cast<TH1D*>
432  (outputListHistos->FindObject("First_Flow_Vtheta_LYZSUM"));
433 
434  pHistR0theta = dynamic_cast<TH1D*>
435  (outputListHistos->FindObject("First_Flow_r0theta_LYZSUM"));
436 
437  pHistQsumforChi = dynamic_cast<TH1F*>
438  (outputListHistos->FindObject("Flow_QsumforChi_LYZSUM"));
439 
440  //Set the histogram pointers and call Finish()
441  if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
442  pHistVtheta && pHistR0theta && pHistQsumforChi ) {
443  this->SetCommonHists(pCommonHist);
444  this->SetCommonHistsRes(pCommonHistResults);
445  this->SetHist1(pLYZHist1);
446  this->SetHistVtheta(pHistVtheta);
447  this->SetHistR0theta(pHistR0theta);
448  this->SetHistQsumforChi(pHistQsumforChi);
449  }
450  else {
451  cout<<"WARNING: Histograms needed to run Finish() firstrun (SUM) are not accessable!"<<endl;
452  }
453  }
454  else {
455  pCommonHist = dynamic_cast<AliFlowCommonHist*>
456  (outputListHistos->FindObject("AliFlowCommonHistLYZ1PROD"));
457  pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
458  (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ1PROD"));
459  //Get the histograms from the output list
460  for(Int_t theta = 0;theta<iNtheta;theta++){
461  TString name = "AliFlowLYZHist1_";
462  name += theta;
463  pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
464  (outputListHistos->FindObject(name));
465  }
466  pHistVtheta = dynamic_cast<TH1D*>
467  (outputListHistos->FindObject("First_Flow_Vtheta_LYZPROD"));
468 
469  pHistR0theta = dynamic_cast<TH1D*>
470  (outputListHistos->FindObject("First_Flow_r0theta_LYZPROD"));
471 
472  pHistQsumforChi = dynamic_cast<TH1F*>
473  (outputListHistos->FindObject("Flow_QsumforChi_LYZPROD"));
474 
475  //Set the histogram pointers and call Finish()
476  if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
477  pHistVtheta && pHistR0theta && pHistQsumforChi ) {
478  this->SetCommonHists(pCommonHist);
479  this->SetCommonHistsRes(pCommonHistResults);
480  this->SetHist1(pLYZHist1);
481  this->SetHistVtheta(pHistVtheta);
482  this->SetHistR0theta(pHistR0theta);
483  this->SetHistQsumforChi(pHistQsumforChi);
484  } else {
485  cout<<"WARNING: Histograms needed to run Finish() firstrun (PROD) are not accessable!"<<endl;
486  }
487  }
488  }
489  else { //second run
490  //Get the common histograms from the output list
491  if (GetUseSum()){
492  pCommonHist = dynamic_cast<AliFlowCommonHist*>
493  (outputListHistos->FindObject("AliFlowCommonHistLYZ2SUM"));
494  pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
495  (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ2SUM"));
496 
497  pHistR0theta = dynamic_cast<TH1D*>
498  (outputListHistos->FindObject("First_Flow_r0theta_LYZSUM"));
499 
500  pHistQsumforChi = dynamic_cast<TH1F*>
501  (outputListHistos->FindObject("Flow_QsumforChi_LYZSUM"));
502 
503  //Get the histograms from the output list
504  for(Int_t theta = 0;theta<iNtheta;theta++){
505  TString nameRP = "AliFlowLYZHist2RP_";
506  nameRP += theta;
507  pLYZHist2RP[theta] = dynamic_cast<AliFlowLYZHist2*>
508  (outputListHistos->FindObject(nameRP));
509  TString namePOI = "AliFlowLYZHist2POI_";
510  namePOI += theta;
511  pLYZHist2POI[theta] = dynamic_cast<AliFlowLYZHist2*>
512  (outputListHistos->FindObject(namePOI));
513  }
514  pHistProReDenom = dynamic_cast<TProfile*>
515  (outputListHistos->FindObject("Second_FlowPro_ReDenom_LYZSUM"));
516  pHistProImDenom = dynamic_cast<TProfile*>
517  (outputListHistos->FindObject("Second_FlowPro_ImDenom_LYZSUM"));
518 
519  TH1D* pHistReDtheta = dynamic_cast<TH1D*>
520  (outputListHistos->FindObject("Second_Flow_ReDtheta_LYZSUM"));
521  TH1D* pHistImDtheta = dynamic_cast<TH1D*>
522  (outputListHistos->FindObject("Second_Flow_ImDtheta_LYZSUM"));
523 
524  pHistProVetaRP = dynamic_cast<TProfile*>
525  (outputListHistos->FindObject("Second_FlowPro_VetaRP_LYZSUM"));
526  pHistProVetaPOI = dynamic_cast<TProfile*>
527  (outputListHistos->FindObject("Second_FlowPro_VetaPOI_LYZSUM"));
528  pHistProVPtRP = dynamic_cast<TProfile*>
529  (outputListHistos->FindObject("Second_FlowPro_VPtRP_LYZSUM"));
530  pHistProVPtPOI = dynamic_cast<TProfile*>
531  (outputListHistos->FindObject("Second_FlowPro_VPtPOI_LYZSUM"));
532 
533 
534  //Set the histogram pointers and call Finish()
535  if (pCommonHist && pCommonHistResults &&
536  pLYZHist2RP[0] && pLYZHist2POI[0] &&
537  pHistR0theta &&
538  pHistProReDenom && pHistProImDenom &&
539  pHistReDtheta && pHistImDtheta &&
540  pHistProVetaRP && pHistProVetaPOI &&
541  pHistProVPtRP && pHistProVPtPOI &&
542  pHistQsumforChi) {
543  this->SetCommonHists(pCommonHist);
544  this->SetCommonHistsRes(pCommonHistResults);
545  this->SetHist2RP(pLYZHist2RP);
546  this->SetHist2POI(pLYZHist2POI);
547  this->SetHistR0theta(pHistR0theta);
548  this->SetHistProReDenom(pHistProReDenom);
549  this->SetHistProImDenom(pHistProImDenom);
550  this->SetHistReDtheta(pHistReDtheta);
551  this->SetHistImDtheta(pHistImDtheta);
552  this->SetHistProVetaRP(pHistProVetaRP);
553  this->SetHistProVetaPOI(pHistProVetaPOI);
554  this->SetHistProVPtRP(pHistProVPtRP);
555  this->SetHistProVPtPOI(pHistProVPtPOI);
556  this->SetHistQsumforChi(pHistQsumforChi);
557  }
558  else {
559  cout<<"WARNING: Histograms needed to run Finish() secondrun (SUM) are not accessable!"<<endl;
560  }
561  }
562  else {
563  pCommonHist = dynamic_cast<AliFlowCommonHist*>
564  (outputListHistos->FindObject("AliFlowCommonHistLYZ2PROD"));
565  pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
566  (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ2PROD"));
567 
568 
569  pHistR0theta = dynamic_cast<TH1D*>
570  (outputListHistos->FindObject("First_Flow_r0theta_LYZPROD"));
571 
572  pHistQsumforChi = dynamic_cast<TH1F*>
573  (outputListHistos->FindObject("Flow_QsumforChi_LYZPROD"));
574 
575  //Get the histograms from the output list
576  for(Int_t theta = 0;theta<iNtheta;theta++){
577  TString nameRP = "AliFlowLYZHist2RP_";
578  nameRP += theta;
579  pLYZHist2RP[theta] = dynamic_cast<AliFlowLYZHist2*>
580  (outputListHistos->FindObject(nameRP));
581  TString namePOI = "AliFlowLYZHist2POI_";
582  namePOI += theta;
583  pLYZHist2POI[theta] = dynamic_cast<AliFlowLYZHist2*>
584  (outputListHistos->FindObject(namePOI));
585  }
586  pHistProReDenom = dynamic_cast<TProfile*>
587  (outputListHistos->FindObject("Second_FlowPro_ReDenom_LYZPROD"));
588  pHistProImDenom = dynamic_cast<TProfile*>
589  (outputListHistos->FindObject("Second_FlowPro_ImDenom_LYZPROD"));
590 
591  pHistProVetaRP = dynamic_cast<TProfile*>
592  (outputListHistos->FindObject("Second_FlowPro_VetaRP_LYZPROD"));
593  pHistProVetaPOI = dynamic_cast<TProfile*>
594  (outputListHistos->FindObject("Second_FlowPro_VetaPOI_LYZPROD"));
595  pHistProVPtRP = dynamic_cast<TProfile*>
596  (outputListHistos->FindObject("Second_FlowPro_VPtRP_LYZPROD"));
597  pHistProVPtPOI = dynamic_cast<TProfile*>
598  (outputListHistos->FindObject("Second_FlowPro_VPtPOI_LYZPROD"));
599 
600  //Set the histogram pointers and call Finish()
601  if (pCommonHist && pCommonHistResults && pLYZHist2RP[0] && pLYZHist2POI[0] &&
602  pHistR0theta && pHistProReDenom && pHistProImDenom && pHistProVetaRP &&
603  pHistProVetaPOI && pHistProVPtRP && pHistProVPtPOI && pHistQsumforChi) {
604  this->SetCommonHists(pCommonHist);
605  this->SetCommonHistsRes(pCommonHistResults);
606  this->SetHist2RP(pLYZHist2RP);
607  this->SetHist2POI(pLYZHist2POI);
608  this->SetHistR0theta(pHistR0theta);
609  this->SetHistProReDenom(pHistProReDenom);
610  this->SetHistProImDenom(pHistProImDenom);
611  this->SetHistProVetaRP(pHistProVetaRP);
612  this->SetHistProVetaPOI(pHistProVetaPOI);
613  this->SetHistProVPtRP(pHistProVPtRP);
614  this->SetHistProVPtPOI(pHistProVPtPOI);
615  this->SetHistQsumforChi(pHistQsumforChi);
616  }
617  else {
618  cout<<"WARNING: Histograms needed to run Finish() secondrun (PROD) are not accessable!"<<endl;
619  }
620  }
621  } //secondrun
622  //outputListHistos->Print();
623  delete [] pLYZHist1;
624  delete [] pLYZHist2RP;
625  delete [] pLYZHist2POI;
626  } //listhistos
627  else {
628  cout << "histogram list pointer is empty in method AliFlowAnalysisWithLeeYangZeros::GetOutputHistograms() " << endl;}
629 }
630 
631  //-----------------------------------------------------------------------
633 {
634  //finish method
635  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
636 
637  //define variables for both runs
638  Double_t dJ01 = 2.405;
640 
641  //set the event number
642  SetEventNumber((int)fCommonHists->GetHistQ()->GetEntries());
643 
644  //Get multiplicity for RP selection
645  Double_t dMultRP = fCommonHists->GetHistMultRP()->GetMean();
646  if (fDebug) cout<<"The average multiplicity is "<<dMultRP<<endl;
647 
648  //set the sum of Q vectors
649  fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
650  SetQ2sum(fHistQsumforChi->GetBinContent(3));
651 
652  if (fFirstRun){
653 
654  //define variables for the first run
655  Double_t dR0 = 0;
656  Double_t dVtheta = 0;
657  Double_t dv = 0;
658  Double_t dV = 0;
659 
660  //reset histograms in case of merged output files
661  fHistR0theta->Reset();
662  fHistVtheta->Reset();
663 
664  for (Int_t theta=0;theta<iNtheta;theta++)
665  {
666  //get the first minimum r0
667  fHist1[theta]->FillGtheta();
668  dR0 = fHist1[theta]->GetR0();
669 
670  //calculate integrated flow
671  if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) { dVtheta = dJ01/dR0; }
672  else { cout<<"r0 is not found! Leaving LYZ analysis."<<endl; return kFALSE; }
673 
674  //for estimating systematic error resulting from d0
675  Double_t dBinsize =0.;
678  Double_t dVplus = -1.;
679  Double_t dVmin = -1.;
680  if (!TMath::AreEqualAbs(dR0+dBinsize, 0., 1e-100)) {dVplus = dJ01/(dR0+dBinsize);}
681  if (!TMath::AreEqualAbs(dR0-dBinsize, 0., 1e-100)) {dVmin = dJ01/(dR0-dBinsize);}
682  //convert V to v
683  Double_t dvplus = -1.;
684  Double_t dvmin= -1.;
685  if (fUseSum){
686  //for SUM: V=v because the Q-vector is scaled by 1/M
687  dv = dVtheta;
688  dvplus = dVplus;
689  dvmin = dVmin; }
690  else {
691  //for PRODUCT: v=V/M
692  if (!TMath::AreEqualAbs(dMultRP, 0., 1e-100)){
693  dv = dVtheta/dMultRP;
694  dvplus = dVplus/dMultRP;
695  dvmin = dVmin/dMultRP; }}
696 
697  if (fDebug) cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
698 
699  //fill the histograms
700  fHistR0theta->SetBinContent(theta+1,dR0);
701  fHistR0theta->SetBinError(theta+1,0.0);
702  fHistVtheta->SetBinContent(theta+1,dVtheta);
703  fHistVtheta->SetBinError(theta+1,0.0);
704  //get average value of fVtheta = fV
705  dV += dVtheta;
706  } //end of loop over theta
707 
708  //get average value of fVtheta = fV
709  dV /=iNtheta;
710  if (!fUseSum) { if (dMultRP!=0.){dV /=dMultRP;}} //scale with multiplicity for PRODUCT
711 
712  //sigma2 and chi
713  Double_t dSigma2 = 0;
714  Double_t dChi= 0;
715  if (fEventNumber!=0) {
716  *fQsum /= fEventNumber;
717  //cout<<"fQsum is "<<fQsum->X()<<" "<<fQsum->Y()<<endl;
718  fQ2sum /= fEventNumber;
719  //cout<<"fQ2sum is "<<fQ2sum<<endl;
720  dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
721  //cout<<"dSigma2 is "<<dSigma2<<endl;
722  if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
723  else dChi = -1.;
724  fCommonHistsRes->FillChi(dChi);
725 
726  cout<<"*************************************"<<endl;
727  cout<<"*************************************"<<endl;
728  cout<<" Integrated flow from "<<endl;
729  if (fUseSum) {
730  cout<<" Lee-Yang Zeroes SUM "<<endl;}
731  else {
732  cout<<" Lee-Yang Zeroes PRODUCT "<<endl;}
733  cout<<endl;
734  cout<<"Chi = "<<dChi<<endl;
735  //cout<<endl;
736  }
737 
738  // recalculate statistical errors on integrated flow
739  //combining 5 theta angles to 1 relative error BP eq. 89
740  Double_t dRelErr2comb = 0.;
741  Int_t iEvts = fEventNumber;
742  if (iEvts!=0) {
743  for (Int_t theta=0;theta<iNtheta;theta++){
744  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
745  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
746  TMath::Cos(dTheta));
747  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
748  TMath::Cos(dTheta));
749  dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
750  TMath::BesselJ1(dJ01)))*
751  (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
752  dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
753  }
754  dRelErr2comb /= iNtheta;
755  }
756  Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
757 
758  //copy content of profile into TH1D and add error
759  Double_t dv2pro = dV; //in the case that fv is equal to fV
760  Double_t dv2Err = dv2pro*dRelErrcomb ;
761  cout<<"dV = "<<dv2pro<<" +- "<<dv2Err<<endl;
762  cout<<endl;
763  cout<<"*************************************"<<endl;
764  cout<<"*************************************"<<endl;
765  fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
766 
767 
768  if (fDebug) cout<<"****histograms filled****"<<endl;
769 
770  return kTRUE;
771  } //firstrun
772 
773 
774  else { //second run to calculate differential flow
775 
776  //declare variables for the second run
777  TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
778  Int_t m = 1;
779  TComplex i = TComplex::I();
780  Double_t dBesselRatio[3] = {1., 1.202, 2.69};
783 
784  Double_t dR0 = 0.;
785  Double_t dVtheta = 0.;
786  Double_t dV = 0.;
787  Double_t dReDenom = 0.;
788  Double_t dImDenom = 0.;
789 
790  //reset histograms in case of merged output files
791  if (fUseSum) {
792  fHistReDtheta->Reset();
793  fHistImDtheta->Reset();
794  }
795  fHistProVetaRP ->Reset();
796  fHistProVetaPOI->Reset();
797  fHistProVPtRP ->Reset();
798  fHistProVPtPOI ->Reset();
799 
800  //scale fHistR0theta by the number of merged files to undo the merging
801  if (!fHistR0theta) { cout<<"Hist pointer R0theta in file does not exist"<<endl; }
802  else {
803  Int_t iEntries = (int)fHistR0theta->GetEntries();
804  if (iEntries > iNtheta){
805  //for each individual file fHistR0theta has iNtheta entries
806  Int_t iFiles = iEntries/iNtheta;
807  cout<<iFiles<<" files were merged!"<<endl;
808  fHistR0theta->Scale(1./iFiles);
809  }
810 
811  //loop over theta
812  for (Int_t theta=0;theta<iNtheta;theta++) {
813  dR0 = fHistR0theta->GetBinContent(theta+1); //histogram starts at bin 1
814  if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
815  if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) dVtheta = dJ01/dR0;
816  dV += dVtheta;
817  // BP Eq. 9 -> Vn^theta = j01/r0^theta
818  if (!fHistProReDenom || !fHistProImDenom) {
819  cout << "Hist pointer fDenom in file does not exist" <<endl;
820  cout<< "Leaving LYZ second pass analysis!" <<endl;
821  return kFALSE;
822  } else {
823  dReDenom = fHistProReDenom->GetBinContent(theta+1);
824  dImDenom = fHistProImDenom->GetBinContent(theta+1);
825  cDenom(dReDenom,dImDenom);
826 
827  //for new method and use by others (only with the sum generating function):
828  if (fUseSum) {
829  cDtheta = dR0*cDenom/dJ01;
830  fHistReDtheta->SetBinContent(theta+1,cDtheta.Re());
831  fHistReDtheta->SetBinError(theta+1,0.0);
832  fHistImDtheta->SetBinContent(theta+1,cDtheta.Im());
833  fHistImDtheta->SetBinError(theta+1,0.0);
834  }
835 
836  cDenom *= TComplex::Power(i, m-1);
837 
838  //v as a function of eta for RP selection
839  for (Int_t be=1;be<=iNbinsEta;be++) {
840  Double_t dReRatioRP = 0.0;
841  Double_t dEtaRP = fHist2RP[theta]->GetBinCenter(be);
842  cNumerRP = fHist2RP[theta]->GetNumerEta(be);
843  if (cNumerRP.Rho()==0) {
844  if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
845  }
846  else if (TMath::AreEqualAbs(cDenom.Rho(), 0, 1e-100)) {
847  if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
848  }
849  else {
850  dReRatioRP = (cNumerRP/cDenom).Re();
851  }
852  Double_t dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
853  fHistProVetaRP->Fill(dEtaRP,dVetaRP);
854  } //loop over bins be
855 
856 
857  //v as a function of eta for POI selection
858  for (Int_t be=1;be<=iNbinsEta;be++) {
859  Double_t dReRatioPOI = 0.0;
860  Double_t dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
861  cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
862  if (cNumerPOI.Rho()==0) {
863  if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
864  }
865  else if (cDenom.Rho()==0) {
866  if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
867  }
868  else {
869  dReRatioPOI = (cNumerPOI/cDenom).Re();
870  }
871  Double_t dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
872  fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
873  } //loop over bins be
874 
875 
876  //v as a function of Pt for RP selection
877  for (Int_t bp=1;bp<=iNbinsPt;bp++) {
878  Double_t dReRatioRP = 0.0;
879  Double_t dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
880  cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
881  if (cNumerRP.Rho()==0) {
882  if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
883  }
884  else if (cDenom.Rho()==0) {
885  if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
886  }
887  else {
888  dReRatioRP = (cNumerRP/cDenom).Re();
889  }
890  Double_t dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
891  fHistProVPtRP->Fill(dPtRP,dVPtRP);
892  } //loop over bins bp
893 
894 
895  //v as a function of Pt for POI selection
896  for (Int_t bp=1;bp<=iNbinsPt;bp++) {
897  Double_t dReRatioPOI = 0.0;
898  Double_t dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
899  cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
900  if (cNumerPOI.Rho()==0) {
901  if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
902  }
903  else if (cDenom.Rho()==0) {
904  if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
905  }
906  else {
907  dReRatioPOI = (cNumerPOI/cDenom).Re();
908  }
909  Double_t dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
910  fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
911  } //loop over bins bp
912 
913  }
914  }
915 
916  }//end of loop over theta
917 
918  //calculate the average of fVtheta = fV
919  dV /= iNtheta;
920  if (!fUseSum) { if (dMultRP!=0.) { dV /=dMultRP; }} //scale by the multiplicity for PRODUCT
921  if (TMath::AreEqualAbs(dV, 0., 1e-100)) { cout<<"dV = 0! Leaving LYZ analysis."<<endl; return kFALSE; }
922 
923  //sigma2 and chi (for statistical error calculations)
924  Double_t dSigma2 = 0;
925  Double_t dChi= 0;
926  if (fEventNumber!=0) {
927  *fQsum /= fEventNumber;
928  //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
929  //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
930  fQ2sum /= fEventNumber;
931  //cout<<"fQ2sum = "<<fQ2sum<<endl;
932  dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
933  if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
934  else dChi = -1.;
935  fCommonHistsRes->FillChi(dChi);
936 
937  // recalculate statistical errors on integrated flow
938  //combining 5 theta angles to 1 relative error BP eq. 89
939  Double_t dRelErr2comb = 0.;
940  Int_t iEvts = fEventNumber;
941  if (iEvts!=0) {
942  for (Int_t theta=0;theta<iNtheta;theta++){
943  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
944  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
945  TMath::Cos(dTheta));
946  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
947  TMath::Cos(dTheta));
948  dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
949  TMath::BesselJ1(dJ01)))*
950  (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
951  dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
952  }
953  dRelErr2comb /= iNtheta;
954  }
955  Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
956  Double_t dVErr = dV*dRelErrcomb ;
958 
959  cout<<"*************************************"<<endl;
960  cout<<"*************************************"<<endl;
961  cout<<" Integrated flow from "<<endl;
962  if (fUseSum) {
963  cout<<" Lee-Yang Zeroes SUM "<<endl;}
964  else {
965  cout<<" Lee-Yang Zeroes PRODUCT "<<endl;}
966  cout<<endl;
967  cout<<"Chi = "<<dChi<<endl;
968  cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
969  //cout<<endl;
970  }
971 
972  //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
973 
974  //v as a function of eta for RP selection
975  for(Int_t b=0;b<iNbinsEta;b++) {
976  Double_t dv2pro = fHistProVetaRP->GetBinContent(b);
978  Double_t dErrdifcomb = 0.; //set error to zero
979  Double_t dErr2difcomb = 0.; //set error to zero
980  //calculate error
981  if (!TMath::AreEqualAbs(dNprime, 0., 1e-100)) {
982  for (Int_t theta=0;theta<iNtheta;theta++) {
983  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
984  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
985  TMath::Cos(dTheta));
986  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
987  TMath::Cos(dTheta));
988  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
989  TMath::BesselJ1(dJ01)))*
990  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
991  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
992  } //loop over theta
993  }
994 
995  if (!TMath::AreEqualAbs(dErr2difcomb, 0., 1e-100)) {
996  dErr2difcomb /= iNtheta;
997  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
998  }
999  else {dErrdifcomb = 0.;}
1000  //fill TH1D
1001  fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
1002  } //loop over bins b
1003 
1004 
1005  //v as a function of eta for POI selection
1006  for(Int_t b=0;b<iNbinsEta;b++) {
1007  Double_t dv2pro = fHistProVetaPOI->GetBinContent(b);
1009  Double_t dErrdifcomb = 0.; //set error to zero
1010  Double_t dErr2difcomb = 0.; //set error to zero
1011  //calculate error
1012  if (dNprime!=0.) {
1013  for (Int_t theta=0;theta<iNtheta;theta++) {
1014  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1015  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1016  TMath::Cos(dTheta));
1017  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1018  TMath::Cos(dTheta));
1019  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1020  TMath::BesselJ1(dJ01)))*
1021  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1022  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1023  } //loop over theta
1024  }
1025 
1026  if (dErr2difcomb!=0.) {
1027  dErr2difcomb /= iNtheta;
1028  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1029  }
1030  else {dErrdifcomb = 0.;}
1031  //fill TH1D
1032  fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
1033  } //loop over bins b
1034 
1035 
1036 
1037  //v as a function of Pt for RP selection
1038  TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
1039  Double_t dVRP = 0.;
1040  Double_t dSum = 0.;
1041  Double_t dErrV =0.;
1042  for(Int_t b=0;b<iNbinsPt;b++) {
1043  Double_t dv2pro = fHistProVPtRP->GetBinContent(b);
1044  Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
1045  Double_t dErrdifcomb = 0.; //set error to zero
1046  Double_t dErr2difcomb = 0.; //set error to zero
1047  //calculate error
1048  if (dNprime!=0.) {
1049  for (Int_t theta=0;theta<iNtheta;theta++) {
1050  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1051  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1052  TMath::Cos(dTheta));
1053  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1054  TMath::Cos(dTheta));
1055  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1056  TMath::BesselJ1(dJ01)))*
1057  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1058  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1059  } //loop over theta
1060  }
1061 
1062  if (dErr2difcomb!=0.) {
1063  dErr2difcomb /= iNtheta;
1064  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1065  //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1066  }
1067  else {dErrdifcomb = 0.;}
1068 
1069  //fill TH1D
1070  fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
1071  //calculate integrated flow for RP selection
1072  if (fHistPtRP){
1073  Double_t dYieldPt = fHistPtRP->GetBinContent(b);
1074  dVRP += dv2pro*dYieldPt;
1075  dSum +=dYieldPt;
1076  dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
1077  } else { cout<<"fHistPtRP is NULL"<<endl; }
1078  } //loop over bins b
1079 
1080  if (!TMath::AreEqualAbs(dSum, 0., 1e-100)) {
1081  dVRP /= dSum; //the pt distribution should be normalised
1082  dErrV /= (dSum*dSum);
1083  dErrV = TMath::Sqrt(dErrV);
1084  }
1085  cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
1086  //cout<<endl;
1088 
1089 
1090  //v as a function of Pt for POI selection
1091  TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
1092  Double_t dVPOI = 0.;
1093  dSum = 0.;
1094  dErrV =0.;
1095 
1096  for(Int_t b=0;b<iNbinsPt;b++) {
1097  Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
1099  Double_t dErrdifcomb = 0.; //set error to zero
1100  Double_t dErr2difcomb = 0.; //set error to zero
1101  //calculate error
1102  if (dNprime!=0.) {
1103  for (Int_t theta=0;theta<iNtheta;theta++) {
1104  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1105  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1106  TMath::Cos(dTheta));
1107  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1108  TMath::Cos(dTheta));
1109  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1110  TMath::BesselJ1(dJ01)))*
1111  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1112  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1113  } //loop over theta
1114  }
1115 
1116  if (dErr2difcomb!=0.) {
1117  dErr2difcomb /= iNtheta;
1118  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1119  //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1120  }
1121  else {dErrdifcomb = 0.;}
1122 
1123  //fill TH1D
1124  fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
1125  //calculate integrated flow for POI selection
1126  if (fHistPtPOI){
1127  Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
1128  dVPOI += dv2pro*dYieldPt;
1129  dSum +=dYieldPt;
1130  dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
1131  } else { cout<<"fHistPtPOI is NULL"<<endl; }
1132  } //loop over bins b
1133 
1134  if (dSum != 0.) {
1135  dVPOI /= dSum; //the pt distribution should be normalised
1136  dErrV /= (dSum*dSum);
1137  dErrV = TMath::Sqrt(dErrV);
1138  }
1139  cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
1140  cout<<endl;
1141  cout<<"*************************************"<<endl;
1142  cout<<"*************************************"<<endl;
1143  fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
1144 
1145  } //secondrun
1146 
1147  //cout<<"----LYZ analysis finished....----"<<endl<<endl;
1148 
1149  return kTRUE;
1150 }
1151 
1152 
1153 //-----------------------------------------------------------------------
1154 
1156 {
1157  // Get event quantities from AliFlowEvent for all particles
1158 
1159  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
1160 
1161  if (!anEvent){
1162  cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1163  return kFALSE;
1164  }
1165 
1166  //define variables
1167  TComplex cExpo, cGtheta, cGthetaNew, cZ;
1170 
1171  //calculate flow
1172  Double_t dOrder = 2.;
1173 
1174  //get the Q vector
1175  AliFlowVector vQ = anEvent->GetQ();
1176  //weight by the multiplicity
1177  Double_t dQX = 0;
1178  Double_t dQY = 0;
1179  if (!TMath::AreEqualAbs(vQ.GetMult(), 0., 1e-100)) {
1180  dQX = vQ.X()/vQ.GetMult();
1181  dQY = vQ.Y()/vQ.GetMult();
1182  }
1183  vQ.Set(dQX,dQY);
1184 
1185  //for chi calculation:
1186  *fQsum += vQ;
1187  fHistQsumforChi->SetBinContent(1,fQsum->X());
1188  fHistQsumforChi->SetBinContent(2,fQsum->Y());
1189  fQ2sum += vQ.Mod2();
1190  fHistQsumforChi->SetBinContent(3,fQ2sum);
1191 
1192  for (Int_t theta=0;theta<iNtheta;theta++) {
1193  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1194 
1195  //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1196  Double_t dQtheta = GetQtheta(vQ, dTheta);
1197 
1198  for (Int_t bin=1;bin<=iNbins;bin++) {
1199  Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED???
1200  if (fUseSum) {
1201  //calculate the sum generating function
1202  cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
1203  cGtheta = TComplex::Exp(cExpo);
1204  }
1205  else {
1206  //calculate the product generating function
1207  cGtheta = GetGrtheta(anEvent, dR, dTheta);
1208  if (cGtheta.Rho2() > 100.) break;
1209  }
1210  //fill real and imaginary part of cGtheta
1211  fHist1[theta]->Fill(dR,cGtheta);
1212  } //loop over bins
1213  } //loop over theta
1214 
1215  return kTRUE;
1216 
1217 }
1218 
1219  //-----------------------------------------------------------------------
1221 {
1222  //for differential flow
1223 
1224  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
1225 
1226  if (!anEvent){
1227  cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1228  return kFALSE;
1229  }
1230 
1231  //define variables
1232  TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
1233  Double_t dR0 = 0.;
1234  Double_t dCosTermRP = 0.;
1235  Double_t dCosTermPOI = 0.;
1236  Double_t m = 1.;
1237  Double_t dOrder = 2.;
1239 
1240  //get the Q vector
1241  AliFlowVector vQ = anEvent->GetQ();
1242  //weight by the multiplicity
1243  Double_t dQX = 0.;
1244  Double_t dQY = 0.;
1245  if (vQ.GetMult() != 0) {
1246  dQX = vQ.X()/vQ.GetMult();
1247  dQY = vQ.Y()/vQ.GetMult();
1248  }
1249  vQ.Set(dQX,dQY);
1250 
1251  //for chi calculation:
1252  *fQsum += vQ;
1253  fHistQsumforChi->SetBinContent(1,fQsum->X());
1254  fHistQsumforChi->SetBinContent(2,fQsum->Y());
1255  fQ2sum += vQ.Mod2();
1256  fHistQsumforChi->SetBinContent(3,fQ2sum);
1257 
1258  for (Int_t theta=0;theta<iNtheta;theta++)
1259  {
1260  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1261 
1262  //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1263  Double_t dQtheta = GetQtheta(vQ, dTheta);
1264 
1265  //denominator for differential v
1266  if (fHistR0theta) {
1267  dR0 = fHistR0theta->GetBinContent(theta+1);
1268  }
1269  else { cout <<"pointer fHistR0theta does not exist" << endl;
1270  }
1271 
1272  if (fUseSum) //sum generating function
1273  {
1274  cExpo(0.,dR0*dQtheta);
1275  cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
1276  //loop over tracks in event
1277  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1278  for (Int_t i=0;i<iNumberOfTracks;i++) {
1279  AliFlowTrackSimple* pTrack = anEvent->GetTrack(i);
1280  if (pTrack) {
1281  Double_t dEta = pTrack->Eta();
1282  Double_t dPt = pTrack->Pt();
1283  Double_t dPhi = pTrack->Phi();
1284  if (pTrack->InRPSelection()) { // RP selection
1285  dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
1286  cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
1287  if (cNumerRP.Rho()==0) { cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
1288  if (fDebug) { cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;}
1289  if (fHist2RP[theta]) {
1290  fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1291  }
1292  }
1293  if (pTrack->InPOISelection()) { //POI selection
1294  dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
1295  cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
1296  if (cNumerPOI.Rho()==0) { cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
1297  if (fDebug) { cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;}
1298  if (fHist2POI[theta]) {
1299  fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1300  }
1301  }
1302  } //if track
1303  else {cerr << "no particle!!!"<<endl;}
1304  } //loop over tracks
1305  } //sum
1306  else { //product generating function
1307  cDenom = GetDiffFlow(anEvent, dR0, theta);
1308  }//product
1309 
1311  fHistProReDenom->Fill(theta,cDenom.Re()); //fill the real part of fDenom
1312  fHistProImDenom->Fill(theta,cDenom.Im()); //fill the imaginary part of fDenom
1313  }
1314  else { cout << "Pointers to cDenom mising" << endl;}
1315 
1316  }//end of loop over theta
1317 
1318  return kTRUE;
1319 
1320 }
1321  //-----------------------------------------------------------------------
1323 {
1324  //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
1325  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
1326 
1327  Double_t dQtheta = 0.;
1328  Double_t dOrder = 2.;
1329 
1330  dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
1331 
1332  return dQtheta;
1333 
1334 }
1335 
1336 
1337 //-----------------------------------------------------------------------
1339 {
1340  // Product Generating Function for LeeYangZeros method
1341  // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1342 
1343  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1344 
1345 
1346  TComplex cG = TComplex::One();
1347  Double_t dOrder = 2.;
1348  Double_t dWgt = 1.;
1349  //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1350 
1351  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1352 
1353  for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1354  {
1355  AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1356  if (pTrack){
1357  if (pTrack->InRPSelection()) {
1358  Double_t dPhi = pTrack->Phi();
1359  Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
1360  TComplex cGi(1., dGIm);
1361  cG *= cGi; //product over all tracks
1362  }
1363  }
1364  else {cerr << "no particle pointer !!!"<<endl;}
1365  }//loop over tracks
1366 
1367  return cG;
1368 
1369 }
1370 
1371 
1372 //-----------------------------------------------------------------------
1374 {
1375  // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1376  // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1377  // Also for v1 mixed harmonics: DF Eq. 5
1378  // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1379 
1380  if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1381 
1382  TComplex cG = TComplex::One();
1383  TComplex cdGr0(0.,0.);
1384  Double_t dOrder = 2.;
1385  Double_t dWgt = 1.;
1386  //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1387 
1388  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1389 
1391  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1392 
1393  //for the denominator (use all RP selected particles)
1394  for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1395  {
1396  AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1397  if (pTrack){
1398  if (pTrack->InRPSelection()) {
1399  Double_t dPhi = pTrack->Phi();
1400  Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
1401  //GetGr0theta
1402  Double_t dGIm = aR0 * dCosTerm;
1403  TComplex cGi(1., dGIm);
1404  TComplex cCosTermComplex(1., aR0*dCosTerm);
1405  cG *= cGi; //product over all tracks
1406  //GetdGr0theta
1407  cdGr0 +=(dCosTerm / cCosTermComplex); //sum over all tracks
1408  }
1409  } //if particle
1410  else {cerr << "no particle!!!"<<endl;}
1411  }//loop over tracks
1412 
1413  //for the numerator
1414  for (Int_t i=0;i<iNumberOfTracks;i++)
1415  {
1416  AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1417  if (pTrack){
1418  Double_t dEta = pTrack->Eta();
1419  Double_t dPt = pTrack->Pt();
1420  Double_t dPhi = pTrack->Phi();
1421  Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1422  TComplex cCosTermComplex(1.,aR0*dCosTerm);
1423  //RP selection
1424  if (pTrack->InRPSelection()) {
1425  TComplex cNumerRP = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
1426  fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1427  }
1428  //POI selection
1429  if (pTrack->InPOISelection()) {
1430  TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
1431  fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1432  }
1433  } //if particle
1434  else {cerr << "no particle pointer!!!"<<endl;}
1435  }//loop over tracks
1436 
1437  TComplex cDenom = cG*cdGr0;
1438  return cDenom;
1439 
1440 }
1441 
1442 //-----------------------------------------------------------------------
1443 
void Fill(Double_t f, TComplex c)
void SetHist2POI(AliFlowLYZHist2 *const aLYZHist2POI[])
TComplex GetGrtheta(AliFlowEventSimple *anEvent, Double_t aR, Double_t aTheta)
double Double_t
Definition: External.C:58
virtual AliFlowVector GetQ(Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
void SetHistProVetaRP(TProfile *const aHistProVetaRP)
void SetHistProVetaPOI(TProfile *const aHistProVetaPOI)
Double_t GetEntriesInEtaBinPOI(Int_t iBin)
Int_t GetNtheta() const
AliFlowTrackSimple * GetTrack(Int_t i)
void SetHistProReDenom(TProfile *const aHistProReDenom)
Double_t GetQtheta(AliFlowVector aQ, Double_t aTheta)
Double_t GetEntriesInPtBinPOI(Int_t iBin)
Double_t GetMult() const
Definition: AliFlowVector.h:46
Double_t GetBinCenter(Int_t i)
Bool_t FillDifferentialFlowPtRP(Int_t aBin, Double_t av, Double_t anError)
void SetCommonHists(AliFlowCommonHist *const aCommonHist)
void SetHistR0theta(TH1D *const aHistR0theta)
Bool_t InRPSelection() const
Bool_t FillControlHistograms(AliFlowEventSimple *anEvent, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
TComplex GetNumerPt(Int_t i)
static AliFlowLYZConstants * GetMaster()
Double_t Phi() const
Bool_t FillIntegratedFlowRP(Double_t aV, Double_t anError)
int Int_t
Definition: External.C:63
Double_t GetBinCenter(Int_t i)
void SetHistProVPtPOI(TProfile *const aHistProVPtPOI)
void Fill(Double_t d1, Double_t d2, TComplex c)
Double_t GetMaxPROD() const
void SetCommonHistsRes(AliFlowCommonHistResults *const aCommonHistResult)
Definition: External.C:212
void SetHistReDtheta(TH1D *const aHistReDtheta)
void SetHist2RP(AliFlowLYZHist2 *const aLYZHist2RP[])
Double_t GetEntriesInPtBinRP(Int_t iBin)
static AliFlowCommonConstants * GetMaster()
Double_t GetEntriesInEtaBinRP(Int_t iBin)
void SetHistVtheta(TH1D *const aHistVtheta)
TComplex GetDiffFlow(AliFlowEventSimple *anEvent, Double_t aR, Int_t theta)
Bool_t FillIntegratedFlow(Double_t aV, Double_t anError)
Bool_t FillDifferentialFlowEtaPOI(Int_t aBin, Double_t av, Double_t anError)
Bool_t FillIntegratedFlowPOI(Double_t aV, Double_t anError)
Double_t GetBinCenterPt(Int_t i)
void SetHistQsumforChi(TH1F *const aHistQsumforChi)
Bool_t Make(AliFlowEventSimple *anEvent)
Bool_t FillDifferentialFlowPtPOI(Int_t aBin, Double_t av, Double_t anError)
Double_t Pt() const
Bool_t FillFromFlowEvent(AliFlowEventSimple *anEvent)
TComplex GetNumerEta(Int_t i)
void GetOutputHistograms(TList *outputListHistos)
Bool_t InPOISelection(Int_t poiType=1) const
Double_t Eta() const
bool Bool_t
Definition: External.C:53
void SetHistProVPtRP(TProfile *const aHistProVPtRP)
void SetHistProImDenom(TProfile *const aHistProImDenom)
void SetHistImDtheta(TH1D *const aHistImDtheta)
Bool_t FillDifferentialFlowEtaRP(Int_t aBin, Double_t av, Double_t anError)
Bool_t SecondFillFromFlowEvent(AliFlowEventSimple *anEvent)
Double_t GetMaxSUM() const
void SetHist1(AliFlowLYZHist1 *const aLYZHist1[])
Int_t NumberOfTracks() const