AliPhysics  f9b5d69 (f9b5d69)
AliFlowAnalysisWithFittingQDistribution.cxx
Go to the documentation of this file.
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  * estimating reference flow by *
18  * fitting q-distribution *
19  * *
20  * author: Ante Bilandzic *
21  * (abilandzic@gmail.com) *
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"
33 #include "TChain.h"
34 #include "TFile.h"
35 #include "TList.h"
36 #include "TF1.h"
37 #include "TH2D.h"
38 #include "TParticle.h"
39 #include "TProfile.h"
40 #include "AliFlowEventSimple.h"
41 #include "AliFlowTrackSimple.h"
43 
44 class TH1;
45 class TGraph;
46 class TPave;
47 class TLatex;
48 class TMarker;
49 class TObjArray;
50 class TList;
51 class TSystem;
52 class TROOT;
53 class AliFlowVector;
54 class TVector;
55 
56 //================================================================================================================
57 
58 using std::endl;
59 using std::cout;
61 
63  fHistList(NULL),
64  fBookOnlyBasicCCH(kTRUE),
65  fCommonHists(NULL),
66  fCommonHistsResults(NULL),
67  fnBinsPhi(0),
68  fPhiMin(0),
69  fPhiMax(0),
70  fPhiBinWidth(0),
71  fnBinsPt(0),
72  fPtMin(0),
73  fPtMax(0),
74  fPtBinWidth(0),
75  fnBinsEta(0),
76  fEtaMin(0),
77  fEtaMax(0),
78  fEtaBinWidth(0),
79  fHarmonic(2),
80  fAnalysisLabel(NULL),
81  fMultiplicityIs(AliFlowCommonConstants::kRP),
82  fWeightsList(NULL),
83  fUsePhiWeights(kFALSE),
84  fUsePtWeights(kFALSE),
85  fUseEtaWeights(kFALSE),
86  fUseParticleWeights(NULL),
87  fPhiWeights(NULL),
88  fPtWeights(NULL),
89  fEtaWeights(NULL),
90  fSumOfParticleWeights(NULL),
91  fqDistribution(NULL),
92  fqMin(0.),
93  fqMax(100.),
94  fqNbins(10000),
95  fStoreqDistributionVsMult(kFALSE),
96  fqDistributionVsMult(NULL),
97  fMinMult(0.),
98  fMaxMult(10000.),
99  fnBinsMult(1000),
100  fFittingParameters(NULL),
101  fTreshold(5.),
102  fvStart(0.05),
103  fvMin(0.0),
104  fvMax(0.25),
105  fSigma2Start(0.75),
106  fSigma2Min(0.5),
107  fSigma2Max(2.5),
108  fFinalResultIsFromSigma2Fitted(kTRUE),
109  fPrintOnTheScreen(kTRUE),
110  fDoFit(kTRUE),
111  fExactNoRPs(0)
112  {
113  // constructor
114 
115  // base list to hold all output objects:
116  fHistList = new TList();
117  fHistList->SetName("cobjFQD");
118  fHistList->SetOwner(kTRUE);
119 
120  // analysis label;
121  fAnalysisLabel = new TString();
122 
123  // list to hold histograms with phi, pt and eta weights:
124  fWeightsList = new TList();
125 
126  // initialize all arrays:
127  this->InitializeArrays();
128 
129  } // end of constructor
130 
131 //================================================================================================================
132 
134 {
135  // destructor
136 
137  delete fHistList;
138 
139 } // end of destructor
140 
141 //================================================================================================================
142 
143 
145 {
146  // Access constants and book everything.
147 
148  //save old value and prevent histograms from being added to directory
149  //to avoid name clashes in case multiple analaysis objects are used
150  //in an analysis
151  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
152  TH1::AddDirectory(kFALSE);
153 
154  // access constants:
155  this->AccessConstants();
156 
157  // booking:
158  this->BookCommonHistograms();
161 
162  // store fitting parameters:
163  this->StoreFittingParameters();
164 
165  // nest lists:
166  fWeightsList->SetName("Weights");
167  fWeightsList->SetOwner(kTRUE);
168  fHistList->Add(fWeightsList);
169 
170  // set harmonic in common control histograms (to be improved (should I do this somewhere else?)):
171  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
172 
173  //restore old status
174  TH1::AddDirectory(oldHistAddStatus);
175 } // end of void AliFlowAnalysisWithFittingQDistribution::Init()
176 
177 //================================================================================================================
178 
180 {
181  // Loop over data only in this method.
182 
183  // a) Check all pointers used in this method;
184  // b) Fill the common control histograms;
185  // c) Loop over data and calculate Q-vector and sum of particle weights;
186  // d) Fill the histogram for q-distribution and sum of particle weights.
187 
188  // a) Check all pointers used in this method:
189  this->CheckPointersUsedInMake();
190 
191  // b) fill the common control histograms:
193 
194  // c) Loop over data and fill histogram for q-distribution:
195  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
196  Double_t dPt = 0.; // transverse momentum
197  Double_t dEta = 0.; // pseudorapidity
198  Double_t wPhi = 1.; // phi weight
199  Double_t wPt = 1.; // pt weight
200  Double_t wEta = 1.; // eta weight
201  Double_t dReQ = 0.; // real part of Q-vector
202  Double_t dImQ = 0.; // imaginary part of Q-vector
203  Int_t nCounterNoRPs = 0; // needed only for shuffling
204  Int_t n = fHarmonic; // shortcut for the harmonic
205  Double_t dSumOfParticleWeights = 0.; // when particle weights are not used dSumOfParticleWeights is equal to multiplicity
206  AliFlowTrackSimple *aftsTrack = NULL;
207  Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks
208  if(fExactNoRPs > 0 && anEvent->GetNumberOfRPs()<fExactNoRPs){return;} // shuffling
209 
210  // Start loop over particles:
211  for(Int_t i=0;i<nPrim;i++)
212  {
213  if(fExactNoRPs > 0 && nCounterNoRPs>fExactNoRPs){continue;} // needed only for shuffling
214  aftsTrack=anEvent->GetTrack(i);
215  if(aftsTrack)
216  {
217  if(!(aftsTrack->InRPSelection())){continue;} // consider only tracks which are RPs
218  nCounterNoRPs++;
219  dPhi = aftsTrack->Phi();
220  dPt = aftsTrack->Pt();
221  dEta = aftsTrack->Eta();
222  if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi weight for this particle:
223  {
224  wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
225  }
226  if(fUsePtWeights && fPtWeights && fPtBinWidth) // determine pt weight for this particle:
227  {
228  wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
229  }
230  if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta weight for this particle:
231  {
232  wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
233  }
234  // Calculate real and imaginary part of Q-vector and sum of particle weights for this event:
235  // Q-vector:
236  dReQ += wPhi*wPt*wEta*TMath::Cos(n*dPhi);
237  dImQ += wPhi*wPt*wEta*TMath::Sin(n*dPhi);
238  // sum of particle weights:
239  dSumOfParticleWeights += wPhi*wPt*wEta; // when particle weights are not used this sum gives # of RPs, i.e. multiplicity
240  } // end of if(aftsTrack)
241  } // end of for(Int_t i=0;i<nPrim;i++)
242 
243  // d) Fill the histogram for q-distribution and sum of particle weights:
244  Double_t q = 0.; // q = Q\sqrt{sum of particle weights}
245  if(dSumOfParticleWeights > 1.)
246  {
247  q = pow(dReQ*dReQ+dImQ*dImQ,0.5)/pow(dSumOfParticleWeights,0.5);
248  fqDistribution->Fill(q,1.);
249  fSumOfParticleWeights->Fill(dSumOfParticleWeights,1.);
251  {
252  // Multiplicity bin of an event (relevant for all histos vs M):
253  Double_t dMultiplicityBin = 0.;
255  {
256  Int_t nNumberOfRPsEBE = anEvent->GetNumberOfRPs(); // number of RPs (i.e. number of reference particles)
257  dMultiplicityBin = nNumberOfRPsEBE+0.5;
259  {
260  Int_t nReferenceMultiplicityEBE = anEvent->GetReferenceMultiplicity(); // reference multiplicity for current event
261  dMultiplicityBin = nReferenceMultiplicityEBE+0.5;
263  {
264  Int_t nNumberOfPOIsEBE = anEvent->GetNumberOfPOIs(); // number of POIs (i.e. number of particles of interest)
265  dMultiplicityBin = nNumberOfPOIsEBE+0.5;
266  }
267  // Fill qDist vs M:
268  fqDistributionVsMult->Fill(dMultiplicityBin,q);
269  } // end of if(fStoreqDistributionVsMult)
270  } // end of if(dSumOfParticleWeights > 1.)
271 
272 } // end of Make()
273 
274 //================================================================================================================
275 
277 {
278  // Calculate the final results.
279 
280  // a) Check all pointers used in this method;
281  // b) Acces common constants;
282  // c) Access fitting paremeters;
283  // d) Do the final fit of q-distribution;
284  // e) Fill common hist results;
285  // f) Print on the screen the final results.
286 
287  // a) Check all pointers used in this method:
288  this->CheckPointersUsedInFinish();
289 
290  // b) Access common constants:
291  this->AccessConstants();
292 
293  // c) Access fitting paremeters:
294  this->AccessFittingParameters();
295 
296  // d) Do the final fit of q-distribution:
297  if(doFit && fDoFit)
298  {
299  this->DoFit(kFALSE); // sigma^2 not fitted (fixed to 0.5)
300  this->DoFit(kTRUE); // sigma^2 fitted
301  }
302 
303  // e) Fill common hist results (by default fill with results obtained with sigma^2 fitted,
304  // alternatively use a setter SetFinalResultIsFromSigma2Fitted(kFALSE)):
306 
307  // f) Print on the screen the final results:
308  if(fPrintOnTheScreen) this->PrintOnTheScreen();
309 
310 } // end of void AliFlowAnalysisWithFittingQDistribution::Finish(Bool_t doFit)
311 
312 //================================================================================================================
313 
315 {
316  // Check all pointers used in method Make().
317 
318  if(!fCommonHists)
319  {
320  cout<<endl;
321  cout<<" WARNING (FQD::Make()): fCommonHists is NULL !!!!"<<endl;
322  cout<<endl;
323  exit(0);
324  }
325  if(!fqDistribution)
326  {
327  cout<<endl;
328  cout<<" WARNING (FQD::Make()): fqDistribution is NULL !!!!"<<endl;
329  cout<<endl;
330  exit(0);
331  }
333  {
334  cout<<endl;
335  cout<<" WARNING (FQD::Make()): fSumOfParticleWeights is NULL !!!!"<<endl;
336  cout<<endl;
337  exit(0);
338  }
340  {
341  cout<<endl;
342  cout<<" WARNING (FQD::Make()): fPhiWeights is NULL !!!!"<<endl;
343  cout<<endl;
344  exit(0);
345  }
346  if(fUsePtWeights && !fPtWeights)
347  {
348  cout<<endl;
349  cout<<" WARNING (FQD::Make()): fPtWeights is NULL !!!!"<<endl;
350  cout<<endl;
351  exit(0);
352  }
354  {
355  cout<<endl;
356  cout<<" WARNING (FQD::Make()): fEtaWeights is NULL !!!!"<<endl;
357  cout<<endl;
358  exit(0);
359  }
360 
361 } // end of void AliFlowAnalysisWithFittingQDistribution::CheckPointersUsedInMake()
362 
363 //================================================================================================================
364 
366 {
367  // Check all pointers used in method Finish().
368 
369  if(!fFittingParameters)
370  {
371  cout<<endl;
372  cout<<" WARNING (FQD::Finish()): fFittingParameters is NULL !!!!"<<endl;
373  cout<<endl;
374  exit(0);
375  }
376  if(!fqDistribution)
377  {
378  cout<<endl;
379  cout<<" WARNING (FQD::Finish()): fqDistribution is NULL !!!!"<<endl;
380  cout<<endl;
381  exit(0);
382  }
384  {
385  cout<<endl;
386  cout<<" WARNING (FQD::Finish()): fSumOfParticleWeights is NULL !!!!"<<endl;
387  cout<<endl;
388  exit(0);
389  }
390  for(Int_t s2F=0;s2F<2;s2F++)
391  {
392  if(!fIntFlow[s2F])
393  {
394  cout<<endl;
395  cout<<" WARNING (FQD::Finish()): "<<Form("fIntFlow[%d] is NULL !!!!",s2F)<<endl;
396  cout<<endl;
397  exit(0);
398  }
399  if(!fSigma2[s2F])
400  {
401  cout<<endl;
402  cout<<" WARNING (FQD::Finish()): "<<Form("fSigma2[%d] is NULL !!!!",s2F)<<endl;
403  cout<<endl;
404  exit(0);
405  }
406  if(!fChi2[s2F])
407  {
408  cout<<endl;
409  cout<<" WARNING (FQD::Finish()): "<<Form("fChi2[%d] is NULL !!!!",s2F)<<endl;
410  cout<<endl;
411  exit(0);
412  }
413  } // end of for(Int_t s2F=0;s2F<2;s2F++)
415  {
416  cout<<endl;
417  cout<<"WARNING (FQD::Finish()): fCommonHistsResults is NULL !!!!"<<endl;
418  cout<<endl;
419  exit(0);
420  }
421  if(!fCommonHists)
422  {
423  cout<<endl;
424  cout<<"WARNING (FQD::Finish()): fCommonHists is NULL !!!!"<<endl;
425  cout<<endl;
426  exit(0);
427  }
428 
429 } // end of void AliFlowAnalysisWithFittingQDistribution::CheckPointersUsedInFinish()
430 
431 //================================================================================================================
432 
434 {
435  // Get pointers to all output histograms (called before Finish()).
436 
437  if(outputListHistos)
438  {
439  // 1.) common control histograms and common histograms for final results:
440  TString commonHistName = "AliFlowCommonHistFQD";
441  commonHistName += fAnalysisLabel->Data();
442  AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject(commonHistName.Data()));
443  if(commonHist) this->SetCommonHists(commonHist);
444 
445  TString commonHistResName = "AliFlowCommonHistResultsFQD";
446  commonHistResName += fAnalysisLabel->Data();
447  AliFlowCommonHistResults *commonHistRes = dynamic_cast<AliFlowCommonHistResults*>
448  (outputListHistos->FindObject(commonHistResName.Data()));
449  if(commonHistRes) this->SetCommonHistsResults(commonHistRes);
450 
451  // 2.) weights:
452  TList *weightsList = dynamic_cast<TList*>(outputListHistos->FindObject("Weights"));
453  if(weightsList){this->SetWeightsList(weightsList);}
454  Bool_t bUsePhiWeights = kFALSE;
455  Bool_t bUsePtWeights = kFALSE;
456  Bool_t bUseEtaWeights = kFALSE;
457  TString fUseParticleWeightsName = "fUseParticleWeightsFQD";
458  fUseParticleWeightsName += fAnalysisLabel->Data();
459  TProfile *useParticleWeights = NULL;
460  if(weightsList)
461  {
462  useParticleWeights = dynamic_cast<TProfile*>(weightsList->FindObject(fUseParticleWeightsName.Data()));
463  }
464  if(useParticleWeights)
465  {
466  this->SetUseParticleWeights(useParticleWeights);
467  bUsePhiWeights = (Int_t)useParticleWeights->GetBinContent(1);
468  this->SetUsePhiWeights(bUsePhiWeights);
469  bUsePtWeights = (Int_t)useParticleWeights->GetBinContent(2);
470  this->SetUsePtWeights(bUsePtWeights);
471  bUseEtaWeights = (Int_t)useParticleWeights->GetBinContent(3);
472  this->SetUseEtaWeights(bUseEtaWeights);
473  }
474 
475  // 3.) distributions and 4.) final results of fitting:
476  TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
477  // q-distribution:
478  TString qDistributionName = "fqDistribution";
479  qDistributionName += fAnalysisLabel->Data();
480  // sum of particle weights:
481  TString sumOfParticleWeightsName = "fSumOfParticleWeights";
482  sumOfParticleWeightsName += fAnalysisLabel->Data();
483  // final results for reference flow:
484  TString intFlowName = "fIntFlowFQD";
485  intFlowName += fAnalysisLabel->Data();
486  // sigma^2:
487  TString sigma2Name = "fSigma2";
488  sigma2Name += fAnalysisLabel->Data();
489  // chi^2:
490  TString chi2Name = "fChi2";
491  chi2Name += fAnalysisLabel->Data();
492  // fitting function:
493  TString fittingFunctionName = "fFittingFunction";
494  fittingFunctionName += fAnalysisLabel->Data();
495 
496  TH1D *qDistribution = NULL;
497  TH1D *sumOfParticleWeights = NULL;
498  TH1D *intFlow[2] = {NULL};
499  TH1D *sigma2[2] = {NULL};
500  TH1D *chi2[2] = {NULL};
501  TF1 *fittingFunction[2] = {NULL};
502 
503  // q-distribution:
504  qDistribution = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s",qDistributionName.Data())));
505  if(qDistribution)
506  {
507  this->SetqDistribution(qDistribution);
508  } else
509  {
510  cout<<"WARNING: qDistribution is NULL in AFAWFQD::GOH() !!!!"<<endl;
511  }
512  // sum of particle weights:
513  sumOfParticleWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s",sumOfParticleWeightsName.Data())));
514  if(sumOfParticleWeights)
515  {
516  this->SetSumOfParticleWeights(sumOfParticleWeights);
517  } else
518  {
519  cout<<"WARNING: sumOfParticleWeights is NULL in AFAWFQD::GOH() !!!!"<<endl;
520  }
521  // final results:
522  for(Int_t f=0;f<2;f++)
523  {
524  // final results for reference flow:
525  intFlow[f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",intFlowName.Data(),sigmaFlag[f].Data())));
526  if(intFlow[f])
527  {
528  this->SetIntFlow(intFlow[f],f);
529  } else
530  {
531  cout<<"WARNING: intFlow[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
532  cout<<"f = "<<f<<endl;
533  }
534  // sigma^2:
535  sigma2[f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",sigma2Name.Data(),sigmaFlag[f].Data())));
536  if(sigma2[f])
537  {
538  this->SetSigma2(sigma2[f],f);
539  } else
540  {
541  cout<<"WARNING: sigma2[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
542  cout<<"f = "<<f<<endl;
543  }
544  // chi^2:
545  chi2[f] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("%s, %s",chi2Name.Data(),sigmaFlag[f].Data())));
546  if(chi2[f])
547  {
548  this->SetChi2(chi2[f],f);
549  } else
550  {
551  cout<<"WARNING: chi2[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
552  cout<<"f = "<<f<<endl;
553  }
554  // fitting functions:
555  fittingFunction[f] = dynamic_cast<TF1*>(outputListHistos->FindObject(Form("%s, %s",fittingFunctionName.Data(),sigmaFlag[f].Data())));
556  if(fittingFunction[f])
557  {
558  this->SetFittingFunction(fittingFunction[f],f);
559  } else
560  {
561  cout<<"WARNING: fittingFunction[f] is NULL in AFAWFQD::GOH() !!!!"<<endl;
562  cout<<"f = "<<f<<endl;
563  }
564  } // end of for(Int_t f=0;f<2;f++)
565 
566  // 5.) fitting parameters:
567  // q-distribution:
568  TString fittingParametersName = "fFittingParameters";
569  fittingParametersName += fAnalysisLabel->Data();
570  TProfile *fittingParameters = NULL;
571  fittingParameters = dynamic_cast<TProfile*>(outputListHistos->FindObject(fittingParametersName.Data()));
572  if(fittingParameters)
573  {
574  this->SetFittingParameters(fittingParameters);
575  } else
576  {
577  cout<<"WARNING:fittingParameters is NULL in AFAWFQD::GOH() !!!!"<<endl;
578  }
579 
580  } else // to if(outputListHistos)
581  {
582  cout<<"WARNING: outputListHistos is NULL in AFAWFQD::GOH() !!!!"<<endl;
583  exit(0);
584  }
585 
586 } // end of void AliFlowAnalysisWithFittingQDistribution::GetOutputHistograms(TList *outputListHistos)
587 
588 //================================================================================================================
589 
591 {
592  //store the final results in output .root file
593  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
594  //output->WriteObject(fHistList, "cobjFQD","SingleKey");
595  fHistList->SetName("cobjFQD");
596  fHistList->SetOwner(kTRUE);
597  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
598  delete output;
599 }
600 
601 //================================================================================================================
602 
604 {
605  //store the final results in output .root file
606  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
607  //output->WriteObject(fHistList, "cobjFQD","SingleKey");
608  fHistList->SetName("cobjFQD");
609  fHistList->SetOwner(kTRUE);
610  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
611  delete output;
612 }
613 
614 //================================================================================================================
615 
617 {
618  //store the final results in output .root file
619  fHistList->SetName("cobjFQD");
620  fHistList->SetOwner(kTRUE);
621  outputFileName->Add(fHistList);
622  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
623 }
624 
625 //================================================================================================================
626 
628 {
629  // Initialize all arrays.
630 
631  for(Int_t s2F=0;s2F<2;s2F++) // sigma^2 not fitted (0) or fitted (1)
632  {
633  fIntFlow[s2F] = NULL;
634  fSigma2[s2F] = NULL;
635  fChi2[s2F] = NULL;
636  fFittingFunction[s2F] = NULL;
637  }
638 
639 } // end of void AliFlowAnalysisWithFittingQDistribution::InitializeArrays()
640 
641 //================================================================================================================
642 
644 {
645  // Book common histograms.
646 
647  // common control histogram:
648  TString commonHistName = "AliFlowCommonHistFQD";
649  commonHistName += fAnalysisLabel->Data();
650  fCommonHists = new AliFlowCommonHist(commonHistName.Data(),commonHistName.Data(),fBookOnlyBasicCCH);
651  fHistList->Add(fCommonHists);
652 
653  // common histograms for final results:
654  TString commonHistResName = "AliFlowCommonHistResultsFQD";
655  commonHistResName += fAnalysisLabel->Data();
656  fCommonHistsResults = new AliFlowCommonHistResults(commonHistResName.Data(),"",fHarmonic);
658 
659 } // end of void AliFlowAnalysisWithFittingQDistribution::BookCommonHistograms()
660 
661 //================================================================================================================
662 
664 {
665  // Book and fill histograms which hold phi, pt and eta weights.
666 
667  if(!fWeightsList)
668  {
669  cout<<"WARNING: fWeightsList is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
670  exit(0);
671  }
672 
673  TString fUseParticleWeightsName = "fUseParticleWeightsFQD";
674  fUseParticleWeightsName += fAnalysisLabel->Data();
675  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
676  fUseParticleWeights->SetLabelSize(0.08);
677  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
678  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
679  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
684  // phi weights:
685  if(fUsePhiWeights)
686  {
687  if(fWeightsList->FindObject("phi_weights"))
688  {
689  fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
690  if(!fPhiWeights){printf("\n WARNING (FQD): !fPhiWeights !!!!\n");exit(0);}
691  if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
692  {
693  cout<<endl;
694  cout<<"WARNING (FQD): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
695  cout<<endl;
696  //exit(0);
697  }
698  } else
699  {
700  cout<<"WARNING: fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
701  exit(0);
702  }
703  } // end of if(fUsePhiWeights)
704  // pt weights:
705  if(fUsePtWeights)
706  {
707  if(fWeightsList->FindObject("pt_weights"))
708  {
709  fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
710  if(!fPtWeights){printf("\n WARNING (FQD): !fPtWeights !!!!\n");exit(0);}
711  if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
712  {
713  cout<<endl;
714  cout<<"WARNING (FQD): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
715  cout<<endl;
716  //exit(0);
717  }
718  } else
719  {
720  cout<<"WARNING: fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
721  exit(0);
722  }
723  } // end of if(fUsePtWeights)
724  // eta weights:
725  if(fUseEtaWeights)
726  {
727  if(fWeightsList->FindObject("eta_weights"))
728  {
729  fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
730  if(!fEtaWeights){printf("\n WARNING (FQD): !fEtaWeights !!!!\n");exit(0);}
731  if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
732  {
733  cout<<endl;
734  cout<<"WARNING (FQD): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
735  cout<<endl;
736  //exit(0);
737  }
738  } else
739  {
740  cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWFQD::BAFWH() !!!!"<<endl;
741  exit(0);
742  }
743  } // end of if(fUseEtaWeights)
744 
745 } // end of AliFlowAnalysisWithFittingQDistribution::BookAndFillWeightsHistograms()
746 
747 //================================================================================================================================
748 
750 {
751  // access needed common constants from AliFlowCommonConstants
752 
765 
766 } // end of void AliFlowAnalysisWithFittingQDistribution::AccessConstants()
767 
768 //================================================================================================================================
769 
771 {
772  // Book all objects relevant for fitting of q-distributions.
773 
774  // Flags:
775  TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
776  // q-distribution:
777  TString fqDistributionName = "fqDistribution";
778  fqDistributionName += fAnalysisLabel->Data();
779  fqDistribution = new TH1D(Form("%s",fqDistributionName.Data()),"q-distribution",fqNbins,fqMin,fqMax);
780  fqDistribution->SetXTitle(Form("q_{%d}=|Q_{%d}|/#sqrt{M}",fHarmonic,fHarmonic));
781  fqDistribution->SetYTitle("Counts");
782  fHistList->Add(fqDistribution);
783  // q-distribution vs multiplicity:
785  {
786  TString fqDistributionVsMultName = "fqDistributionVsMult";
787  fqDistributionVsMultName += fAnalysisLabel->Data();
788  fqDistributionVsMult = new TH2D(Form("%s",fqDistributionVsMultName.Data()),"q-distribution vs M",fnBinsMult,fMinMult,fMaxMult,fqNbins,fqMin,fqMax);
789  fqDistributionVsMult->GetYaxis()->SetTitle(Form("q_{%d}=|Q_{%d}|/#sqrt{M}",fHarmonic,fHarmonic));
791  {
792  fqDistributionVsMult->GetXaxis()->SetTitle("# RPs");
794  {
795  fqDistributionVsMult->GetXaxis()->SetTitle("Reference multiplicity (from ESD)");
797  {
798  fqDistributionVsMult->GetXaxis()->SetTitle("# POIs");
799  }
801  } // end of if(fStoreqDistributionVsMult)
802  // Sum of particle weights:
803  TString fSumOfParticleWeightsName = "fSumOfParticleWeights";
804  fSumOfParticleWeightsName += fAnalysisLabel->Data();
805  fSumOfParticleWeights = new TH1D(Form("%s",fSumOfParticleWeightsName.Data()),"Sum of particle weights",10000,0,100000); // (to be improved - harwired limits and number of bins)
806  fSumOfParticleWeights->SetXTitle("#sum_{i=1}^{N} w_{i}");
807  fSumOfParticleWeights->GetXaxis()->SetTitleSize(0.03);
808  fSumOfParticleWeights->SetYTitle("Counts");
810  // Final results for integrated flow:
811  TString fIntFlowName = "fIntFlowFQD";
812  fIntFlowName += fAnalysisLabel->Data();
813  // sigma^2:
814  TString fSigma2Name = "fSigma2";
815  fSigma2Name += fAnalysisLabel->Data();
816  // chi^2:
817  TString fChi2Name = "fChi2";
818  fChi2Name += fAnalysisLabel->Data();
819  // fitting function:
820  TString fittingFunctionName = "fFittingFunction";
821  fittingFunctionName += fAnalysisLabel->Data();
822 
823  for(Int_t f=0;f<2;f++) // sigma^2 not fitted (0) or fitted (1)
824  {
825  // final results for integrated flow:
826  fIntFlow[f] = new TH1D(Form("%s, %s",fIntFlowName.Data(),sigmaFlag[f].Data()),Form("Reference Flow (%s)",sigmaFlag[f].Data()),1,0,1);
827  fIntFlow[f]->SetLabelSize(0.08);
828  fHistList->Add(fIntFlow[f]);
829  // sigma^2:
830  fSigma2[f] = new TH1D(Form("%s, %s",fSigma2Name.Data(),sigmaFlag[f].Data()),Form("#sigma^{2} (%s)",sigmaFlag[f].Data()),1,0,1);
831  fSigma2[f]->SetLabelSize(0.08);
832  (fSigma2[f]->GetXaxis())->SetBinLabel(1,"#sigma^{2}");
833  fHistList->Add(fSigma2[f]);
834  // chi^2:
835  fChi2[f] = new TH1D(Form("%s, %s",fChi2Name.Data(),sigmaFlag[f].Data()),Form("#chi^{2} (%s)",sigmaFlag[f].Data()),1,0,1);
836  fChi2[f]->SetLabelSize(0.08);
837  (fChi2[f]->GetXaxis())->SetLabelOffset(0.01);
838  (fChi2[f]->GetXaxis())->SetBinLabel(1,"#chi^{2}");
839  fHistList->Add(fChi2[f]);
840  // fitting functions:
841  fFittingFunction[f] = new TF1(Form("%s, %s",fittingFunctionName.Data(),sigmaFlag[f].Data()),"[2]*(x/[1])*exp(-(x*x+[0]*[0])/(2.*[1]))*TMath::BesselI0(x*[0]/[1])");
842  fHistList->Add(fFittingFunction[f]);
843  } // end of for(Int_t f=0;f<2;f++) // sigma^2 not fitted or fitted
844  // Book profile fFittingParameters which will hold all fitting parameters:
845  TString fFittingParametersName = "fFittingParameters";
846  fFittingParametersName += fAnalysisLabel->Data();
847  fFittingParameters = new TProfile(fFittingParametersName.Data(),"Parameters for fitting q-distribution",13,0,13);
848  fFittingParameters->SetLabelSize(0.05);
849  fFittingParameters->GetXaxis()->SetBinLabel(1,"treshold");
850  fFittingParameters->GetXaxis()->SetBinLabel(2,Form("starting v_{%d}",fHarmonic));
851  fFittingParameters->GetXaxis()->SetBinLabel(3,Form("min. v_{%d}",fHarmonic));
852  fFittingParameters->GetXaxis()->SetBinLabel(4,Form("max. v_{%d}",fHarmonic));
853  fFittingParameters->GetXaxis()->SetBinLabel(5,"starting #sigma^{2}");
854  fFittingParameters->GetXaxis()->SetBinLabel(6,"min. #sigma^{2}");
855  fFittingParameters->GetXaxis()->SetBinLabel(7,"max. #sigma^{2}");
856  fFittingParameters->GetXaxis()->SetBinLabel(8,"Final result is from #sigma^{2} fitted?");
857  fFittingParameters->GetXaxis()->SetBinLabel(9,"Print results on the screen?");
858  fFittingParameters->GetXaxis()->SetBinLabel(10,"fStoreqDistributionVsMult");
859  fFittingParameters->GetXaxis()->SetBinLabel(11,"fMultiplicityIs");
860  fFittingParameters->GetXaxis()->SetBinLabel(12,"fDoFit");
861  fFittingParameters->GetXaxis()->SetBinLabel(13,"fExactNoRPs");
863 
864 } // end of void AliFlowAnalysisWithFittingQDistribution::BookEverythingForDistributions()
865 
866 //================================================================================================================================
867 
869 {
870  // Do the final fit of q-distribution.
871 
872  Int_t s2F = (Int_t)(sigma2Fitted); // shortcut
873  Double_t AvM = fSumOfParticleWeights->GetMean(1); // average multiplicity
874  //Int_t nEvts = (Int_t)fSumOfParticleWeights->GetEntries(); // number of events:
875 
876  // Start fitting from the bin with at least fTreshold entries,
877  // finish fitting at the bin with at least fTreshold entries:
878  Int_t binMin = fqDistribution->FindFirstBinAbove(fTreshold);
879  Int_t binMax = fqDistribution->FindLastBinAbove(fTreshold);
880  Double_t binWidth = fqDistribution->GetBinWidth(4); // assuming that all bins have the same width
881  if(TMath::Abs(binWidth) < 1.e-44)
882  {
883  cout<<endl;
884  cout<<"WARNING (FQD): binWidth == 0 in AFAWFQD::DoFit()"<<endl;
885  cout<<endl;
886  exit(0);
887  }
888  Double_t qmin = (binMin-1)*binWidth;
889  Double_t qmax = (binMax)*binWidth;
890  Double_t ent = 0.; // number of entries between binMin and binMax:
891  for(Int_t b=binMin;b<=binMax;b++)
892  {
893  ent += fqDistribution->GetBinContent(b);
894  }
895 
896  Double_t norm = binWidth*ent; // norm (assuming that all bins have the same width)
897  // Fitting function:
898  fFittingFunction[s2F]->SetRange(qmin,qmax);
899  fFittingFunction[s2F]->SetParNames("v*sqrt{sum of particle weights}","sigma^2","norm");
900  fFittingFunction[s2F]->SetParameters(fvStart*pow(AvM,0.5),fSigma2Start,norm);
901  fFittingFunction[s2F]->SetParLimits(0,fvMin*pow(AvM,0.5),fvMax*pow(AvM,0.5));
902  if(s2F == 0)
903  {
904  fFittingFunction[s2F]->FixParameter(1,0.5);
905  } else
906  {
907  fFittingFunction[s2F]->SetParLimits(1,fSigma2Min,fSigma2Max);
908  }
909  fFittingFunction[s2F]->FixParameter(2,norm);
910  // Fitting is done here:
911  // Remark: do fit only if # of entries > 50 - this is only a pragmatics fix to avoid TMinuit crash (to be improved)
912  if(ent > 50)
913  {
914  fqDistribution->Fit(fFittingFunction[s2F]->GetName(),"NQ","",qmin,qmax);
915  }
916 
917  // Final results:
918  Double_t v = 0.; // reference flow
919  Double_t vError = 0.; // error of reference flow
920  Double_t sigma2 = 0.; // sigma^2
921  Double_t sigma2Error = 0.; // error of sigma^2
922  Double_t chi2 = 0; // chi^2 from Minuit
923  // Reference flow:
924  if(AvM)
925  {
926  v = fFittingFunction[s2F]->GetParameter(0)/pow(AvM,0.5);
927  vError = fFittingFunction[s2F]->GetParError(0)/pow(AvM,0.5);
928  fIntFlow[s2F]->SetBinContent(1,v); // s2F is shortcut for "sigma^2 fitted"
929  fIntFlow[s2F]->SetBinError(1,vError); // s2F is shortcut for "sigma^2 fitted"
930  } else
931  {
932  cout<<endl;
933  cout<<"WARNING (FQD): AvM == 0 in AFAWFQD::DoFit()"<<endl;
934  cout<<endl;
935  }
936  // sigma^2::
937  if(s2F == 0) // sigma^2 not fitted, but fixed to 0.5
938  {
939  sigma2 = 0.5;
940  fSigma2[0]->SetBinContent(1,sigma2);
941  fSigma2[0]->SetBinError(1,0.);
942  } else // sigma^2 fitted
943  {
944  sigma2 = fFittingFunction[s2F]->GetParameter(1);
945  sigma2Error = fFittingFunction[s2F]->GetParError(1);
946  fSigma2[1]->SetBinContent(1,sigma2);
947  fSigma2[1]->SetBinError(1,sigma2Error);
948  }
949  // chi^2:
950  chi2 = fFittingFunction[s2F]->GetChisquare();
951  fChi2[s2F]->SetBinContent(1,chi2);
952 
953 } // end of void AliFlowAnalysisWithFittingQDistribution::DoFit()
954 
955 //================================================================================================================================
956 
958 {
959  // Fill common result histogram for reference flow and resolution.
960 
961  // Remark: by default the result obtained with sigma^2 fitted is being stored.
962  // In order to store the result obtained with sigma^2 fixed use a setter SetFinalResultIsFromSigma2Fitted(kFALSE).
963 
964  Int_t s2F = (Int_t)sigma2Fitted; // shortcut
965 
966  // Reference flow:
967  Double_t v = fIntFlow[s2F]->GetBinContent(1);
968  Double_t vError = fIntFlow[s2F]->GetBinError(1);
970  // Resolution:
971  Double_t AvM = fSumOfParticleWeights->GetMean(1);
972  Double_t chi2 = AvM*pow(v,2.); // chi^2
973  if(chi2>=0.)
974  {
975  fCommonHistsResults->FillChi(pow(chi2,0.5));
976  }
977 
978 } // end of void AliFlowAnalysisWithFittingQDistribution::FillCommonHistResultsIntFlow(Bool_t useParticleWeights, Bool_t sigma2NotFixed)
979 
980 //================================================================================================================================
981 
983 {
984  // Print the final results on the screen.
985 
986  // shortcut for the harmonic:
987  Int_t n = 0;
989  {
990  n = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1);
991  }
992 
993  // printing:
994  cout<<endl;
995  cout<<"***************************************"<<endl;
996  cout<<"***************************************"<<endl;
997  cout<<" reference flow by fitting "<<endl;
998  cout<<" q-distribution: "<<endl;
1000  {
1001  cout<<" (with weights) "<<endl;
1002  } else
1003  {
1004  cout<<" (without weights) "<<endl;
1005  }
1006  cout<<endl;
1007  cout<<"1.) sigma^2 not fitted: "<<endl;
1008  cout<<" v_"<<n<<"{FQD} = "<<fIntFlow[0]->GetBinContent(1)<<" +/- "<<fIntFlow[0]->GetBinError(1)<<endl;
1009  cout<<" sigma^2 = 0.5 +/- 0 "<<endl;
1010  cout<<" chi^2 = "<<fChi2[0]->GetBinContent(1)<<" (Minuit)"<<endl;
1011  cout<<endl;
1012  if(TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMin)<1.e-10 ||
1013  TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMax)<1.e-10)
1014  {
1015  cout<<" WARNING: value of v_"<<n<<"{FQD}"<<" is on the boundary"<<endl;
1016  cout<<" of fitting interval. Redo the fit."<< endl;
1017  cout<<endl;
1018  }
1019  cout<<"2.) sigma^2 fitted: "<<endl;
1020  cout<<" v_"<<n<<"{FQD} = "<<fIntFlow[1]->GetBinContent(1)<<" +/- "<<fIntFlow[1]->GetBinError(1)<<endl;
1021  cout<<" sigma^2 = "<<fSigma2[1]->GetBinContent(1)<<" +/- "<<fSigma2[1]->GetBinError(1)<<endl;
1022  cout<<" chi^2 = "<<fChi2[1]->GetBinContent(1)<<" (Minuit)"<<endl;
1023  cout<<endl;
1024  if(TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMin)<1.e-10 ||
1025  TMath::Abs(fIntFlow[0]->GetBinContent(1)-fvMax)<1.e-10)
1026  {
1027  cout<<" WARNING: value of v_"<<n<<"{FQD}"<<" is on the boundary"<<endl;
1028  cout<<" of fitting interval. Redo the fit."<< endl;
1029  cout<<endl;
1030  }
1031  if(TMath::Abs(fSigma2[1]->GetBinContent(1)-fSigma2Min)<1.e-10 ||
1032  TMath::Abs(fSigma2[1]->GetBinContent(1)-fSigma2Max)<1.e-10)
1033  {
1034  cout<<" WARNING: value of sigma^2 is on the boundary"<<endl;
1035  cout<<" of fitting interval. Redo the fit."<< endl;
1036  cout<<endl;
1037  }
1038  cout<<" nEvts = "<<fSumOfParticleWeights->GetEntries()<<", AvM = "<<fSumOfParticleWeights->GetMean()<<endl;
1039  cout<<endl;
1040  cout<<"***************************************"<<endl;
1041  cout<<"***************************************"<<endl;
1042  cout<<endl;
1043 
1044 } // end of void AliFlowAnalysisWithFittingQDistribution::PrintOnTheScreen()
1045 
1046 //================================================================================================================================
1047 
1049 {
1050  // Store fitting parameters for the fit of q-distribution in profile fFittingParameters.
1051 
1052  if(!fFittingParameters)
1053  {
1054  cout<<endl;
1055  cout<<"WARNING (FQD): fFittingParameters is NULL in AFAWFQD::SFP() !!!!"<<endl;
1056  cout<<endl;
1057  exit(0);
1058  }
1059  fFittingParameters->Reset();
1060  fFittingParameters->Fill(0.5,fTreshold);
1061  fFittingParameters->Fill(1.5,fvStart);
1062  fFittingParameters->Fill(2.5,fvMin);
1063  fFittingParameters->Fill(3.5,fvMax);
1064  fFittingParameters->Fill(4.5,fSigma2Start);
1065  fFittingParameters->Fill(5.5,fSigma2Min);
1066  fFittingParameters->Fill(6.5,fSigma2Max);
1070  // which multiplicity was used:
1071  if(fMultiplicityIs==AliFlowCommonConstants::kRP) // # of Reference Particles
1072  {
1073  fFittingParameters->Fill(10.5,0); // 0 = # of Reference Particles
1075  {
1076  fFittingParameters->Fill(10.5,1); // 1 = ref. mult. from ESD
1078  {
1079  fFittingParameters->Fill(10.5,2); // 2 = # of Particles of Interest
1080  }
1081  fFittingParameters->Fill(11.5,fDoFit);
1082  fFittingParameters->Fill(12.5,fExactNoRPs);
1083 
1084 } // end of void AliFlowAnalysisWithFittingQDistribution::StoreFittingParameters()
1085 
1086 //================================================================================================================================
1087 
1089 {
1090  // Access fitting parameters for the fit of q-distribution.
1091 
1092  fTreshold = fFittingParameters->GetBinContent(1);
1093  fvStart = fFittingParameters->GetBinContent(2);
1094  fvMin = fFittingParameters->GetBinContent(3);
1095  fvMax = fFittingParameters->GetBinContent(4);
1096  fSigma2Start = fFittingParameters->GetBinContent(5);
1097  fSigma2Min = fFittingParameters->GetBinContent(6);
1098  fSigma2Max = fFittingParameters->GetBinContent(7);
1100  fPrintOnTheScreen = (Bool_t)fFittingParameters->GetBinContent(9);
1101  fStoreqDistributionVsMult = (Bool_t)fFittingParameters->GetBinContent(10);
1102  fDoFit = (Bool_t)fFittingParameters->GetBinContent(12);
1103  fExactNoRPs = (Int_t)fFittingParameters->GetBinContent(13);
1104 
1105 } // end of void AliFlowAnalysisWithFittingQDistribution::AccessFittingParameters()
double Double_t
Definition: External.C:58
AliFlowTrackSimple * GetTrack(Int_t i)
Int_t GetNumberOfRPs() const
void SetCommonHistsResults(AliFlowCommonHistResults *const chr)
void SetIntFlow(TH1D *const intFlow, Int_t sigmaFitted)
Bool_t InRPSelection() const
Bool_t FillControlHistograms(AliFlowEventSimple *anEvent, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
Double_t Phi() const
int Int_t
Definition: External.C:63
Definition: External.C:228
Definition: External.C:212
static AliFlowCommonConstants * GetMaster()
Int_t GetNumberOfPOIs(Int_t i=1) const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Bool_t FillIntegratedFlow(Double_t aV, Double_t anError)
Int_t GetReferenceMultiplicity() const
Double_t Pt() const
Double_t Eta() const
bool Bool_t
Definition: External.C:53
TProfile * GetHarmonic()
void SetSigma2(TH1D *const sigma2, Int_t sigmaFitted)
Definition: External.C:196
Int_t NumberOfTracks() const