bug fix
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithFittingQDistribution.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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  * integrated flow estimate by  *
18  *   fitting q-distribution     * 
19  *                              *
20  * author: Ante Bilandzic       * 
21  *          (anteb@nikhef.nl)   *
22  *                              *  
23  *  based on the macro written  *
24  *     by Sergei Voloshin       *
25  *******************************/  
26
27 #define AliFlowAnalysisWithFittingQDistribution_cxx
28
29 #include "Riostream.h"
30 #include "AliFlowCommonConstants.h"
31 #include "AliFlowCommonHist.h"
32 #include "AliFlowCommonHistResults.h"
33 #include "TChain.h"
34 #include "TFile.h"
35 #include "TList.h" 
36 #include "TF1.h"
37 #include "TLegend.h"
38 #include "TParticle.h"
39 #include "TProfile.h"
40 #include "AliFlowEventSimple.h"
41 #include "AliFlowTrackSimple.h"
42 #include "AliFlowAnalysisWithFittingQDistribution.h"
43
44
45 class TH1;
46 class TGraph;
47 class TPave;
48 class TLatex;
49 class TMarker;
50 class TObjArray;
51 class TList;
52 class TCanvas;
53 class TSystem;
54 class TROOT;
55 class AliFlowVector;
56 class TVector;
57
58 //================================================================================================================
59
60 ClassImp(AliFlowAnalysisWithFittingQDistribution)
61
62 AliFlowAnalysisWithFittingQDistribution::AliFlowAnalysisWithFittingQDistribution():  
63  fHistList(NULL),
64  fCommonHists(NULL),
65  fCommonHistsResults(NULL),
66  fnBinsPhi(0),
67  fPhiMin(0),
68  fPhiMax(0),
69  fPhiBinWidth(0),
70  fnBinsPt(0),
71  fPtMin(0),
72  fPtMax(0),
73  fPtBinWidth(0),
74  fnBinsEta(0),
75  fEtaMin(0),
76  fEtaMax(0),
77  fEtaBinWidth(0),
78  fHarmonic(2),
79  fAnalysisLabel(NULL),
80  fWeightsList(NULL),
81  fUsePhiWeights(kFALSE),
82  fUsePtWeights(kFALSE),
83  fUseEtaWeights(kFALSE),
84  fUseParticleWeights(NULL),
85  fPhiWeights(NULL),
86  fPtWeights(NULL),
87  fEtaWeights(NULL),
88  // fitting parameters with default values:
89  fFittingParameters(NULL), 
90  fTreshold(5),
91  fvStart(0.05),
92  fvMin(0.0),
93  fvMax(0.25),
94  fSigma2Start(0.75),
95  fSigma2Min(0.5), // (should be kept fixed at 0.5 according to theorists...)
96  fSigma2Max(2.5),
97  fPlotResults(kFALSE),
98  // rest:
99  fLegend(NULL)
100  {
101   // constructor 
102   
103   // base list to hold all output objects:
104   fHistList = new TList();
105   fHistList->SetName("cobjFQD");
106   fHistList->SetOwner(kTRUE);
107   
108   // analysis label;
109   fAnalysisLabel = new TString();
110  
111   // list to hold histograms with phi, pt and eta weights:      
112   fWeightsList = new TList();
113
114   // initialize all arrays:  
115   this->InitializeArrays();
116
117  } // end of constructor
118  
119  
120 //================================================================================================================
121  
122
123 AliFlowAnalysisWithFittingQDistribution::~AliFlowAnalysisWithFittingQDistribution()
124 {
125  // desctructor
126  delete fHistList; 
127 }
128
129
130 //================================================================================================================
131
132
133 void AliFlowAnalysisWithFittingQDistribution::Init()
134 {
135  // access constants and book everything
136  
137  // access constants:
138  this->AccessConstants();
139  
140  // booking:
141  this->BookCommonHistograms();
142  this->BookAndFillWeightsHistograms();
143  this->BookEverythingForDistributions();
144  
145  // store fitting parameters:  
146  this->StoreFittingParameters();
147  
148  // nest lists:
149  fWeightsList->SetName("Weights");
150  fWeightsList->SetOwner(kTRUE);   
151  fHistList->Add(fWeightsList);
152  
153  // set harmonic in common control histograms (to be improved (should I do this somewhere else?)):
154  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic); 
155  
156 } // end of void AliFlowAnalysisWithFittingQDistribution::Init()
157
158
159 //================================================================================================================
160
161
162 void AliFlowAnalysisWithFittingQDistribution::Make(AliFlowEventSimple* anEvent)
163 {
164  // loop over data
165  
166  // a) fill the common control histograms
167  // b) loop over data and calculate non-weighted and weighted Q-vector and sum of particle weights
168  // c) fill histograms for distribution
169  // d) reset e-b-e quantities
170  
171  // fill the common control histograms
172  fCommonHists->FillControlHistograms(anEvent); 
173  
174  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
175  Double_t dPt  = 0.; // transverse momentum
176  Double_t dEta = 0.; // pseudorapidity
177
178  Double_t wPhi = 1.; // phi weight
179  Double_t wPt  = 1.; // pt weight
180  Double_t wEta = 1.; // eta weight
181  
182  Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
183                                            // nRP   = # of particles used to determine the reaction plane;
184                                            // nPOI  = # of particles of interest for a detailed flow analysis;
185                                            // rest  = # of particles which are not niether RPs nor POIs.  
186  
187  Int_t n = fHarmonic; // shortcut for the harmonic
188  
189  Double_t dReQ[2] = {0.}; // real part of Q-vector [0=particle weights not used, 1=particle weights used]
190  Double_t dImQ[2] = {0.}; // imaginary part of Q-vector [0=particle weights not used, 1=particle weights used]
191  Double_t dSumOfParticleWeights[2]; // [0=particle weights not used, 1=particle weights used] 
192                                                                                                                                
193  AliFlowTrackSimple *aftsTrack = NULL;                                          
194                                            
195  for(Int_t i=0;i<nPrim;i++) 
196  { 
197   aftsTrack=anEvent->GetTrack(i);
198   if(aftsTrack)
199   {
200    if(!(aftsTrack->InRPSelection())) continue; // consider only tracks which are RPs
201      
202    dPhi = aftsTrack->Phi();
203    dPt  = aftsTrack->Pt();
204    dEta = aftsTrack->Eta();
205   
206    if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi weight for this particle:
207    {
208     wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
209    }
210    if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt weight for this particle:
211    {
212     wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); 
213    }            
214    if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta weight for this particle: 
215    {
216     wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); 
217    } 
218     
219    // calculate real and imaginary part of non-weighted and weighted Q-vector and sum of particle weights for this event:
220    for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++) // pW not used or used
221    {
222     // Q-vector:
223     dReQ[pW]+=pow(wPhi*wPt*wEta,pW)*TMath::Cos(n*dPhi); 
224     dImQ[pW]+=pow(wPhi*wPt*wEta,pW)*TMath::Sin(n*dPhi);
225     // sum of particle weights:
226     dSumOfParticleWeights[pW] += pow(wPhi*wPt*wEta,pW); // if pW = 0, this sum gives nRP
227    } 
228    
229   } // end of if(aftsTrack)
230  } // end of for(Int_t i=0;i<nPrim;i++)                                      
231                                            
232  // calculate q = Q\sqrt{sum of particle weights}:
233  // Remark: if particle weights are unit than sum of particle weights = multiplicity
234  Double_t q=0;                                          
235  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++) // pW not used or used
236  {
237   if(dSumOfParticleWeights[pW])
238   {
239    q = pow(dReQ[pW]*dReQ[pW]+dImQ[pW]*dImQ[pW],0.5)/pow(dSumOfParticleWeights[pW],0.5);
240    // fill histograms:
241    fqDistribution[pW]->Fill(q,1.);
242    fSumOfParticleWeights[pW]->Fill(dSumOfParticleWeights[pW],1.);
243   }
244  } 
245  
246  // reset e-b-e quantities:
247  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++) // pW not used or used
248  {
249   dReQ[pW] = 0.;
250   dImQ[pW] = 0.;
251   dSumOfParticleWeights[pW] = 0.;
252  }
253
254 }//end of Make()
255
256
257 //================================================================================================================
258
259
260 void AliFlowAnalysisWithFittingQDistribution::GetOutputHistograms(TList *outputListHistos) 
261 {
262  // get pointers to all output histograms (called before Finish()) 
263  
264  if(outputListHistos)   
265  {   
266   // 1.) common control histograms and common histograms for final results:
267   TString commonHistName = "AliFlowCommonHistFQD";
268   commonHistName += fAnalysisLabel->Data();
269   AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject(commonHistName.Data()));
270   if(commonHist) this->SetCommonHists(commonHist); 
271   
272   TString commonHistResName = "AliFlowCommonHistResultsFQD";
273   commonHistResName += fAnalysisLabel->Data();
274   AliFlowCommonHistResults *commonHistRes = dynamic_cast<AliFlowCommonHistResults*>
275                                             (outputListHistos->FindObject(commonHistResName.Data()));
276   if(commonHistRes) this->SetCommonHistsResults(commonHistRes); 
277   
278   // 2.) weights: 
279   TList *weightsList = dynamic_cast<TList*>(outputListHistos->FindObject("Weights"));
280   if(weightsList) this->SetWeightsList(weightsList);
281   Bool_t bUsePhiWeights = kFALSE;
282   Bool_t bUsePtWeights = kFALSE;
283   Bool_t bUseEtaWeights = kFALSE;
284   TString fUseParticleWeightsName = "fUseParticleWeightsFQD";
285   fUseParticleWeightsName += fAnalysisLabel->Data();
286   TProfile *useParticleWeights = dynamic_cast<TProfile*>(weightsList->FindObject(fUseParticleWeightsName.Data()));
287   if(useParticleWeights)
288   {
289    this->SetUseParticleWeights(useParticleWeights);  
290    bUsePhiWeights = (Int_t)useParticleWeights->GetBinContent(1);
291    bUsePtWeights = (Int_t)useParticleWeights->GetBinContent(2);
292    bUseEtaWeights = (Int_t)useParticleWeights->GetBinContent(3);
293   }
294   
295   // 3.) distributions and 4.) final results of fitting:
296   TString pWeightsFlag[2] = {"pWeights not used","pWeights used"};
297   TString sigmaFlag[2] = {"#sigma^{2} not fixed","#sigma^{2} fixed"};
298   
299   // q-distribution:
300   TString qDistributionName = "fqDistribution";
301   qDistributionName += fAnalysisLabel->Data();
302   // sum of particle weights:
303   TString sumOfParticleWeightsName = "fSumOfParticleWeightsName"; 
304   sumOfParticleWeightsName += fAnalysisLabel->Data();
305   // final results for integrated flow:
306   TString intFlowName = "fIntFlowFQD";
307   intFlowName += fAnalysisLabel->Data();
308   // sigma^2:
309   TString sigma2Name = "fSigma2";
310   sigma2Name += fAnalysisLabel->Data();
311   
312   TH1D *qDistribution[2] = {NULL};
313   TH1D *sumOfParticleWeights[2] = {NULL};
314   TH1D *intFlow[2][2] = {{NULL}};
315   TH1D *sigma2[2][2] = {{NULL}};
316    
317   for(Int_t pW=0;pW<1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights);pW++)
318   {
319    // q-distribution:
320    qDistribution[pW] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",qDistributionName.Data(),pWeightsFlag[pW].Data())));
321    if(qDistribution[pW])
322    {
323     this->SetqDistribution(qDistribution[pW],pW);
324    } else
325      {
326       cout<<"WARNING: qDistribution[pW] is NULL in AFAWFQD::GOH() !!!!"<<endl;
327       cout<<"pW = "<<pW<<endl;
328      }
329    // sum of particle weights:
330    sumOfParticleWeights[pW] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",sumOfParticleWeightsName.Data(),pWeightsFlag[pW].Data())));
331    if(sumOfParticleWeights[pW])
332    {
333     this->SetSumOfParticleWeights(sumOfParticleWeights[pW],pW);
334    } else
335      {
336       cout<<"WARNING: sumOfParticleWeights[pW] is NULL in AFAWFQD::GOH() !!!!"<<endl;
337       cout<<"pW = "<<pW<<endl;
338      }
339    // final results:
340    for(Int_t f=0;f<2;f++)
341    {
342     // final results for integrated flow:
343     intFlow[pW][f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s, %s",intFlowName.Data(),pWeightsFlag[pW].Data(),sigmaFlag[f].Data())));
344     if(intFlow[pW][f])
345     {
346      this->SetIntFlow(intFlow[pW][f],pW,f);
347     } else 
348       {
349        cout<<"WARNING: intFlow[pW][f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
350        cout<<"pW = "<<pW<<endl;
351        cout<<"f  = "<<f<<endl;
352       }
353     // sigma^2:
354     sigma2[pW][f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s, %s",sigma2Name.Data(),pWeightsFlag[pW].Data(),sigmaFlag[f].Data())));
355     if(sigma2[pW][f])
356     {
357      this->SetSigma2(sigma2[pW][f],pW,f);
358     } else 
359       {
360        cout<<"WARNING: sigma2[pW][f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
361        cout<<"pW = "<<pW<<endl;
362        cout<<"f  = "<<f<<endl;
363       } 
364    } // end of for(Int_t f=0;f<2;f++)
365   } // end of for(Int_t pW=0;pW<1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights);pW++)
366   
367   // 5.) fitting parameters:
368   // q-distribution:
369   TString fittingParametersName = "fFittingParameters";
370   fittingParametersName += fAnalysisLabel->Data();
371   TProfile *fittingParameters = NULL;
372   fittingParameters = dynamic_cast<TProfile*>(outputListHistos->FindObject(fittingParametersName.Data()));
373   if(fittingParameters)
374   {
375    this->SetFittingParameters(fittingParameters);
376   } else
377     {
378      cout<<"WARNING:fittingParameters is NULL in AFAWFQD::GOH() !!!!"<<endl;
379     }
380   
381  } else // to if(outputListHistos)
382    {
383     cout<<"WARNING: outputListHistos is NULL in AFAWFQD::GOH() !!!!"<<endl;
384     exit(0);
385    } 
386     
387    
388 } // end of void AliFlowAnalysisWithFittingQDistribution::GetOutputHistograms(TList *outputListHistos) 
389
390
391 //================================================================================================================
392
393
394 void AliFlowAnalysisWithFittingQDistribution::Finish(Bool_t doFit)
395 {
396  // calculate the final results
397  
398  // a) acces the constants and fitting paremeters;
399  // b) access the flags;
400  // c) do final fit;
401  // d) fill common hist results;
402  // e) print on the screen the final results.
403  
404  // access the constants and fitting paremeters:
405  this->AccessConstants();
406  this->AccessFittingParameters();
407  
408  // access the flags: 
409  fUsePhiWeights = (Int_t)fUseParticleWeights->GetBinContent(1); 
410  fUsePtWeights = (Int_t)fUseParticleWeights->GetBinContent(2); 
411  fUseEtaWeights = (Int_t)fUseParticleWeights->GetBinContent(3);
412
413  // to be improved (moved somewhere else):
414  if(fPlotResults)
415  {
416   fLegend = new TLegend(0.6,0.55,0.85,0.7); 
417  }
418
419  // do final fit:             
420  if(doFit) 
421  {
422   // particle weights not used:
423   // a) sigma^2 not fixed:
424   this->DoFit(kFALSE,kTRUE);
425   // b) sigma^2 fixed to 0.5:
426   this->DoFit(kFALSE,kFALSE);
427   // particle weights used:
428   if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)      
429   {
430    // a) sigma^2 not fixed:
431    this->DoFit(kTRUE,kTRUE);  
432    // b) sigma^2 fixed:
433    this->DoFit(kTRUE,kFALSE);  
434   }
435   
436   // fill common hist results (by default fill results obtained with sigma^2 not fixed):
437   if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
438   {
439    this->FillCommonHistResultsIntFlow(kTRUE,kFALSE); 
440   } else
441     {
442      this->FillCommonHistResultsIntFlow(kFALSE,kFALSE);    
443     } 
444   
445   // print final results on the screen:
446   this->PrintFinalResultsForIntegratedFlow();  
447  } // end of if(doFit)
448    
449 } // end of void AliFlowAnalysisWithFittingQDistribution::Finish(Bool_t doFit)
450
451
452 //================================================================================================================
453
454
455 void AliFlowAnalysisWithFittingQDistribution::WriteHistograms(TString* outputFileName)
456 {
457  //store the final results in output .root file
458  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
459  //output->WriteObject(fHistList, "cobjFQD","SingleKey");
460  fHistList->SetName("cobjFQD");
461  fHistList->SetOwner(kTRUE);
462  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
463  delete output;
464 }
465
466 //================================================================================================================
467
468 void AliFlowAnalysisWithFittingQDistribution::WriteHistograms(TString outputFileName)
469 {
470  //store the final results in output .root file
471  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
472  //output->WriteObject(fHistList, "cobjFQD","SingleKey");
473  fHistList->SetName("cobjFQD");
474  fHistList->SetOwner(kTRUE);
475  fHistList->Write(fHistList->GetName(), TObject::kSingleKey); 
476  delete output;
477 }
478
479
480 //================================================================================================================
481
482
483 void AliFlowAnalysisWithFittingQDistribution::InitializeArrays()
484 {
485  // initialize all arrays
486  
487  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++) // pW not used or used
488  {
489   fSumOfParticleWeights[pW] = NULL;
490   fqDistribution[pW] = NULL; 
491   for(Int_t f=0;f<2;f++) // sigma^2 not fixed or fixed
492   {
493    fIntFlow[pW][f] = NULL;
494    fSigma2[pW][f] = NULL;
495   }
496  } 
497
498 } // end of void AliFlowAnalysisWithFittingQDistribution::InitializeArrays()
499
500
501 //================================================================================================================
502
503
504 void AliFlowAnalysisWithFittingQDistribution::BookCommonHistograms()
505 {
506  // book common histograms
507  
508  // common control histogram: 
509  TString commonHistName = "AliFlowCommonHistFQD";
510  commonHistName += fAnalysisLabel->Data();
511  fCommonHists = new AliFlowCommonHist(commonHistName.Data());
512  fHistList->Add(fCommonHists);  
513
514  // common histograms for final results:
515  TString commonHistResName = "AliFlowCommonHistResultsFQD";
516  commonHistResName += fAnalysisLabel->Data();
517  fCommonHistsResults = new AliFlowCommonHistResults(commonHistResName.Data());
518  fHistList->Add(fCommonHistsResults); 
519
520 } // end of void AliFlowAnalysisWithFittingQDistribution::BookCommonHistograms(
521
522
523 //================================================================================================================
524
525
526 void AliFlowAnalysisWithFittingQDistribution::BookAndFillWeightsHistograms()
527 {
528  // book and fill histograms which hold phi, pt and eta weights
529
530  if(!fWeightsList)
531  {
532   cout<<"WARNING: fWeightsList is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
533   exit(0);  
534  }
535     
536  TString fUseParticleWeightsName = "fUseParticleWeightsFQD";
537  fUseParticleWeightsName += fAnalysisLabel->Data();
538  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
539  fUseParticleWeights->SetLabelSize(0.06);
540  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
541  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
542  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
543  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
544  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
545  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
546  fWeightsList->Add(fUseParticleWeights); 
547   
548  if(fUsePhiWeights)
549  {
550   if(fWeightsList->FindObject("phi_weights"))
551   {
552    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
553    if(fPhiWeights->GetBinWidth(1) != fPhiBinWidth)
554    {
555     cout<<"WARNING: fPhiWeights->GetBinWidth(1) != fPhiBinWidth in AFAWFQD::BAFWH() !!!!        "<<endl;
556     cout<<"         This indicates inconsistent binning in phi histograms throughout the code."<<endl;
557     exit(0);
558    }
559   } else 
560     {
561      cout<<"WARNING: fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
562      exit(0);
563     }
564  } // end of if(fUsePhiWeights)
565  
566  if(fUsePtWeights) 
567  {
568   if(fWeightsList->FindObject("pt_weights"))
569   {
570    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
571    if(fPtWeights->GetBinWidth(1) != fPtBinWidth)
572    {
573     cout<<"WARNING: fPtWeights->GetBinWidth(1) != fPtBinWidth in AFAWFQD::BAFWH() !!!!         "<<endl;
574     cout<<"         This indicates insconsistent binning in pt histograms throughout the code."<<endl;
575     exit(0);
576    }
577   } else 
578     {
579      cout<<"WARNING: fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
580      exit(0);
581     }
582  } // end of if(fUsePtWeights)    
583
584  if(fUseEtaWeights) 
585  {
586   if(fWeightsList->FindObject("eta_weights"))
587   {
588    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
589    if(fEtaWeights->GetBinWidth(1) != fEtaBinWidth)
590    {
591     cout<<"WARNING: fEtaWeights->GetBinWidth(1) != fEtaBinWidth in AFAWFQD::BAFWH() !!!!        "<<endl;
592     cout<<"         This indicates insconsistent binning in eta histograms throughout the code."<<endl;
593     exit(0);
594    }
595   } else 
596     {
597      cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
598      exit(0);
599     }
600  } // end of if(fUseEtaWeights)
601  
602 } // end of AliFlowAnalysisWithFittingQDistribution::BookAndFillWeightsHistograms()
603
604
605 //================================================================================================================================
606
607
608 void AliFlowAnalysisWithFittingQDistribution::AccessConstants()
609 {
610  // access needed common constants from AliFlowCommonConstants
611  
612  fnBinsPhi = AliFlowCommonConstants::GetNbinsPhi();
613  fPhiMin = AliFlowCommonConstants::GetPhiMin();      
614  fPhiMax = AliFlowCommonConstants::GetPhiMax();
615  if(fnBinsPhi) fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;  
616  fnBinsPt = AliFlowCommonConstants::GetNbinsPt();
617  fPtMin = AliFlowCommonConstants::GetPtMin();        
618  fPtMax = AliFlowCommonConstants::GetPtMax();
619  if(fnBinsPt) fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;  
620  fnBinsEta = AliFlowCommonConstants::GetNbinsEta();
621  fEtaMin = AliFlowCommonConstants::GetEtaMin();      
622  fEtaMax = AliFlowCommonConstants::GetEtaMax();
623  if(fnBinsEta) fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;  
624  
625 } // end of void AliFlowAnalysisWithFittingQDistribution::AccessConstants()
626
627
628 //================================================================================================================================
629
630
631 void AliFlowAnalysisWithFittingQDistribution::BookEverythingForDistributions()
632 {
633  // book histograms for distributions
634  
635  TString pWeightsFlag[2] = {"pWeights not used","pWeights used"};
636  TString sigmaFlag[2] = {"#sigma^{2} not fixed","#sigma^{2} fixed"};
637  // q-distribution:
638  TString fqDistributionName = "fqDistribution";
639  fqDistributionName += fAnalysisLabel->Data();
640  // sum of particle weights:
641  TString fSumOfParticleWeightsName = "fSumOfParticleWeightsName";
642  fSumOfParticleWeightsName += fAnalysisLabel->Data();
643  // final results for integrated flow:
644  TString fIntFlowName = "fIntFlowFQD";
645  fIntFlowName += fAnalysisLabel->Data();
646  // sigma^2:
647  TString fSigma2Name = "fSigma2";
648  fSigma2Name += fAnalysisLabel->Data();
649  
650  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++) // pW not used or used
651  {
652   // q-distribution:
653   fqDistribution[pW] = new TH1D(Form("%s, %s",fqDistributionName.Data(),pWeightsFlag[pW].Data()),"q-distribution",10000,0,1000);
654   fqDistribution[pW]->SetXTitle("q_{n}=Q_{n}/#sqrt{M}");
655   fqDistribution[pW]->SetYTitle("Counts");
656   fHistList->Add(fqDistribution[pW]);
657   // sum of particle weights: 
658   fSumOfParticleWeights[pW] = new TH1D(Form("%s, %s",fSumOfParticleWeightsName.Data(),pWeightsFlag[pW].Data()),"Sum of particle weights",1000,0,10000);
659   fSumOfParticleWeights[pW]->SetXTitle("#sum_{i=1}^{N} w_{i}");
660   fSumOfParticleWeights[pW]->SetYTitle("Counts");
661   fHistList->Add(fSumOfParticleWeights[pW]);
662   
663   for(Int_t f=0;f<2;f++) // sigma^2 not fixed or fixed
664   {
665    // final results for integrated flow:
666    fIntFlow[pW][f] = new TH1D(Form("%s, %s, %s",fIntFlowName.Data(),pWeightsFlag[pW].Data(),sigmaFlag[f].Data()),"Integrated Flow",1,0,1);
667    fIntFlow[pW][f]->SetLabelSize(0.08);
668    (fIntFlow[pW][f]->GetXaxis())->SetBinLabel(1,"v_{n}");
669    fHistList->Add(fIntFlow[pW][f]);
670    // sigma^2:
671    fSigma2[pW][f] = new TH1D(Form("%s, %s, %s",fSigma2Name.Data(),pWeightsFlag[pW].Data(),sigmaFlag[f].Data()),"#sigma^{2}",1,0,1);
672    fSigma2[pW][f]->SetLabelSize(0.08);
673    (fSigma2[pW][f]->GetXaxis())->SetBinLabel(1,"#sigma^{2}");
674    fHistList->Add(fSigma2[pW][f]);
675   } // end of for(Int_t f=0;f<2;f++) // sigma^2 not fixed or fixed
676   
677  } // end of for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++) // pW not used or used
678  
679  // book profile fFittingParameters whiuch will hols fitting parameters:
680  TString fFittingParametersName = "fFittingParameters";
681  fFittingParametersName += fAnalysisLabel->Data(); 
682  fFittingParameters = new TProfile(fFittingParametersName.Data(),"Parameters for fitting q-distribution",8,0,8);
683  fFittingParameters->SetLabelSize(0.05);
684  (fFittingParameters->GetXaxis())->SetBinLabel(1,"Treshold");
685  (fFittingParameters->GetXaxis())->SetBinLabel(2,"starting v_{n}");
686  (fFittingParameters->GetXaxis())->SetBinLabel(3,"min. v_{n}");
687  (fFittingParameters->GetXaxis())->SetBinLabel(4,"max. v_{n}");
688  (fFittingParameters->GetXaxis())->SetBinLabel(5,"starting #sigma^{2}");
689  (fFittingParameters->GetXaxis())->SetBinLabel(6,"min. #sigma^{2}");
690  (fFittingParameters->GetXaxis())->SetBinLabel(7,"max. #sigma^{2}");
691  (fFittingParameters->GetXaxis())->SetBinLabel(8,"plot or not?");
692  fHistList->Add(fFittingParameters);
693  
694 } // end of void AliFlowAnalysisWithFittingQDistribution::BookEverythingForDistributions()
695
696
697 //================================================================================================================================
698
699
700 void AliFlowAnalysisWithFittingQDistribution::DoFit(Bool_t useParticleWeights, Bool_t sigma2NotFixed)
701 {
702  // do final fit to q-distribution
703  
704  // shortcuts for flags:
705  Int_t pW = (Int_t)(useParticleWeights);
706  Int_t s2NF = 1-(Int_t)(sigma2NotFixed);
707  
708  for(Int_t f=0;f<2;f++)
709  {
710   if(!(fqDistribution[pW] && fSumOfParticleWeights[pW] && fIntFlow[pW][f] && fSigma2[pW][f])) 
711   { 
712    cout<<"WARNING: fqDistribution[pW] && fSumOfParticleWeights[pW] && fIntFlow[pW][f] && fSigma2[pw][f] is NULL in AFAWFQD::DoFit() !!!!"<<endl;
713    cout<<"pW = "<<pW<<endl;
714    cout<<"f  = "<<f<<endl;
715    exit(0);
716   }
717  }
718  
719  // average multiplicity and number of events:
720  Double_t AvM = fSumOfParticleWeights[pW]->GetMean(1);
721  //Int_t nEvts = (Int_t)fSumOfParticleWeights[pW]->GetEntries();
722  
723  // for fitting take into account only bins with at least fTreshold entries:
724  Int_t binMin = fqDistribution[pW]->FindFirstBinAbove(fTreshold); // to be improved (add setter for this)  
725  Int_t binMax = fqDistribution[pW]->FindLastBinAbove(fTreshold); // to be improved (add setter for this) 
726  Double_t binWidth = fqDistribution[pW]->GetBinWidth(4); // assuming that all bins have the same width 
727  if(binWidth == 0) 
728  {
729   cout<<"WARNING: binWidth == 0 in AFAWFQD::DoFit()"<<endl;
730   exit(0);
731  }
732  Double_t qmin = (binMin-1)*binWidth; 
733  Double_t qmax = (binMax)*binWidth;
734  Double_t ent = 0.; // number of entries between binMin and binMax:
735  for(Int_t b=binMin;b<=binMax;b++)
736  {
737   ent += fqDistribution[pW]->GetBinContent(b);
738  }
739  Double_t norm = binWidth*ent; // norm (assuming that all bins have the same width)
740
741  // fitting function:
742  TF1 *fittingFun = new TF1("fittingFun","[2]*(x/[1])*exp(-(x*x+[0]*[0])/(2.*[1]))*TMath::BesselI0(x*[0]/[1])",qmin,qmax); 
743  
744  fittingFun->SetParNames("v*sqrt{sum of particle weights}","sigma^2","norm");
745  fittingFun->SetParameters(fvStart*pow(AvM,0.5),fSigma2Start,norm); // to be improved (add setter for starting v)         
746  fittingFun->SetParLimits(0,fvMin*pow(AvM,0.5),fvMax*pow(AvM,0.5)); // to be improved (add setters for vmin and vmax)
747  
748  if(s2NF == 0)
749  {
750   fittingFun->SetParLimits(1,fSigma2Min,fSigma2Max); // to be improved (add setters for sigma^2_min and sigma^2_max)      
751  } else
752    {
753     fittingFun->FixParameter(1,0.5);
754    }
755  fittingFun->FixParameter(2,norm);  
756
757  // fitting (do it only if # of entries >50): // to be improved (this is only a pragmatics fix to avoid TMinuit crash)
758  if(ent > 50)
759  {
760   fqDistribution[pW]->Fit("fittingFun","NQ","",qmin,qmax);
761  }
762  // results:
763  Double_t v = 0.; // integrated flow
764  Double_t vError = 0.; // error of integrated flow 
765  Double_t sigma2 = 0.; // sigma^2
766  Double_t sigma2Error = 0.; // error of sigma^2
767  
768  if(AvM)
769  { 
770   v = fittingFun->GetParameter(0)/pow(AvM,0.5);
771   vError = fittingFun->GetParError(0)/pow(AvM,0.5);
772   // store the results:
773   fIntFlow[pW][s2NF]->SetBinContent(1,v);
774   fIntFlow[pW][s2NF]->SetBinError(1,vError);
775  }
776  if(s2NF == 0)
777  {
778   sigma2 = fittingFun->GetParameter(1);
779   sigma2Error = fittingFun->GetParError(1);
780   // store the results:
781   fSigma2[pW][0]->SetBinContent(1,sigma2);  
782   fSigma2[pW][0]->SetBinError(1,sigma2Error);  
783  } else // sigma^2 not fitted, but fixed to 0.5
784    {
785     sigma2 = 0.5;
786     fSigma2[pW][1]->SetBinContent(1,sigma2);  
787     fSigma2[pW][1]->SetBinError(1,0.);
788    }
789  
790  if(fPlotResults)
791  {
792   if(s2NF == 0)
793   { 
794    fqDistribution[pW]->SetFillColor(16);  
795    fqDistribution[pW]->SetTitle("Fitted q-distribution");
796    fqDistribution[pW]->Draw("");
797    fLegend->AddEntry(fqDistribution[pW],"q-distribution","f"); 
798    fittingFun->SetLineColor(2); // 2 = red color   
799    fittingFun->Draw("SAME");
800    fLegend->AddEntry("fittingFun","#sigma^{2} fitted","l");
801   } else
802     {
803      // to be improved (perhaps there is a better way to implement this?)
804      TF1 *fittingFunTemp = (TF1*)fittingFun->Clone("fittingFunTemp");
805      fittingFunTemp->SetLineColor(4); // 4 = blue color
806      fittingFunTemp->Draw("SAME"); 
807      fLegend->AddEntry("fittingFunTemp","#sigma^{2} fixed","l");
808      fLegend->Draw("SAME");      
809     } 
810  } // end of if(fPlotResults)
811
812 } // end of void AliFlowAnalysisWithFittingQDistribution::DoFit(Bool_t useParticleWeights)
813
814
815 //================================================================================================================================ 
816
817
818 void AliFlowAnalysisWithFittingQDistribution::FillCommonHistResultsIntFlow(Bool_t useParticleWeights, Bool_t sigma2NotFixed)
819 {
820  // fill in AliFlowCommonHistResults histograms relevant for 'NONAME' integrated flow (to be improved (name))
821  
822  // shortcuts for the flags:
823  Int_t pW = (Int_t)(useParticleWeights); // 0 = pWeights not useed, 1 = pWeights used
824  Int_t s2NF = (Int_t)(sigma2NotFixed); // 0 = sigma^2 not fixed, 1 = sigma^2 fixed to 0.5
825  
826  if(!fIntFlow[pW][s2NF])
827  {
828   cout<<"WARNING: fIntFlow[pW][s2NF] is NULL in AFAWFQD::FCHRIF() !!!!"<<endl;
829   cout<<"pW   = "<<pW<<endl;
830   cout<<"s2NF = "<<s2NF<<endl;
831   exit(0); 
832  }  
833  
834  if(!fSumOfParticleWeights[pW])
835  {
836   cout<<"WARNING: fSumOfParticleWeights[pW] is NULL in AFAWFQD::FCHRIF() !!!!"<<endl;
837   cout<<"pW = "<<pW<<endl;
838   exit(0);
839  }
840  
841  if(!(fCommonHistsResults))
842  {
843   cout<<"WARNING: fCommonHistsResults is NULL in AFAWFQD::FCHRIF() !!!!"<<endl; 
844   exit(0);
845  }
846   
847  // fill integrated flow:
848  Double_t v = fIntFlow[pW][s2NF]->GetBinContent(1); 
849  Double_t vError = fIntFlow[pW][s2NF]->GetBinError(1);
850  
851  fCommonHistsResults->FillIntegratedFlow(v,vError);   
852  
853  // fill chi:
854  Double_t AvM = fSumOfParticleWeights[pW]->GetMean(1);
855  Double_t chi = AvM*pow(v,2); 
856  if(chi>=0)
857  {
858   fCommonHistsResults->FillChi(pow(chi,0.5));   
859   fCommonHistsResults->FillChiRP(pow(chi,0.5));   
860  }
861    
862 } // end of void AliFlowAnalysisWithFittingQDistribution::FillCommonHistResultsIntFlow(Bool_t useParticleWeights, Bool_t sigma2NotFixed) 
863
864
865 //================================================================================================================================ 
866
867
868 void AliFlowAnalysisWithFittingQDistribution::PrintFinalResultsForIntegratedFlow()
869 {
870  // print the final results for integrated flow on the screen
871  
872  // shortcuts: pW   = particle weights 
873  //            s2NF = sigma^2 not fixed  
874  
875  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++)
876  {
877   if(!fSumOfParticleWeights[pW])
878   {
879    cout<<"WARNING: fSumOfParticleWeights[pW] is NULL in AFAWFQD::FCHRIF() !!!!"<<endl;
880    cout<<"pW = "<<pW<<endl;
881    exit(0);
882   }
883   for(Int_t s2NF=0;s2NF<2;s2NF++)
884   {
885    if(!fIntFlow[pW][s2NF])
886    {
887     cout<<"WARNING: fIntFlow[pW][s2NF] is NULL in AFAWFQD::FCHRIF() !!!!"<<endl;
888     cout<<"pW   = "<<pW<<endl;
889     cout<<"s2NF = "<<s2NF<<endl;
890     exit(0); 
891    }
892   }  
893  }  
894  
895  if(!(fCommonHistsResults))
896  {
897   cout<<"WARNING: fCommonHistsResults is NULL in AFAWFQD::FCHRIF() !!!!"<<endl; 
898   exit(0);
899  }
900  
901  // shortcut for the harmonic:
902  Int_t n = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1); 
903
904  // printing:
905  cout<<" "<<endl;
906  cout<<"***************************************"<<endl;
907  cout<<"***************************************"<<endl;
908  cout<<"      integrated flow by fitting "<<endl;
909  cout<<"           q-distribution:      "<<endl;
910  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
911  {
912   cout<<"           (with weights)       "<<endl;
913  } else
914    {
915     cout<<"          (without weights)       "<<endl;
916    }   
917  cout<<endl;
918
919  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
920  {
921   cout<<"  v_"<<n<<"{FQD} = "<<fIntFlow[1][0]->GetBinContent(1)<<" +/- "<<fIntFlow[1][0]->GetBinError(1)
922     <<" (sigma^2 fitted = "<<fSigma2[1][0]->GetBinContent(1)<<")"<<endl; 
923   cout<<"  v_"<<n<<"{FQD} = "<<fIntFlow[1][1]->GetBinContent(1)<<" +/- "<<fIntFlow[1][1]->GetBinError(1)
924     <<" (sigma^2 fixed = 0.5"<<")"<<endl; 
925   cout<<" "<<endl; 
926   cout<<"      nEvts = "<<fSumOfParticleWeights[1]->GetEntries()<<", AvM = "<<fSumOfParticleWeights[1]->GetMean()<<endl;
927   cout<<" "<<endl;
928  } else
929    { 
930     cout<<"  v_"<<n<<"{FQD} = "<<fIntFlow[0][0]->GetBinContent(1)<<" +/- "<<fIntFlow[0][0]->GetBinError(1)
931      <<" (sigma^2 fitted = "<<fSigma2[0][0]->GetBinContent(1)<<")"<<endl; 
932     cout<<"  v_"<<n<<"{FQD} = "<<fIntFlow[0][1]->GetBinContent(1)<<" +/- "<<fIntFlow[0][1]->GetBinError(1)
933      <<" (sigma^2 fixed = 0.5"<<")"<<endl; 
934     cout<<" "<<endl;  
935     cout<<"      nEvts = "<<fSumOfParticleWeights[0]->GetEntries()<<", AvM = "<<fSumOfParticleWeights[0]->GetMean()<<endl;
936     cout<<" "<<endl;
937    }
938     
939  cout<<"***************************************"<<endl;
940  cout<<"***************************************"<<endl; 
941  
942 } // end of void AliFlowAnalysisWithFittingQDistribution::PrintFinalResultsForIntegratedFlow()
943
944
945 //================================================================================================================================ 
946
947
948 void AliFlowAnalysisWithFittingQDistribution::StoreFittingParameters()
949 {
950  // store fitting parameters in profile fFittingParameters
951  
952  // Binning of fFittingParameters is organized as follows:
953  // 1st bin: fTreshold
954  // 2nd bin: fvStart
955  // 3rd bin: fvMin
956  // 4th bin: fvMax
957  // 5th bin: fSigma2Start
958  // 6th bin: fSigma2Min
959  // 7th bin: fSigma2Max
960  
961  if(!fFittingParameters)
962  {
963   cout<<"WARNING: fFittingParameters is NULL in AFAWFQD::SFP() !!!!"<<endl;
964   exit(0);
965  }
966  
967  fFittingParameters->Reset();
968  fFittingParameters->Fill(0.5,fTreshold);
969  fFittingParameters->Fill(1.5,fvStart);
970  fFittingParameters->Fill(2.5,fvMin);
971  fFittingParameters->Fill(3.5,fvMax);
972  fFittingParameters->Fill(4.5,fSigma2Start);
973  fFittingParameters->Fill(5.5,fSigma2Min);
974  fFittingParameters->Fill(6.5,fSigma2Max);
975  fFittingParameters->Fill(7.5,fPlotResults);
976  
977 } // end of void AliFlowAnalysisWithFittingQDistribution::StoreFittingParameters()
978
979
980 //================================================================================================================================ 
981
982
983 void AliFlowAnalysisWithFittingQDistribution::AccessFittingParameters()
984 {
985  // access fitting parameters:
986  
987  if(!fFittingParameters)
988  {
989   cout<<"WARNING: fFittingParameters is NULL in AFAWFQD::AFP() !!!!"<<endl;
990   exit(0);
991  }
992  
993  fTreshold = fFittingParameters->GetBinContent(1);
994  fvStart = fFittingParameters->GetBinContent(2);
995  fvMin = fFittingParameters->GetBinContent(3);
996  fvMax = fFittingParameters->GetBinContent(4);
997  fSigma2Start = fFittingParameters->GetBinContent(5);
998  fSigma2Min = fFittingParameters->GetBinContent(6);
999  fSigma2Max = fFittingParameters->GetBinContent(7);
1000  fPlotResults = (Bool_t) fFittingParameters->GetBinContent(8);
1001  
1002 } // end of void AliFlowAnalysisWithFittingQDistribution::AccessFittingParameters()