]>
Commit | Line | Data |
---|---|---|
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 | |
38 | TFile f("Output.root"); | |
ab886e13 | 39 | AliPerformancePtCalib *compObj = (AliPerformancePtCalib*)coutput->FindObject("AliPerformancePtCalib"); |
40 | //or | |
41 | AliPerformancePtCalibMC *compObj = (AliPerformancePtCalibMC*)coutput->FindObject("AliPerformancePtCalibMC"); | |
42 | ||
43 | // set phi and theta bins for fitting and analyse comparison data | |
44 | compObj->SetProjBinsTheta(thetaBins,nThetaBins,minPhi, maxPhi); | |
45 | compObj->SetProjBinsPhi(phiBins,nPhiBins,minTheta,maxTheta); | |
46 | compObj->SetMakeFitOption(kFALSE,exclRange,fitRange); | |
47 | compObj->SetDoRebin(rebin); | |
ba06aaec | 48 | compObj->Analyse(); |
ab886e13 | 49 | //for details see functions of class AliPerformancePtCalib* |
ba06aaec | 50 | |
51 | // the output histograms/graphs will be stored in the folder "folderRes" | |
52 | compObj->GetAnalysisFolder()->ls("*"); | |
53 | ||
54 | // user can save whole comparison object (or only folder with anlysed histograms) | |
55 | // in the seperate output file (e.g.) | |
56 | TFile fout("Analysed_InvPt.root","recreate"); | |
57 | compObj->Write(); // compObj->GetAnalysisFolder()->Write(); | |
58 | fout.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 | ||
74 | ClassImp(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 | 93 | Double_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 | 111 | Double_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 | 122 | Double_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 | //_____________________________________________________________________________________________________________________________________________ | |
4823e4d4 | 141 | AliPerfAnalyzeInvPt::AliPerfAnalyzeInvPt(const Char_t* name, const Char_t* title): |
142 | TNamed(name,title), | |
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 | } |
acb9d358 | 182 | |
ba06aaec | 183 | //______________________________________________________________________________________________________________________________________ |
184 | void AliPerfAnalyzeInvPt::InitGraphs(Double_t *binsXTheta,Double_t *fitParamTheta,Double_t *errFitParamTheta,Double_t *binsXPhi,Double_t *fitParamPhi,Double_t *errFitParamPhi){ | |
acb9d358 | 185 | |
ba06aaec | 186 | // initialize graphs for storing fit parameter |
acb9d358 | 187 | |
ba06aaec | 188 | fGrMinPosTheta = new TGraphErrors(fNThetaBins,binsXTheta,fitParamTheta,0, errFitParamTheta); |
189 | fGrMinPosTheta->SetMarkerStyle(20); | |
190 | fGrMinPosTheta->SetMarkerColor(2); | |
191 | fGrMinPosTheta->SetLineColor(2); | |
ab886e13 | 192 | fGrMinPosTheta->GetYaxis()->SetTitle("min pos (GeV/c)^{-1}"); |
ba06aaec | 193 | fGrMinPosTheta->GetXaxis()->SetTitle("#theta bin no."); |
194 | fGrMinPosTheta->GetYaxis()->SetTitleOffset(1.2); | |
ab886e13 | 195 | fGrMinPosTheta->SetTitle("#theta bins"); |
ba06aaec | 196 | |
197 | fGrMinPosPhi = new TGraphErrors(fNPhiBins,binsXPhi,fitParamPhi,0,errFitParamPhi); | |
198 | fGrMinPosPhi->SetMarkerStyle(20); | |
199 | fGrMinPosPhi->SetMarkerColor(4); | |
200 | fGrMinPosPhi->SetLineColor(4); | |
ab886e13 | 201 | fGrMinPosPhi->GetYaxis()->SetTitle("min pos (GeV/c)^{-1}"); |
ba06aaec | 202 | fGrMinPosPhi->GetXaxis()->SetTitle("#phi bin no."); |
203 | fGrMinPosPhi->GetYaxis()->SetTitleOffset(1.2); | |
ab886e13 | 204 | fGrMinPosPhi->SetTitle("#phi bins"); |
acb9d358 | 205 | } |
206 | ||
207 | //______________________________________________________________________________________________________________________________________ | |
208 | ||
209 | void AliPerfAnalyzeInvPt::InitFitFcn(){ | |
ba06aaec | 210 | // fit functions |
211 | fFitMinPos = new TF1("fFitMinPos", Polynomial,-4.0,4.0,5); | |
212 | fFitMinPos->SetLineColor(4); | |
213 | fFitMinPos->SetLineWidth(1); | |
214 | fFitMinPos->SetParameter(0,1.0); | |
215 | fFitMinPos->SetParameter(1,0.0); | |
216 | fFitMinPos->SetParameter(2,0.0); | |
217 | ||
218 | fFitMinPosRejP = new TF1("fFitMinPosRejP",PolynomialRejP,-4.0,4.0,6); | |
219 | fFitMinPosRejP->SetLineColor(2); | |
220 | fFitMinPosRejP->SetLineWidth(1); | |
221 | fFitMinPosRejP->SetParameter(0,1.0); | |
222 | fFitMinPosRejP->SetParameter(1,0.0); | |
223 | fFitMinPosRejP->SetParameter(2,0.0); | |
acb9d358 | 224 | |
225 | ||
ba06aaec | 226 | fFitInvGauss = new TF1("fFitInvGauss", InvGauss,-4.0,4.0,6); |
227 | fFitInvGauss->SetLineColor(4); | |
228 | fFitInvGauss->SetLineWidth(1); | |
229 | fFitInvGauss->SetParameter(2,1.0); | |
acb9d358 | 230 | |
ba06aaec | 231 | fFitInvGaussRejP = new TF1("fFitInvGaussRejP", InvGaussRejP,-4.0,4.0,7); |
232 | fFitInvGaussRejP->SetLineColor(2); | |
233 | fFitInvGaussRejP->SetLineWidth(1); | |
234 | fFitInvGaussRejP->SetParameter(2,1.0); | |
acb9d358 | 235 | |
ba06aaec | 236 | } |
acb9d358 | 237 | //______________________________________________________________________________________________________________________________________ |
ba06aaec | 238 | void AliPerfAnalyzeInvPt::StartAnalysis(const TH2F *histThetaInvPt, const TH2F *histPhiInvPt, TObjArray* aFolderObj){ |
ab886e13 | 239 | //start analysis: fitting charge/pt spectra |
acb9d358 | 240 | |
241 | ||
242 | ||
243 | if(!histThetaInvPt) { | |
244 | Printf("warning: no 1/pt histogram to analyse in theta bins!"); | |
18e588e9 | 245 | return; |
acb9d358 | 246 | } |
247 | if(!histPhiInvPt) { | |
248 | Printf("warning: no 1/pt histogram to analyse in phit bins!"); | |
18e588e9 | 249 | return; |
acb9d358 | 250 | } |
251 | ||
ab886e13 | 252 | 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 |
253 | 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 | 254 | |
255 | Int_t nThetaBins = 9; | |
256 | Int_t nPhiBins = 13; | |
257 | Double_t range = 1.0; // fit range | |
258 | Double_t exclRange = 0.13; // range of point rejection in fit | |
259 | Bool_t fitGaus=kFALSE; // fit with Gaussian or with polynomial (kFALSE) | |
260 | ||
ba06aaec | 261 | if(!fNThetaBins){ |
acb9d358 | 262 | fNThetaBins = nThetaBins; |
263 | for(Int_t k = 0;k<fNThetaBins;k++){ | |
264 | fThetaBins[k] = thetaBins[k]; | |
265 | } | |
266 | Printf("warning: for theta bins no user intput! bins are set to default values."); | |
267 | } | |
268 | ||
ba06aaec | 269 | if(!fNPhiBins){ |
acb9d358 | 270 | fNPhiBins = nPhiBins; |
271 | for(Int_t k = 0;k<fNPhiBins;k++){ | |
272 | fPhiBins[k] = phiBins[k]; | |
273 | } | |
274 | ||
275 | Printf("warning: for phi bins no user intput! bins are set to default values."); | |
276 | ||
277 | } | |
278 | ||
ba06aaec | 279 | if(!fRange){ |
acb9d358 | 280 | fRange = range; |
281 | fExclRange = exclRange; | |
282 | fFitGaus = fitGaus; | |
283 | Printf("warning: no fit range is set. default fitting conditions are used."); | |
284 | } | |
285 | ||
286 | ||
ab886e13 | 287 | //fit parameter arrays |
acb9d358 | 288 | Double_t fitParamTheta[100],fitParamPhi[100]; |
289 | Double_t errFitParamTheta[100], errFitParamPhi[100]; | |
290 | Double_t binsXTheta[100],binsXPhi[100]; | |
ab886e13 | 291 | |
292 | //copy histograms | |
ba06aaec | 293 | fHistH2InvPtTheta = (TH2F*)histThetaInvPt->Clone("invPtVsTheta"); |
294 | fHistH2InvPtPhi = (TH2F*)histPhiInvPt->Clone("invPtVsPhi"); | |
ab886e13 | 295 | |
296 | // initialize fit functions | |
acb9d358 | 297 | InitFitFcn(); |
ab886e13 | 298 | |
299 | // canvas for 2D histograms (input) | |
300 | TCanvas *thephiCan =new TCanvas("thephiCan","theta and phi vs invPt",800,400); | |
301 | thephiCan->Divide(2,1); | |
302 | thephiCan->cd(1); | |
303 | fHistH2InvPtTheta->SetTitle("#theta vs charge/pt for selected #phi range"); | |
304 | fHistH2InvPtTheta->GetXaxis()->SetTitle("charge/p_{t} (GeV/c)^{-1}"); | |
305 | fHistH2InvPtTheta->GetYaxis()->SetTitle("#theta (rad)"); | |
306 | fHistH2InvPtTheta->DrawCopy("col1z"); | |
307 | ||
308 | thephiCan->cd(2); | |
309 | fHistH2InvPtPhi->SetTitle("#phi vs charge/pt for selected #theta range"); | |
310 | fHistH2InvPtPhi->GetXaxis()->SetTitle("charge/p_{t} (GeV/c)^{-1}"); | |
311 | fHistH2InvPtPhi->GetYaxis()->SetTitle("#phi (rad)"); | |
312 | fHistH2InvPtPhi->DrawCopy("col1z"); | |
313 | ||
314 | // canvas for 1D histograms (projections) | |
acb9d358 | 315 | TCanvas *theCan =new TCanvas("theCan","invPt theta bins",1200,900); |
316 | theCan->Divide((fNThetaBins+2)/2,2); | |
317 | TCanvas *phiCan =new TCanvas("phiCan","invPt phi bins",1200,900); | |
318 | phiCan->Divide((fNPhiBins+2)/2,2); | |
ab886e13 | 319 | |
acb9d358 | 320 | Int_t countPad = 1; |
321 | ||
ab886e13 | 322 | // analyse charge/pt in bins of theta |
acb9d358 | 323 | for(Int_t i=0;i<fNThetaBins;i++){ |
ab886e13 | 324 | |
325 | //set name for each projection | |
acb9d358 | 326 | TString name = "fit_theta_"; |
327 | name +=i; | |
ab886e13 | 328 | fHistH2InvPtTheta->SetName(name.Data()); |
329 | ||
330 | //find bins for projection | |
acb9d358 | 331 | Int_t firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]); |
332 | if(i>0) firstBin +=1; | |
333 | Int_t lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i+1]); | |
334 | if( i == fNThetaBins-1) { | |
335 | firstBin= fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[0]); | |
336 | lastBin = fHistH2InvPtTheta->GetYaxis()->FindBin(fThetaBins[i]); | |
337 | } | |
ab886e13 | 338 | |
339 | //make projection | |
340 | fHistFitTheta[i] = (TH1D*)fHistH2InvPtTheta->ProjectionX("_px",firstBin,lastBin,"e"); | |
acb9d358 | 341 | |
342 | Char_t titleTheta[50]; | |
18e588e9 | 343 | if(TMath::Abs(i-(fNThetaBins-1))<0.5) snprintf(titleTheta,50,"charge/pt (GeV/c) integrated over #theta"); |
344 | else snprintf(titleTheta,50,"charge/pt (GeV/c) for #theta range: %1.3f - %1.3f",fThetaBins[i],fThetaBins[i+1]); | |
acb9d358 | 345 | |
346 | fHistFitTheta[i]->SetTitle(titleTheta); | |
347 | ||
348 | Double_t invPtMinPos = 0; | |
349 | Double_t invPtMinPosErr = 0; | |
350 | Double_t invPtMinPosImpr = 0; | |
351 | Double_t invPtMinPosErrImpr = 0; | |
352 | ||
6bcacaf1 | 353 | if(fDoRebin) fHistFitTheta[i]->Rebin(fRebin); |
ab886e13 | 354 | |
ba06aaec | 355 | //start fitting |
356 | if(!fFitGaus){ | |
ab886e13 | 357 | Printf("making polynomial fit in charge/pt in theta bins"); |
acb9d358 | 358 | theCan->cd(countPad); |
359 | MakeFit(fHistFitTheta[i],fFitMinPos, invPtMinPos,invPtMinPosErr, fExclRange,fRange); | |
360 | MakeFitBetter(fHistFitTheta[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos,fExclRange,fRange); | |
361 | ||
ab886e13 | 362 | //plot projection and fit |
acb9d358 | 363 | fHistFitTheta[i]->DrawCopy(); |
364 | fFitMinPos->DrawCopy("L,same"); | |
365 | fFitMinPosRejP->DrawCopy("L,same"); | |
366 | } | |
367 | else{ | |
ab886e13 | 368 | Printf("making gauss fit in charge/pt in theta bins"); |
acb9d358 | 369 | theCan->cd(countPad); |
370 | MakeFitInvGauss(fHistFitTheta[i],fFitInvGauss, invPtMinPos, invPtMinPosErr,fExclRange,fRange); | |
371 | MakeFitInvGaussBetter(fHistFitTheta[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange); | |
ab886e13 | 372 | |
373 | //plot projection and fit | |
acb9d358 | 374 | fHistFitTheta[i]->DrawCopy(); |
375 | fFitInvGauss->DrawCopy("L,same"); | |
376 | fFitInvGaussRejP->DrawCopy("L,same"); | |
377 | } | |
ab886e13 | 378 | //add objects for analysis folder |
acb9d358 | 379 | aFolderObj->Add(fHistFitTheta[i]); |
ab886e13 | 380 | |
381 | //store fit parameter | |
acb9d358 | 382 | fitParamTheta[i] = invPtMinPosImpr; |
383 | errFitParamTheta[i] = invPtMinPosErrImpr; | |
ab886e13 | 384 | |
385 | //count bins and pad number for displaying | |
acb9d358 | 386 | binsXTheta[i] = i+1.0; |
387 | countPad++; | |
388 | } | |
389 | ||
390 | ||
391 | countPad = 1; | |
392 | ||
393 | ||
ab886e13 | 394 | // analyse charge/pt in bins of phi |
acb9d358 | 395 | |
396 | for(Int_t i=0;i<fNPhiBins;i++){ | |
ab886e13 | 397 | |
398 | //set name for each projection | |
acb9d358 | 399 | TString name = "fit_phi_"; |
400 | name +=i; | |
acb9d358 | 401 | fHistH2InvPtPhi->SetName(name.Data()); |
ab886e13 | 402 | |
403 | //find bins for projection | |
acb9d358 | 404 | Int_t firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]); |
405 | if(i>0) firstBin +=1; | |
406 | Int_t lastBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i+1]); | |
b90a2249 | 407 | if(TMath::Abs(i-(fNPhiBins-1))<0.5){ |
acb9d358 | 408 | firstBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[0]); |
409 | lastBin = fHistH2InvPtPhi->GetYaxis()->FindBin(fPhiBins[i]); | |
410 | } | |
ab886e13 | 411 | |
412 | //make projection | |
413 | fHistFitPhi[i] = (TH1D*) fHistH2InvPtPhi->ProjectionX("_px",firstBin,lastBin,"e"); | |
acb9d358 | 414 | |
415 | Char_t titlePhi[50]; | |
18e588e9 | 416 | if(TMath::Abs(i-(fNPhiBins-1))<0.5) snprintf(titlePhi,50,"charge/pt (GeV/c) integrated over #phi"); |
417 | else snprintf(titlePhi,50,"charge/pt (GeV/c) for #phi range: %1.3f - %1.3f",fPhiBins[i],fPhiBins[i+1]); | |
acb9d358 | 418 | |
419 | fHistFitPhi[i]->SetTitle(titlePhi); | |
420 | ||
421 | Double_t invPtMinPos = 0; | |
422 | Double_t invPtMinPosErr = 0; | |
423 | Double_t invPtMinPosImpr = 0; | |
424 | Double_t invPtMinPosErrImpr = 0; | |
ab886e13 | 425 | if(fDoRebin) fHistFitPhi[i]->Rebin(fRebin); |
426 | ||
427 | //start fitting | |
ba06aaec | 428 | if(!fFitGaus){ |
ab886e13 | 429 | Printf("making polynomial fit in charge/pt in phi bins"); |
acb9d358 | 430 | phiCan->cd(countPad); |
431 | MakeFit(fHistFitPhi[i],fFitMinPos, invPtMinPos, invPtMinPosErr,fExclRange,fRange); | |
432 | MakeFitBetter(fHistFitPhi[i],fFitMinPosRejP, invPtMinPosImpr, invPtMinPosErrImpr,invPtMinPos,fExclRange,fRange); | |
ab886e13 | 433 | |
434 | //plot projection and fit | |
acb9d358 | 435 | fHistFitPhi[i]->DrawCopy(); |
436 | fFitMinPos->DrawCopy("L,same"); | |
437 | fFitMinPosRejP->DrawCopy("L,same"); | |
438 | ||
439 | } | |
440 | else { | |
ab886e13 | 441 | Printf("making gauss fit in charge/pt in phi bins"); |
acb9d358 | 442 | phiCan->cd(countPad); |
443 | MakeFitInvGauss(fHistFitPhi[i],fFitInvGauss, invPtMinPos, invPtMinPosErr, exclRange,fRange); | |
444 | MakeFitInvGaussBetter(fHistFitPhi[i],fFitInvGaussRejP, invPtMinPosImpr, invPtMinPosErrImpr, invPtMinPos, fExclRange,fRange); | |
ab886e13 | 445 | |
446 | //plot projection and fit | |
acb9d358 | 447 | fHistFitPhi[i]->DrawCopy(); |
448 | fFitInvGauss->DrawCopy("L,same"); | |
449 | fFitInvGaussRejP->DrawCopy("L,same"); | |
450 | } | |
ab886e13 | 451 | //add objects for analysis folder |
acb9d358 | 452 | aFolderObj->Add(fHistFitPhi[i]); |
ab886e13 | 453 | |
454 | //store fit parameter | |
acb9d358 | 455 | fitParamPhi[i] = invPtMinPosImpr; |
456 | errFitParamPhi[i] = invPtMinPosErrImpr; | |
ab886e13 | 457 | |
458 | //count bins and pad number for displaying | |
acb9d358 | 459 | binsXPhi[i] = i+1.0; |
460 | countPad++; | |
461 | } | |
ab886e13 | 462 | |
463 | //initialize graphs for displaying the fit parameter | |
ba06aaec | 464 | InitGraphs(binsXTheta,fitParamTheta,errFitParamTheta,binsXPhi,fitParamPhi,errFitParamPhi); |
acb9d358 | 465 | |
ba06aaec | 466 | //plot fit values = minimum positions of charge/pt |
acb9d358 | 467 | TCanvas *canFitVal = new TCanvas("canFitVal","min pos histos",800,400); |
468 | canFitVal->Divide(2,1); | |
469 | ||
470 | canFitVal->cd(1); | |
6bcacaf1 | 471 | fGrMinPosTheta->Draw("AP"); |
acb9d358 | 472 | canFitVal->cd(2); |
6bcacaf1 | 473 | fGrMinPosPhi->Draw("AP"); |
acb9d358 | 474 | |
ab886e13 | 475 | 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 | 476 | |
477 | //add objects to folder | |
acb9d358 | 478 | aFolderObj->Add(fGrMinPosTheta); |
479 | aFolderObj->Add(fGrMinPosPhi); | |
480 | aFolderObj->Add(fHistH2InvPtTheta); | |
481 | aFolderObj->Add(fHistH2InvPtPhi); | |
482 | ||
483 | ||
484 | } | |
485 | ||
486 | ||
487 | //____________________________________________________________________________________________________________________________________________ | |
488 | ||
ab886e13 | 489 | void AliPerfAnalyzeInvPt::MakeFit(TH1D *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range) |
acb9d358 | 490 | { |
ab886e13 | 491 | // fit charge/pt and extract minimum position |
ba06aaec | 492 | |
acb9d358 | 493 | fitpb->SetRange(-range,range); |
494 | fitpb->SetParLimits(1,-0.05,0.05); | |
495 | fitpb->SetParameter(0,1.0); | |
496 | fitpb->FixParameter(4,excl); | |
497 | fitpb->FixParameter(2,0.0); | |
498 | ||
ba06aaec | 499 | hproy->Fit(fitpb,"RM"); |
acb9d358 | 500 | mean = fitpb->GetParameter(1); |
501 | errMean = fitpb->GetParError(1); | |
ba06aaec | 502 | if(!mean) errMean = 0.0; |
acb9d358 | 503 | } |
504 | //____________________________________________________________________________________________________________________________________________ | |
ab886e13 | 505 | void AliPerfAnalyzeInvPt::MakeFitBetter(TH1D *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range) |
acb9d358 | 506 | { |
ba06aaec | 507 | // adjust fit range to minimum position of AliPerfAnalyzeInvPt::MakeFit and fit 1/pt and extract new minimum position |
508 | ||
acb9d358 | 509 | fitpb2->FixParameter(5,excl); |
510 | fitpb2->FixParameter(4,f); | |
511 | // fitpb2->FixParameter(2,0.0); | |
ba06aaec | 512 | fitpb2->SetRange(-range+f,range+f); |
acb9d358 | 513 | fitpb2->SetParLimits(1,-0.05,0.05); |
514 | fitpb2->SetParameter(0,1.0); | |
ba06aaec | 515 | hproy->Fit(fitpb2,"RM"); |
acb9d358 | 516 | mean = fitpb2->GetParameter(1); |
517 | errMean = fitpb2->GetParError(1); | |
ba06aaec | 518 | if(!mean) errMean = 0.0; |
acb9d358 | 519 | |
520 | } | |
521 | //____________________________________________________________________________________________________________________________________________ | |
ab886e13 | 522 | void AliPerfAnalyzeInvPt::MakeFitInvGauss(TH1D *hproy, TF1 * fitpb, Double_t &mean, Double_t &errMean, Double_t &excl,Double_t &range) |
acb9d358 | 523 | { |
ab886e13 | 524 | // fit charge/pt and extract minimum position |
ba06aaec | 525 | |
acb9d358 | 526 | fitpb->FixParameter(6,excl); |
527 | fitpb->SetRange(-range,range); | |
528 | fitpb->SetParameter(0,-1.0); | |
529 | fitpb->SetParameter(2,1.0); | |
530 | fitpb->SetParameter(3,25000.0); | |
531 | fitpb->SetParameter(4,-1.0); | |
532 | fitpb->SetParLimits(1,-0.02,0.02); | |
ba06aaec | 533 | hproy->Fit(fitpb,"RM"); |
acb9d358 | 534 | mean = fitpb->GetParameter(1); |
535 | errMean = fitpb->GetParError(1); | |
ba06aaec | 536 | if(!mean) errMean = 0.0; |
acb9d358 | 537 | |
538 | } | |
539 | //____________________________________________________________________________________________________________________________________________ | |
ab886e13 | 540 | void AliPerfAnalyzeInvPt::MakeFitInvGaussBetter(TH1D *hproy, TF1 * fitpb2, Double_t &mean, Double_t &errMean, Double_t &f, Double_t &excl, Double_t &range) |
acb9d358 | 541 | { |
ba06aaec | 542 | // adjust fit range to minimum position of AliPerfAnalyzeInvPt::MakeFitInvGauss and fit 1/pt and extract new minimum position |
543 | ||
acb9d358 | 544 | fitpb2->FixParameter(7,excl); |
545 | fitpb2->FixParameter(6,f); | |
546 | fitpb2->SetRange(-range+f,range+f); | |
547 | fitpb2->SetParameter(0,-1.0); | |
548 | fitpb2->SetParameter(2,1.0); | |
549 | fitpb2->SetParameter(3,25000.0); | |
550 | fitpb2->SetParameter(4,-1.0); | |
551 | fitpb2->SetParLimits(1,-0.02,0.02); | |
ba06aaec | 552 | hproy->Fit(fitpb2,"RM"); |
acb9d358 | 553 | mean = fitpb2->GetParameter(1); |
554 | errMean = fitpb2->GetParError(1); | |
ba06aaec | 555 | if(!mean) errMean = 0.0; |
acb9d358 | 556 | |
557 | } | |
558 | ||
559 | ||
560 | //____________________________________________________________________________________________________________________________________________ | |
acb9d358 | 561 | |
ba06aaec | 562 | // set variables for analysis |
563 | void AliPerfAnalyzeInvPt::SetProjBinsPhi(const Double_t *phiBinArray,const Int_t nphBins){ | |
acb9d358 | 564 | |
ba06aaec | 565 | // set theta bins for Analyse() |
566 | //set phi bins as array and set number of this array which is equal to number of bins analysed | |
567 | //the last analysed bin will always be the projection from first to last bin in the array | |
568 | if(nphBins){ | |
569 | fNPhiBins = nphBins; | |
acb9d358 | 570 | |
ba06aaec | 571 | for(Int_t k = 0;k<fNPhiBins;k++){ |
572 | fPhiBins[k] = phiBinArray[k]; | |
573 | } | |
574 | Printf("AliPerfAnalyzeInvPt::SetProjBinsPhi:number of bins in phi set to %i",fNPhiBins); | |
acb9d358 | 575 | } |
ba06aaec | 576 | else Printf("AliPerfAnalyzeInvPt::SetProjBinsPhi: number of bins in theta NOT set. Default values are used."); |
acb9d358 | 577 | } |
578 | //____________________________________________________________________________________________________________________________________________ | |
ba06aaec | 579 | void AliPerfAnalyzeInvPt::SetProjBinsTheta(const Double_t *thetaBinArray, const Int_t nthBins){ |
580 | // set theta bins for Analyse() | |
581 | //set theta bins as array and set number of this array which is equal to number of bins analysed | |
582 | //the last analysed bin will always be the projection from first to last bin in the array | |
583 | if(nthBins){ | |
584 | fNThetaBins = nthBins; | |
585 | for(Int_t k = 0;k<fNThetaBins;k++){ | |
586 | fThetaBins[k] = thetaBinArray[k]; | |
587 | } | |
588 | Printf("AliPerfAnalyzeInvPt::SetProjBinsTheta: number of bins in theta set to %i",fNThetaBins); | |
acb9d358 | 589 | } |
ba06aaec | 590 | else Printf("AliPerfAnalyzeInvPt::SetProjBinsTheta: number of bins in theta NOT set. Default values are used."); |
acb9d358 | 591 | } |
592 | //____________________________________________________________________________________________________________________________________________ | |
593 | void AliPerfAnalyzeInvPt::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){ | |
594 | ||
ba06aaec | 595 | //set the fit options: |
596 | //for usage of gaussian function instead of polynomial (default) set setGausFit=kTRUE | |
597 | //set the range of rejection of points around 0 via exclusionR | |
598 | //set the fit range around 0 with fitR | |
599 | if(fitR){ | |
600 | fFitGaus = setGausFit; | |
601 | fExclRange = exclusionR; | |
602 | fRange = fitR; | |
acb9d358 | 603 | |
ab886e13 | 604 | if(fFitGaus) Printf("AliPerfAnalyzeInvPt::SetMakeFitOption: set MakeGausFit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange); |
605 | 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 | 606 | } |
607 | else Printf(" AliPerfAnalyzeInvPt::SetMakeFitOption: no user input. Set standard polynomial fit with fit range 1.0 and exclusion range 0.13."); | |
acb9d358 | 608 | } |