]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bug fix
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithFittingQDistribution.cxx
CommitLineData
ce4a88f5 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"
ee0860e8 37#include "TLegend.h"
ce4a88f5 38#include "TParticle.h"
39#include "TProfile.h"
40#include "AliFlowEventSimple.h"
41#include "AliFlowTrackSimple.h"
42#include "AliFlowAnalysisWithFittingQDistribution.h"
43
ee0860e8 44
ce4a88f5 45class TH1;
46class TGraph;
47class TPave;
48class TLatex;
49class TMarker;
50class TObjArray;
51class TList;
52class TCanvas;
53class TSystem;
54class TROOT;
55class AliFlowVector;
56class TVector;
57
58//================================================================================================================
59
60ClassImp(AliFlowAnalysisWithFittingQDistribution)
61
62AliFlowAnalysisWithFittingQDistribution::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),
ee0860e8 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)
ce4a88f5 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
123AliFlowAnalysisWithFittingQDistribution::~AliFlowAnalysisWithFittingQDistribution()
124{
125 // desctructor
126 delete fHistList;
127}
128
129
130//================================================================================================================
131
132
133void 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
ee0860e8 145 // store fitting parameters:
146 this->StoreFittingParameters();
147
ce4a88f5 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
162void 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
260void 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++)
ee0860e8 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
ce4a88f5 381 } else // to if(outputListHistos)
382 {
383 cout<<"WARNING: outputListHistos is NULL in AFAWFQD::GOH() !!!!"<<endl;
384 exit(0);
385 }
ee0860e8 386
ce4a88f5 387
388} // end of void AliFlowAnalysisWithFittingQDistribution::GetOutputHistograms(TList *outputListHistos)
389
390
391//================================================================================================================
392
393
394void AliFlowAnalysisWithFittingQDistribution::Finish(Bool_t doFit)
395{
396 // calculate the final results
397
ee0860e8 398 // a) acces the constants and fitting paremeters;
ce4a88f5 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
ee0860e8 404 // access the constants and fitting paremeters:
ce4a88f5 405 this->AccessConstants();
ee0860e8 406 this->AccessFittingParameters();
ce4a88f5 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
ee0860e8 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
ce4a88f5 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
455void 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
468void 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
483void 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
504void 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
526void 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
608void 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
631void 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
ee0860e8 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
ce4a88f5 694} // end of void AliFlowAnalysisWithFittingQDistribution::BookEverythingForDistributions()
695
696
697//================================================================================================================================
698
699
700void 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);
ee0860e8 706 Int_t s2NF = 1-(Int_t)(sigma2NotFixed);
ce4a88f5 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
ee0860e8 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)
ce4a88f5 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");
ee0860e8 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
d6d06a8e 748 if(s2NF == 0)
ce4a88f5 749 {
ee0860e8 750 fittingFun->SetParLimits(1,fSigma2Min,fSigma2Max); // to be improved (add setters for sigma^2_min and sigma^2_max)
ce4a88f5 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 }
d6d06a8e 776 if(s2NF == 0)
ce4a88f5 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 {
ee0860e8 785 sigma2 = 0.5;
786 fSigma2[pW][1]->SetBinContent(1,sigma2);
ce4a88f5 787 fSigma2[pW][1]->SetBinError(1,0.);
788 }
789
ee0860e8 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
ce4a88f5 812} // end of void AliFlowAnalysisWithFittingQDistribution::DoFit(Bool_t useParticleWeights)
813
814
815//================================================================================================================================
816
817
818void 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
ee0860e8 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
ce4a88f5 841 if(!(fCommonHistsResults))
842 {
843 cout<<"WARNING: fCommonHistsResults is NULL in AFAWFQD::FCHRIF() !!!!"<<endl;
844 exit(0);
845 }
846
ee0860e8 847 // fill integrated flow:
ce4a88f5 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
ee0860e8 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
ce4a88f5 862} // end of void AliFlowAnalysisWithFittingQDistribution::FillCommonHistResultsIntFlow(Bool_t useParticleWeights, Bool_t sigma2NotFixed)
863
864
865//================================================================================================================================
866
867
868void 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//================================================================================================================================
ee0860e8 946
947
948void 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
983void 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()