]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TPC/AliPerfAnalyzeInvPt.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity fix
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliPerfAnalyzeInvPt.cxx
CommitLineData
ab886e13 1
2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
ba06aaec 17//------------------------------------------------------------------------------
18// Implementation of AliPerfAnalyzeInvPt class. It analyzes the output of
ab886e13 19// AliPerformancePtCalib.cxx and AliPerformancePtCalibMC.cxx: via
20// AliPerformancePtCalib::Analyse() or AliPerformancePtCalibMC::Analyse()
21//
22// Projection of charge/pt vs theta and vs phi resp. histoprams will be fitted with either
23// polynomial or gaussian fit function to extract minimum position of charge/pt.
ba06aaec 24// Fit options and theta, phi bins can be set by user.
25// Attention: use the Set* functions of AliPerformancePtCalib.h and AliPerformancePtCalibMC.h
26// when running AliPerformancePtCalib*::Analyse()
27// The result of the analysis (histograms/graphs) are stored in the folder which is
28// a data member of AliPerformancePtCalib*.
29//
30// Author: S.Schuchmann 11/13/2009
31//------------------------------------------------------------------------------
32
33/*
34
ab886e13 35// after running the performance task of AliPerformancePtCalib* or AliPerformancePtCalibMC*
36// read the file, and get component:
ba06aaec 37
38TFile f("Output.root");
ab886e13 39AliPerformancePtCalib *compObj = (AliPerformancePtCalib*)coutput->FindObject("AliPerformancePtCalib");
40//or
41AliPerformancePtCalibMC *compObj = (AliPerformancePtCalibMC*)coutput->FindObject("AliPerformancePtCalibMC");
42
43// set phi and theta bins for fitting and analyse comparison data
44compObj->SetProjBinsTheta(thetaBins,nThetaBins,minPhi, maxPhi);
45compObj->SetProjBinsPhi(phiBins,nPhiBins,minTheta,maxTheta);
46compObj->SetMakeFitOption(kFALSE,exclRange,fitRange);
47compObj->SetDoRebin(rebin);
ba06aaec 48compObj->Analyse();
ab886e13 49//for details see functions of class AliPerformancePtCalib*
ba06aaec 50
51// the output histograms/graphs will be stored in the folder "folderRes"
52compObj->GetAnalysisFolder()->ls("*");
53
54// user can save whole comparison object (or only folder with anlysed histograms)
55// in the seperate output file (e.g.)
56TFile fout("Analysed_InvPt.root","recreate");
57compObj->Write(); // compObj->GetAnalysisFolder()->Write();
58fout.Close();
59
60*/
61
acb9d358 62#include "TNamed.h"
ba06aaec 63#include "TF1.h"
acb9d358 64#include "TH1F.h"
65#include "TH2F.h"
ba06aaec 66#include "TH1D.h"
ab886e13 67#include "TMath.h"
acb9d358 68#include "TGraphErrors.h"
ba06aaec 69#include "TCanvas.h"
70#include "TObjArray.h"
acb9d358 71
acb9d358 72#include "AliPerfAnalyzeInvPt.h"
73
74ClassImp(AliPerfAnalyzeInvPt)
75
76
77// fit functions
78//____________________________________________________________________________________________________________________________________________
ba06aaec 79 Double_t AliPerfAnalyzeInvPt::Polynomial(Double_t *x, const Double_t *par)
acb9d358 80{
ba06aaec 81 // fit function for fitting of 1/pt with a polynomial of 4th order, rejecting points between +/-par[4]
acb9d358 82
6bcacaf1 83
acb9d358 84 if (x[0] > -par[4] && x[0] < par[4]) {
85 TF1::RejectPoint();
86 return 0;
87 }
6bcacaf1 88
acb9d358 89 return par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4);
6bcacaf1 90
acb9d358 91}
92//____________________________________________________________________________________________________________________________________________
ba06aaec 93Double_t AliPerfAnalyzeInvPt::PolynomialRejP(Double_t *x, const Double_t *par)
acb9d358 94{
ba06aaec 95 // fit function for fitting of 1/pt with a polynomial of 4th order to improve result (fit range and rejection zone is adjusted to first guess of minimum position)
6bcacaf1 96
acb9d358 97 Double_t pos = par[5];
98 Double_t neg = - par[5];
99 pos += par[4];
100 neg += par[4];
6bcacaf1 101
acb9d358 102 if (x[0] > neg && x[0] < pos) {
103 TF1::RejectPoint();
104 return 0;
105 }
6bcacaf1 106
acb9d358 107 return par[2]+par[0]*pow((x[0]-par[1]),2)+par[3]*pow((x[0]-par[1]),4);
6bcacaf1 108
acb9d358 109}
110//____________________________________________________________________________________________________________________________________________
ba06aaec 111Double_t AliPerfAnalyzeInvPt::InvGauss(Double_t *x, const Double_t *par)
acb9d358 112{
ba06aaec 113 // fit function for fitting of 1/pt with gaussian, rejecting points between +/-par[4]
acb9d358 114 if (x[0] > -par[6] && x[0] < par[6]) {
115 TF1::RejectPoint();
116 return 0;
117 }
118
ba06aaec 119 return par[3]+par[0]*(TMath::Exp(-0.5*(TMath::Power((x[0]-par[1])/par[2], 2.0)))+par[4]*pow((x[0]-par[1]),2)+par[5]*pow((x[0]-par[1]),4)) ;
acb9d358 120}
121//____________________________________________________________________________________________________________________________________________
ba06aaec 122Double_t AliPerfAnalyzeInvPt::InvGaussRejP(Double_t *x, const Double_t *par)
acb9d358 123{
ba06aaec 124 // fit function for fitting of 1/pt with gaussian to improve result (fit range and rejection zone is adjusted to first guess of minimum position)
acb9d358 125 Double_t pos = par[7];//0.12;
126 Double_t neg = - par[7];//0.12;
127 pos += par[6];
128 neg += par[6];
129
130 if (x[0] > neg && x[0] < pos) {
131 TF1::RejectPoint();
132 return 0;
133 }
134
135
136 return par[3]+par[0]*TMath::Exp(-0.5*(TMath::Power((x[0]-par[1])/par[2], 2.0)))+par[4]*pow((x[0]-par[1]),2)+par[5]*pow((x[0]-par[1]),4) ;
137}
138
139
140//_____________________________________________________________________________________________________________________________________________
141AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt():
142 TNamed("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt"),
67354362 143 fThetaBins(),
144 fPhiBins(),
acb9d358 145 fNThetaBins(0),
146 fNPhiBins(0),
147 fRange(0),
148 fExclRange(0),
149 fFitGaus(0) ,
6bcacaf1 150 fDoRebin(0),
151 fRebin(0),
67354362 152 fHistFitTheta(),
153 fHistFitPhi(),
ab886e13 154 //histograms,graphs and functions
acb9d358 155 fHistH2InvPtTheta(0),
156 fHistH2InvPtPhi(0),
157 fGrMinPosTheta(0),
158 fGrMinPosPhi(0),
159 fFitMinPos(0),
160 fFitMinPosRejP(0),
161 fFitInvGauss(0),
162 fFitInvGaussRejP(0)
163{
164 // Default constructor
165
ab886e13 166 fFitGaus = kFALSE;// flag for gaussian fit
167 fNThetaBins = 0; // theta bins for projections
168 fNPhiBins = 0; //phi bins for projections
169 fRange = 0; //fit range
170 fExclRange = 0; //fit exclusion range
171 fDoRebin = kFALSE; // flag for rebin
172 fRebin = 0; // bins for rebin
acb9d358 173
acb9d358 174 for(Int_t i=0;i<100;i++){
175
176 fHistFitTheta[i] = NULL;
177 fHistFitPhi[i] = NULL;
178 }
ba06aaec 179
180
acb9d358 181}
182//_____________________________________________________________________________________________________________________________________________
c6d45264 183AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt(const Char_t* name="AliAnalyzeInvPt",const Char_t* title="AliAnalyzeInvPt"):
ba06aaec 184 TNamed(name, title),
67354362 185 fThetaBins(),
186 fPhiBins(),
acb9d358 187 fNThetaBins(0),
188 fNPhiBins(0),
189 fRange(0),
190 fExclRange(0),
191 fFitGaus(0) ,
6bcacaf1 192 fDoRebin(0),
193 fRebin(0),
67354362 194 fHistFitTheta(),
195 fHistFitPhi(),
ab886e13 196 //histograms,graphs and functions
acb9d358 197 fHistH2InvPtTheta(0),
198 fHistH2InvPtPhi(0),
199 fGrMinPosTheta(0),
200 fGrMinPosPhi(0),
201 fFitMinPos(0),
202 fFitMinPosRejP(0),
203 fFitInvGauss(0),
204 fFitInvGaussRejP(0)
ba06aaec 205{
ab886e13 206 fFitGaus = kFALSE;// flag for gaussian fit
207 fNThetaBins = 0; // theta bins for projections
208 fNPhiBins = 0; //phi bins for projections
209 fRange = 0; //fit range
210 fExclRange = 0; //fit exclusion range
211 fDoRebin = kFALSE; // flag for rebin
212 fRebin = 0; // bins for rebin
ba06aaec 213
ba06aaec 214 for(Int_t i=0;i<100;i++){
acb9d358 215
ba06aaec 216 fHistFitTheta[i] = NULL;
217 fHistFitPhi[i] = NULL;
acb9d358 218 }
ba06aaec 219
220}
acb9d358 221
222
ba06aaec 223//______________________________________________________________________________________________________________________________________
224void AliPerfAnalyzeInvPt::InitGraphs(Double_t *binsXTheta,Double_t *fitParamTheta,Double_t *errFitParamTheta,Double_t *binsXPhi,Double_t *fitParamPhi,Double_t *errFitParamPhi){
acb9d358 225
ba06aaec 226 // initialize graphs for storing fit parameter
acb9d358 227
ba06aaec 228 fGrMinPosTheta = new TGraphErrors(fNThetaBins,binsXTheta,fitParamTheta,0, errFitParamTheta);
229 fGrMinPosTheta->SetMarkerStyle(20);
230 fGrMinPosTheta->SetMarkerColor(2);
231 fGrMinPosTheta->SetLineColor(2);
ab886e13 232 fGrMinPosTheta->GetYaxis()->SetTitle("min pos (GeV/c)^{-1}");
ba06aaec 233 fGrMinPosTheta->GetXaxis()->SetTitle("#theta bin no.");
234 fGrMinPosTheta->GetYaxis()->SetTitleOffset(1.2);
ab886e13 235 fGrMinPosTheta->SetTitle("#theta bins");
ba06aaec 236
237 fGrMinPosPhi = new TGraphErrors(fNPhiBins,binsXPhi,fitParamPhi,0,errFitParamPhi);
238 fGrMinPosPhi->SetMarkerStyle(20);
239 fGrMinPosPhi->SetMarkerColor(4);
240 fGrMinPosPhi->SetLineColor(4);
ab886e13 241 fGrMinPosPhi->GetYaxis()->SetTitle("min pos (GeV/c)^{-1}");
ba06aaec 242 fGrMinPosPhi->GetXaxis()->SetTitle("#phi bin no.");
243 fGrMinPosPhi->GetYaxis()->SetTitleOffset(1.2);
ab886e13 244 fGrMinPosPhi->SetTitle("#phi bins");
acb9d358 245}
246
247//______________________________________________________________________________________________________________________________________
248
249void AliPerfAnalyzeInvPt::InitFitFcn(){
ba06aaec 250 // fit functions
251 fFitMinPos = new TF1("fFitMinPos", Polynomial,-4.0,4.0,5);
252 fFitMinPos->SetLineColor(4);
253 fFitMinPos->SetLineWidth(1);
254 fFitMinPos->SetParameter(0,1.0);
255 fFitMinPos->SetParameter(1,0.0);
256 fFitMinPos->SetParameter(2,0.0);
257
258 fFitMinPosRejP = new TF1("fFitMinPosRejP",PolynomialRejP,-4.0,4.0,6);
259 fFitMinPosRejP->SetLineColor(2);
260 fFitMinPosRejP->SetLineWidth(1);
261 fFitMinPosRejP->SetParameter(0,1.0);
262 fFitMinPosRejP->SetParameter(1,0.0);
263 fFitMinPosRejP->SetParameter(2,0.0);
acb9d358 264
265
ba06aaec 266 fFitInvGauss = new TF1("fFitInvGauss", InvGauss,-4.0,4.0,6);
267 fFitInvGauss->SetLineColor(4);
268 fFitInvGauss->SetLineWidth(1);
269 fFitInvGauss->SetParameter(2,1.0);
acb9d358 270
ba06aaec 271 fFitInvGaussRejP = new TF1("fFitInvGaussRejP", InvGaussRejP,-4.0,4.0,7);
272 fFitInvGaussRejP->SetLineColor(2);
273 fFitInvGaussRejP->SetLineWidth(1);
274 fFitInvGaussRejP->SetParameter(2,1.0);
acb9d358 275
ba06aaec 276}
acb9d358 277//______________________________________________________________________________________________________________________________________
ba06aaec 278void AliPerfAnalyzeInvPt::StartAnalysis(const TH2F *histThetaInvPt, const TH2F *histPhiInvPt, TObjArray* aFolderObj){
ab886e13 279 //start analysis: fitting charge/pt spectra
acb9d358 280
281
282
283 if(!histThetaInvPt) {
284 Printf("warning: no 1/pt histogram to analyse in theta bins!");
18e588e9 285 return;
acb9d358 286 }
287 if(!histPhiInvPt) {
288 Printf("warning: no 1/pt histogram to analyse in phit bins!");
18e588e9 289 return;
acb9d358 290 }
291
ab886e13 292 Double_t thetaBins[9] = {0.77,0.97,1.17,1.37,1.57,1.77,1.97,2.17,2.37}; // default theta bins
293 Double_t phiBins[13] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.5}; // default phi bins
acb9d358 294
295 Int_t nThetaBins = 9;
296 Int_t nPhiBins = 13;
297 Double_t range = 1.0; // fit range
298 Double_t exclRange = 0.13; // range of point rejection in fit
299 Bool_t fitGaus=kFALSE; // fit with Gaussian or with polynomial (kFALSE)
300
ba06aaec 301 if(!fNThetaBins){
acb9d358 302 fNThetaBins = nThetaBins;
303 for(Int_t k = 0;k<fNThetaBins;k++){
304 fThetaBins[k] = thetaBins[k];
305 }
306 Printf("warning: for theta bins no user intput! bins are set to default values.");
307 }
308
ba06aaec 309 if(!fNPhiBins){
acb9d358 310 fNPhiBins = nPhiBins;
311 for(Int_t k = 0;k<fNPhiBins;k++){
312 fPhiBins[k] = phiBins[k];
313 }
314
315 Printf("warning: for phi bins no user intput! bins are set to default values.");
316
317 }
318
ba06aaec 319 if(!fRange){
acb9d358 320 fRange = range;
321 fExclRange = exclRange;
322 fFitGaus = fitGaus;
323 Printf("warning: no fit range is set. default fitting conditions are used.");
324 }
325
326
ab886e13 327 //fit parameter arrays
acb9d358 328 Double_t fitParamTheta[100],fitParamPhi[100];
329 Double_t errFitParamTheta[100], errFitParamPhi[100];
330 Double_t binsXTheta[100],binsXPhi[100];
ab886e13 331
332 //copy histograms
ba06aaec 333 fHistH2InvPtTheta = (TH2F*)histThetaInvPt->Clone("invPtVsTheta");
334 fHistH2InvPtPhi = (TH2F*)histPhiInvPt->Clone("invPtVsPhi");
ab886e13 335
336 // initialize fit functions
acb9d358 337 InitFitFcn();
ab886e13 338
339 // canvas for 2D histograms (input)
340 TCanvas *thephiCan =new TCanvas("thephiCan","theta and phi vs invPt",800,400);
341 thephiCan->Divide(2,1);
342 thephiCan->cd(1);
343 fHistH2InvPtTheta->SetTitle("#theta vs charge/pt for selected #phi range");
344 fHistH2InvPtTheta->GetXaxis()->SetTitle("charge/p_{t} (GeV/c)^{-1}");
345 fHistH2InvPtTheta->GetYaxis()->SetTitle("#theta (rad)");
346 fHistH2InvPtTheta->DrawCopy("col1z");
347
348 thephiCan->cd(2);
349 fHistH2InvPtPhi->SetTitle("#phi vs charge/pt for selected #theta range");
350 fHistH2InvPtPhi->GetXaxis()->SetTitle("charge/p_{t} (GeV/c)^{-1}");
351 fHistH2InvPtPhi->GetYaxis()->SetTitle("#phi (rad)");
352 fHistH2InvPtPhi->DrawCopy("col1z");
353
354 // canvas for 1D histograms (projections)
acb9d358 355 TCanvas *theCan =new TCanvas("theCan","invPt theta bins",1200,900);
356 theCan->Divide((fNThetaBins+2)/2,2);
357 TCanvas *phiCan =new TCanvas("phiCan","invPt phi bins",1200,900);
358 phiCan->Divide((fNPhiBins+2)/2,2);
ab886e13 359
acb9d358 360 Int_t countPad = 1;
361
ab886e13 362 // analyse charge/pt in bins of theta
acb9d358 363 for(Int_t i=0;i<fNThetaBins;i++){
ab886e13 364
365 //set name for each projection
acb9d358 366 TString name = "fit_theta_";
367 name +=i;
ab886e13 368 fHistH2InvPtTheta->SetName(name.Data());
369
370 //find bins for projection
acb9d358 371 Int_t firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]);
372 if(i>0) firstBin +=1;
373 Int_t lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i+1]);
374 if( i == fNThetaBins-1) {
375 firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[0]);
376 lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]);
377 }
ab886e13 378
379 //make projection
380 fHistFitTheta[i] = (TH1D*)fHistH2InvPtTheta->ProjectionX("_px",firstBin,lastBin,"e");
acb9d358 381
382 Char_t titleTheta[50];
18e588e9 383 if(TMath::Abs(i-(fNThetaBins-1))<0.5) snprintf(titleTheta,50,"charge/pt (GeV/c) integrated over #theta");
384 else snprintf(titleTheta,50,"charge/pt (GeV/c) for #theta range: %1.3f - %1.3f",fThetaBins[i],fThetaBins[i+1]);
acb9d358 385
386 fHistFitTheta[i]->SetTitle(titleTheta);
387
388 Double_t invPtMinPos = 0;
389 Double_t invPtMinPosErr = 0;
390 Double_t invPtMinPosImpr = 0;
391 Double_t invPtMinPosErrImpr = 0;
392
6bcacaf1 393 if(fDoRebin) fHistFitTheta[i]->Rebin(fRebin);
ab886e13 394
ba06aaec 395 //start fitting
396 if(!fFitGaus){
ab886e13 397 Printf("making polynomial fit in charge/pt in theta bins");
acb9d358 398 theCan->cd(countPad);
399 MakeFit(fHistFitTheta[i],fFitMinPos, invPtMinPos,invPtMinPosErr, fExclRange,fRange);
400 MakeFitBetter(fHistFitTheta[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos,fExclRange,fRange);
401
ab886e13 402 //plot projection and fit
acb9d358 403 fHistFitTheta[i]->DrawCopy();
404 fFitMinPos->DrawCopy("L,same");
405 fFitMinPosRejP->DrawCopy("L,same");
406 }
407 else{
ab886e13 408 Printf("making gauss fit in charge/pt in theta bins");
acb9d358 409 theCan->cd(countPad);
410 MakeFitInvGauss(fHistFitTheta[i],fFitInvGauss, invPtMinPos, invPtMinPosErr,fExclRange,fRange);
411 MakeFitInvGaussBetter(fHistFitTheta[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange);
ab886e13 412
413 //plot projection and fit
acb9d358 414 fHistFitTheta[i]->DrawCopy();
415 fFitInvGauss->DrawCopy("L,same");
416 fFitInvGaussRejP->DrawCopy("L,same");
417 }
ab886e13 418 //add objects for analysis folder
acb9d358 419 aFolderObj->Add(fHistFitTheta[i]);
ab886e13 420
421 //store fit parameter
acb9d358 422 fitParamTheta[i] = invPtMinPosImpr;
423 errFitParamTheta[i] = invPtMinPosErrImpr;
ab886e13 424
425 //count bins and pad number for displaying
acb9d358 426 binsXTheta[i] = i+1.0;
427 countPad++;
428 }
429
430
431 countPad = 1;
432
433
ab886e13 434 // analyse charge/pt in bins of phi
acb9d358 435
436 for(Int_t i=0;i<fNPhiBins;i++){
ab886e13 437
438 //set name for each projection
acb9d358 439 TString name = "fit_phi_";
440 name +=i;
acb9d358 441 fHistH2InvPtPhi->SetName(name.Data());
ab886e13 442
443 //find bins for projection
acb9d358 444 Int_t firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]);
445 if(i>0) firstBin +=1;
446 Int_t lastBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i+1]);
b90a2249 447 if(TMath::Abs(i-(fNPhiBins-1))<0.5){
acb9d358 448 firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[0]);
449 lastBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]);
450 }
ab886e13 451
452 //make projection
453 fHistFitPhi[i] = (TH1D*) fHistH2InvPtPhi->ProjectionX("_px",firstBin,lastBin,"e");
acb9d358 454
455 Char_t titlePhi[50];
18e588e9 456 if(TMath::Abs(i-(fNPhiBins-1))<0.5) snprintf(titlePhi,50,"charge/pt (GeV/c) integrated over #phi");
457 else snprintf(titlePhi,50,"charge/pt (GeV/c) for #phi range: %1.3f - %1.3f",fPhiBins[i],fPhiBins[i+1]);
acb9d358 458
459 fHistFitPhi[i]->SetTitle(titlePhi);
460
461 Double_t invPtMinPos = 0;
462 Double_t invPtMinPosErr = 0;
463 Double_t invPtMinPosImpr = 0;
464 Double_t invPtMinPosErrImpr = 0;
ab886e13 465 if(fDoRebin) fHistFitPhi[i]->Rebin(fRebin);
466
467 //start fitting
ba06aaec 468 if(!fFitGaus){
ab886e13 469 Printf("making polynomial fit in charge/pt in phi bins");
acb9d358 470 phiCan->cd(countPad);
471 MakeFit(fHistFitPhi[i],fFitMinPos, invPtMinPos, invPtMinPosErr,fExclRange,fRange);
472 MakeFitBetter(fHistFitPhi[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange);
ab886e13 473
474 //plot projection and fit
acb9d358 475 fHistFitPhi[i]->DrawCopy();
476 fFitMinPos->DrawCopy("L,same");
477 fFitMinPosRejP->DrawCopy("L,same");
478
479 }
480 else {
ab886e13 481 Printf("making gauss fit in charge/pt in phi bins");
acb9d358 482 phiCan->cd(countPad);
483 MakeFitInvGauss(fHistFitPhi[i],fFitInvGauss, invPtMinPos, invPtMinPosErr, exclRange,fRange);
484 MakeFitInvGaussBetter(fHistFitPhi[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos, fExclRange,fRange);
ab886e13 485
486 //plot projection and fit
acb9d358 487 fHistFitPhi[i]->DrawCopy();
488 fFitInvGauss->DrawCopy("L,same");
489 fFitInvGaussRejP->DrawCopy("L,same");
490 }
ab886e13 491 //add objects for analysis folder
acb9d358 492 aFolderObj->Add(fHistFitPhi[i]);
ab886e13 493
494 //store fit parameter
acb9d358 495 fitParamPhi[i] = invPtMinPosImpr;
496 errFitParamPhi[i] = invPtMinPosErrImpr;
ab886e13 497
498 //count bins and pad number for displaying
acb9d358 499 binsXPhi[i] = i+1.0;
500 countPad++;
501 }
ab886e13 502
503 //initialize graphs for displaying the fit parameter
ba06aaec 504 InitGraphs(binsXTheta,fitParamTheta,errFitParamTheta,binsXPhi,fitParamPhi,errFitParamPhi);
acb9d358 505
ba06aaec 506 //plot fit values = minimum positions of charge/pt
acb9d358 507 TCanvas *canFitVal = new TCanvas("canFitVal","min pos histos",800,400);
508 canFitVal->Divide(2,1);
509
510 canFitVal->cd(1);
6bcacaf1 511 fGrMinPosTheta->Draw("AP");
acb9d358 512 canFitVal->cd(2);
6bcacaf1 513 fGrMinPosPhi->Draw("AP");
acb9d358 514
ab886e13 515 Printf("AliPerfAnalyzeInvPt::StartAnalysis: ***** NOTE: last bin is always fit result of integral over all angle BINS which have been set by user! angle ranges of 2D projections are kept.*****");
ba06aaec 516
517 //add objects to folder
acb9d358 518 aFolderObj->Add(fGrMinPosTheta);
519 aFolderObj->Add(fGrMinPosPhi);
520 aFolderObj->Add(fHistH2InvPtTheta);
521 aFolderObj->Add(fHistH2InvPtPhi);
522
523
524}
525
526
527//____________________________________________________________________________________________________________________________________________
528
ab886e13 529void AliPerfAnalyzeInvPt::MakeFit(TH1D *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range)
acb9d358 530{
ab886e13 531 // fit charge/pt and extract minimum position
ba06aaec 532
acb9d358 533 fitpb->SetRange(-range,range);
534 fitpb->SetParLimits(1,-0.05,0.05);
535 fitpb->SetParameter(0,1.0);
536 fitpb->FixParameter(4,excl);
537 fitpb->FixParameter(2,0.0);
538
ba06aaec 539 hproy->Fit(fitpb,"RM");
acb9d358 540 mean = fitpb->GetParameter(1);
541 errMean = fitpb->GetParError(1);
ba06aaec 542 if(!mean) errMean = 0.0;
acb9d358 543}
544//____________________________________________________________________________________________________________________________________________
ab886e13 545void AliPerfAnalyzeInvPt::MakeFitBetter(TH1D *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range)
acb9d358 546{
ba06aaec 547 // adjust fit range to minimum position of AliPerfAnalyzeInvPt::MakeFit and fit 1/pt and extract new minimum position
548
acb9d358 549 fitpb2->FixParameter(5,excl);
550 fitpb2->FixParameter(4,f);
551 // fitpb2->FixParameter(2,0.0);
ba06aaec 552 fitpb2->SetRange(-range+f,range+f);
acb9d358 553 fitpb2->SetParLimits(1,-0.05,0.05);
554 fitpb2->SetParameter(0,1.0);
ba06aaec 555 hproy->Fit(fitpb2,"RM");
acb9d358 556 mean = fitpb2->GetParameter(1);
557 errMean = fitpb2->GetParError(1);
ba06aaec 558 if(!mean) errMean = 0.0;
acb9d358 559
560}
561//____________________________________________________________________________________________________________________________________________
ab886e13 562void AliPerfAnalyzeInvPt::MakeFitInvGauss(TH1D *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range)
acb9d358 563{
ab886e13 564 // fit charge/pt and extract minimum position
ba06aaec 565
acb9d358 566 fitpb->FixParameter(6,excl);
567 fitpb->SetRange(-range,range);
568 fitpb->SetParameter(0,-1.0);
569 fitpb->SetParameter(2,1.0);
570 fitpb->SetParameter(3,25000.0);
571 fitpb->SetParameter(4,-1.0);
572 fitpb->SetParLimits(1,-0.02,0.02);
ba06aaec 573 hproy->Fit(fitpb,"RM");
acb9d358 574 mean = fitpb->GetParameter(1);
575 errMean = fitpb->GetParError(1);
ba06aaec 576 if(!mean) errMean = 0.0;
acb9d358 577
578}
579//____________________________________________________________________________________________________________________________________________
ab886e13 580void AliPerfAnalyzeInvPt::MakeFitInvGaussBetter(TH1D *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range)
acb9d358 581{
ba06aaec 582 // adjust fit range to minimum position of AliPerfAnalyzeInvPt::MakeFitInvGauss and fit 1/pt and extract new minimum position
583
acb9d358 584 fitpb2->FixParameter(7,excl);
585 fitpb2->FixParameter(6,f);
586 fitpb2->SetRange(-range+f,range+f);
587 fitpb2->SetParameter(0,-1.0);
588 fitpb2->SetParameter(2,1.0);
589 fitpb2->SetParameter(3,25000.0);
590 fitpb2->SetParameter(4,-1.0);
591 fitpb2->SetParLimits(1,-0.02,0.02);
ba06aaec 592 hproy->Fit(fitpb2,"RM");
acb9d358 593 mean = fitpb2->GetParameter(1);
594 errMean = fitpb2->GetParError(1);
ba06aaec 595 if(!mean) errMean = 0.0;
acb9d358 596
597}
598
599
600//____________________________________________________________________________________________________________________________________________
acb9d358 601
ba06aaec 602// set variables for analysis
603void AliPerfAnalyzeInvPt::SetProjBinsPhi(const Double_t *phiBinArray,const Int_t nphBins){
acb9d358 604
ba06aaec 605 // set theta bins for Analyse()
606 //set phi bins as array and set number of this array which is equal to number of bins analysed
607 //the last analysed bin will always be the projection from first to last bin in the array
608 if(nphBins){
609 fNPhiBins = nphBins;
acb9d358 610
ba06aaec 611 for(Int_t k = 0;k<fNPhiBins;k++){
612 fPhiBins[k] = phiBinArray[k];
613 }
614 Printf("AliPerfAnalyzeInvPt::SetProjBinsPhi:number of bins in phi set to %i",fNPhiBins);
acb9d358 615 }
ba06aaec 616 else Printf("AliPerfAnalyzeInvPt::SetProjBinsPhi: number of bins in theta NOT set. Default values are used.");
acb9d358 617}
618//____________________________________________________________________________________________________________________________________________
ba06aaec 619void AliPerfAnalyzeInvPt::SetProjBinsTheta(const Double_t *thetaBinArray, const Int_t nthBins){
620 // set theta bins for Analyse()
621 //set theta bins as array and set number of this array which is equal to number of bins analysed
622 //the last analysed bin will always be the projection from first to last bin in the array
623 if(nthBins){
624 fNThetaBins = nthBins;
625 for(Int_t k = 0;k<fNThetaBins;k++){
626 fThetaBins[k] = thetaBinArray[k];
627 }
628 Printf("AliPerfAnalyzeInvPt::SetProjBinsTheta: number of bins in theta set to %i",fNThetaBins);
acb9d358 629 }
ba06aaec 630 else Printf("AliPerfAnalyzeInvPt::SetProjBinsTheta: number of bins in theta NOT set. Default values are used.");
acb9d358 631}
632//____________________________________________________________________________________________________________________________________________
633void AliPerfAnalyzeInvPt::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
634
ba06aaec 635 //set the fit options:
636 //for usage of gaussian function instead of polynomial (default) set setGausFit=kTRUE
637 //set the range of rejection of points around 0 via exclusionR
638 //set the fit range around 0 with fitR
639 if(fitR){
640 fFitGaus = setGausFit;
641 fExclRange = exclusionR;
642 fRange = fitR;
acb9d358 643
ab886e13 644 if(fFitGaus) Printf("AliPerfAnalyzeInvPt::SetMakeFitOption: set MakeGausFit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange);
645 else Printf("AliPerfAnalyzeInvPt::SetMakeFitOption: set standard polynomial fit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange);
ba06aaec 646 }
647 else Printf(" AliPerfAnalyzeInvPt::SetMakeFitOption: no user input. Set standard polynomial fit with fit range 1.0 and exclusion range 0.13.");
acb9d358 648}