Belle II Software development
XT.h
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include "TProfile.h"
9#include "TCanvas.h"
10#include "TF1.h"
11#include "TH1D.h"
12#include "Math/ChebyshevPol.h"
13#include "iostream"
14
18// cppcheck-suppress constParameter
19Double_t pol5pol1(Double_t* x, Double_t* par)
20{
21 Double_t xx = x[0];
22 Double_t x6, x2;
23 Double_t f, ctp;
24 if (xx < par[6]) {
25 f = par[0] + par[1] * xx + par[2] * xx * xx + par[3] * xx * xx * xx + par[4] * xx * xx * xx * xx + par[5] * xx * xx * xx * xx * xx;
26 } else {
27 x6 = par[6];
28 x2 = xx - x6;
29 ctp = par[0] + par[1] * x6 + par[2] * x6 * x6 + par[3] * x6 * x6 * x6 + par[4] * x6 * x6 * x6 * x6 + par[5] * x6 * x6 * x6 * x6 *
30 x6;
31 f = par[7] * x2 + ctp;
32 }
33 return f;
34}
35
39// cppcheck-suppress constParameter
40Double_t Cheb5pol1(Double_t* x, Double_t* par)
41{
42 Double_t xx = x[0];
43 Double_t x6, x2;
44 Double_t f, ctp;
45
46 if (xx < par[6]) {
47 f = ROOT::Math::Chebyshev5(xx, par[0], par[1], par[2], par[3], par[4], par[5]);
48 } else {
49 x6 = par[6];
50 x2 = xx - x6;
51 ctp = ROOT::Math::Chebyshev5(x6, par[0], par[1], par[2], par[3], par[4], par[5]);
52 f = par[7] * x2 + ctp;
53 }
54 return f;
55}
60class XT {
61public:
65 explicit XT(TProfile* h1)
66 {
67 m_h1 = (TProfile*)h1->Clone();
68 m_h1->SetDirectory(0);
69 }
73 XT(TProfile* h1, int mode)
74 {
75 m_h1 = (TProfile*)h1->Clone();
76 m_h1->SetDirectory(0);
77 m_mode = mode;
78 }
79
83 XT(TH1D* h1, int mode)
84 {
85 m_h1 = (TProfile*)h1->Clone();
86 m_h1->SetDirectory(0);
87 m_mode = mode;
88 }
89
93 void setP6(double p6)
94 {
95 m_XTParam[6] = p6;
96 }
97
103 void setMode(int mode)
104 {
105 m_mode = mode;
106 }
110 void BField(bool bfield) {m_BField = bfield;}
111
115 void setXTParams(const double p[8])
116 {
117 for (int i = 0; i < 8; ++i) {m_XTParam[i] = p[i];}
118 m_tmax = p[6] + 250;
119 }
123 void setXTParams(double p0, double p1, double p2, double p3,
124 double p4, double p5, double p6, double p7)
125 {
126 m_XTParam[0] = p0; m_XTParam[1] = p1; m_XTParam[2] = p2;
127 m_XTParam[3] = p3; m_XTParam[4] = p4; m_XTParam[5] = p5;
128 m_XTParam[6] = p6; m_XTParam[7] = p7;
129 m_tmax = p6 + 250;
130 }
131
135 void setFitRange(double tmin, double tmax)
136 {
137 m_tmin = tmin;
138 m_tmax = tmax;
139 }
144 {
145 m_minRequiredEntry = min;
146 }
150 void setDebug(bool debug)
151 {
152 m_debug = debug;
153 }
158 {
159 return m_fitflag;
160 }
161
165 double getProb()
166 {
167 return m_Prob;
168 }
172 void getFittedXTParams(double pa[8])
173 {
174 for (int i = 0; i < 8; ++i) {pa[i] = m_FittedXTParams[i];}
175 }
179 TF1* getXTFunction(int mode)
180 {
181 if (mode == 0) return xtpol5;
182 else return xtCheb5;
183 }
188 {
189 if (m_mode == 0) return xtpol5;
190 else return xtCheb5;
191 };
195 TProfile* getFittedHisto() {return m_h1;}
196
200 void FitXT()
201 {
202 if (m_mode == 0) FitPol5();
203 else FitChebyshev();
204 }
205
209 void FitXT(int mode)
210 {
211 if (mode == 0) FitPol5();
212 else FitChebyshev();
213 }
217 void FitPol5();
221 void FitChebyshev();
222private:
223
224 TProfile* m_h1;
225 TF1* xtpol5 = new TF1("xtpol5", pol5pol1, 0.0, 700, 8);
226 TF1* xtCheb5 = new TF1("xtCheb5", Cheb5pol1, 0.0, 700, 8);
228 int m_mode = 1;
229 bool m_debug = true;
230 bool m_draw = false;
231 bool m_BField = true;
233 double m_XTParam[8] = {};
234 double m_FittedXTParams[8] = {};
244 int m_fitflag = 0;
245 double m_Prob = 0;
246 double m_tmin = 20;
247 double m_tmax = m_XTParam[6] + 250;
248};
249
251{
252 double max_dif = 0.12;
253 double max_dif2 = 0.05;
254 double par[8];
255 m_h1->Fit("pol1", "MQ", "", m_tmin, 50);
256 TF1* f1 = (TF1*)m_h1->GetFunction("pol1");
257 double p0 = f1->GetParameter(0);
258 double p1 = f1->GetParameter(1);
259 double f10 = f1->Eval(10);
260 /****************************/
261 int out = 0; /*how many time outer part change fit limit*/
262 xtpol5->SetParameters(p0, p1, 0, 0, 0, 0, m_XTParam[6], 0);
263 double p6default = m_XTParam[6];
264 if (m_BField) {
265 xtpol5->SetParLimits(7, 0.0, 0.001);
266 xtpol5->SetParLimits(1, 0.0, 0.01);
267 } else {
268 xtpol5->SetParLimits(7, 0.0, 0.01);
269 xtpol5->SetParLimits(1, 0.0, 0.01);
270 }
271 // xtpol5->SetParLimits(6, p6default - 30, p6default + 30);
272 for (int i = 0; i < 10; ++i) {
273 m_fitflag = 1;
274 std::cout << "Fitting" << std::endl;
275 double stat = m_h1->Fit(xtpol5, "M", "0", m_tmin, m_tmax);
276 xtpol5->GetParameters(par);
277
278 /*Eval outer region,*/
279 double fp6 = xtpol5->Eval(par[6]);
280 double fbehindp6 = xtpol5->Eval(par[6] - 12) - 0.01;
281 if (fp6 < fbehindp6 || fp6 > 1) { /*may be change to good value*/
282 m_fitflag = 2;
283 out += 1;
284 xtpol5->SetParameters(p0, p1, 0, 0, 0, 0, p6default, 0);
285 xtpol5->SetParLimits(6, p6default - 10, p6default + 10 - out / 2);
286 m_tmax += 2;
287 if (m_tmax < p6default + 30) {
288 m_tmax = p6default + 30;
289 }
290 }
291 /* EVal inner region*/
292 /* compare p0 of linear fit vs p0 of xtfun fit and eval point t=10)*/
293 if (fabs(par[0] - p0) > max_dif || fabs(f10 - xtpol5->Eval(10)) > max_dif2) {
294 m_fitflag = 3;
295 if (i == 9) std::cout << "ERROR XT FIT inner part" << std::endl;
296 xtpol5->SetParameters(p0, p1, 0, 0, 0, 0, p6default, 0);
297 xtpol5->SetParLimits(1, 0, 0.08);
298 m_tmin -= 0.5;
299
300 if (m_tmin < 14) {
301 m_tmin = 14; std::cout << "ERROR: tmin so small, fit iter: " << std::endl;
302 }
303 }
304
305 if (stat == 0) {
306 m_fitflag = 0;
307 if (m_debug) {std::cout << "FIT Failure; Layer: " << std::endl;}
308 }
309
310 if (m_debug) {printf("P6default= %3.2f, p6 fit = %3.2f", p6default, par[6]);}
311
312 if (m_fitflag == 1) {
313 if (m_debug) std::cout << "Fit success" << std::endl;
314 xtpol5->GetParameters(m_FittedXTParams);
315 m_Prob = xtpol5->GetProb();
316 break;
317 }
318
319 } //end loop of fitting
320 if (m_draw) {
321 TString hname = m_h1->GetName();
322 TString name = hname + ".pdf";
323 TCanvas* c1 = new TCanvas("c1", "", 800, 600);
324 m_h1->Draw();
325 xtpol5->SetLineColor(kBlack);
326 xtpol5->DrawF1(0, 400, "same");
327 c1->SaveAs(name);
328
329 }
330}
331
333{
334 if (m_h1->GetEntries() < m_minRequiredEntry) {
335 m_fitflag = -1;
336 return;
337 }
338 // m_tmax = m_XTParam[6] + 100;
339 //xtCheb5->SetParameters(0.0, 0.005, 0., 0., 0., 0., m_XTParam[6], 0.001);
340 double p[6];
341 double par[8];
342 xtCheb5->SetParLimits(7, 0., 0.001);
343 m_h1->Fit("chebyshev5", "QME", "", m_tmin, m_XTParam[6]);
344 m_h1->GetFunction("chebyshev5")->GetParameters(p);
345 xtCheb5->SetParameters(p[0], p[1], p[2], p[3], p[4], p[5], m_XTParam[6], 0.000);
346
347 double stat;
348 for (int i = 0; i < 10; ++i) {
349 stat = m_h1->Fit("xtCheb5", "MQ", "0", m_tmin, m_tmax);
350 if (stat == 0) {
351 xtCheb5->SetParameters(p[0], p[1], p[2], p[3], p[4], p[5], m_XTParam[6] - 20, 0.000);
352 m_tmax -= 10;
353 continue;
354 }
355 xtCheb5->GetParameters(par);
356 /*Eval outer region,*/
357 double fp6 = xtCheb5->Eval(par[6]);
358 double fbehindp6 = xtCheb5->Eval(par[6] - 10) - 0.005;
359 if (fp6 < fbehindp6 || fp6 > 1) { /*may be change to good value*/
360 m_fitflag = 2;
361 // out += 1;
362 xtCheb5->SetParameters(p[0], p[1], p[2], p[3], p[4], p[5], par[6] - 20, 0.000);
363 xtCheb5->SetParLimits(6, par[6] - 50, par[6] - 10);
364 m_tmax -= 10;
365 // if (m_tmax < p6default + 30) {
366 // m_tmax = p6default + 30;
367 // }
368 } else {
369 break;
370 }
371 // m_tmax +=10;
372 //if (stat != 0) break;
373 }
374 xtCheb5->GetParameters(m_FittedXTParams);
375 m_Prob = xtCheb5->GetProb();
376 if (stat == 0)
377 m_fitflag = 0;
378 else
379 m_fitflag = 1;
380 // xtCheb5->DrawF1(0,400,"same");
381 if (m_draw) {
382 TString hname = m_h1->GetName();
383 TString name = hname + ".pdf";
384 TCanvas* c1 = new TCanvas("c1", "", 800, 600);
385 m_h1->Draw();
386 xtCheb5->SetLineColor(kBlack);
387 xtCheb5->DrawF1(0, 400, "same");
388 c1->SaveAs(name);
389 }
390}
Class to perform fitting for each xt function.
Definition: XT.h:60
bool m_BField
With magnetic field or not.
Definition: XT.h:231
void setXTParams(const double p[8])
Set Parameters for fit.
Definition: XT.h:115
int m_minRequiredEntry
Minimum entry required for each histo.
Definition: XT.h:232
double m_FittedXTParams[8]
Fitted parameters.
Definition: XT.h:234
XT(TProfile *h1)
Initialized with TProfile histogram.
Definition: XT.h:65
void setP6(double p6)
Set Parameter 6 for polynomial fit.
Definition: XT.h:93
void setFitRange(double tmin, double tmax)
Set Fit range.
Definition: XT.h:135
void FitChebyshev()
Fit xt histogram incase 5th order Chebeshev polynomial is used.
Definition: XT.h:332
void FitXT()
Do fitting.
Definition: XT.h:200
double getProb()
Get the chi2 probability.
Definition: XT.h:165
TF1 * getXTFunction()
Get XT function.
Definition: XT.h:187
XT(TH1D *h1, int mode)
Initialized with TH1D histogram and mode.
Definition: XT.h:83
int m_mode
XT mode, 0 is for 5th order polynomial, 1 is Chebshev polynomial.
Definition: XT.h:228
int m_fitflag
Fit Flag =-1: low statitic =1: good =0: Fit failure =2: Error Outer =3: Error Inner part;.
Definition: XT.h:244
bool m_debug
Print debug durring fitting or not.
Definition: XT.h:229
TF1 * xtpol5
5th order polynomial function
Definition: XT.h:225
void setDebug(bool debug)
Set Debug.
Definition: XT.h:150
void FitXT(int mode)
Do fitting.
Definition: XT.h:209
double m_Prob
Chi2 prob of fitting.
Definition: XT.h:245
void setMode(int mode)
Set XT mode.
Definition: XT.h:103
TF1 * xtCheb5
5th order Cheb.
Definition: XT.h:226
void setSmallestEntryRequired(int min)
Set minimum number of entry required for fit.
Definition: XT.h:143
double m_tmin
lower boundary of fit range
Definition: XT.h:246
double m_tmax
upper boundary of fit range
Definition: XT.h:247
void getFittedXTParams(double pa[8])
get fit parameters.
Definition: XT.h:172
void setXTParams(double p0, double p1, double p2, double p3, double p4, double p5, double p6, double p7)
Set Initial parameters for fitting.
Definition: XT.h:123
TProfile * m_h1
Histogram of xt relation.
Definition: XT.h:224
TF1 * getXTFunction(int mode)
Get XT function.
Definition: XT.h:179
bool m_draw
Draw and store png plot of each histo or not.
Definition: XT.h:230
double m_XTParam[8]
Parameter fo xt.
Definition: XT.h:233
XT(TProfile *h1, int mode)
Initialized with TProfile histogram and mode.
Definition: XT.h:73
TProfile * getFittedHisto()
Get histogram.
Definition: XT.h:195
void FitPol5()
Fit xt histogram incase 5th order polynomial is used.
Definition: XT.h:250
int getFitStatus()
get fitted flag.
Definition: XT.h:157
void BField(bool bfield)
set to use BField
Definition: XT.h:110