Belle II Software  release-05-01-25
XTFunction.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Makoto Uchida *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include "TProfile.h"
11 #include "TCanvas.h"
12 #include "TF1.h"
13 #include "TH1D.h"
14 #include "Math/ChebyshevPol.h"
15 #include "iostream"
16 
17 
18 namespace Belle2 {
23  namespace CDC {
24 
29  Double_t pol5pol1(Double_t* x, Double_t* par)
30  {
31  Double_t xx = x[0];
32  Double_t x6, x2;
33  Double_t f, ctp;
34  if (xx < par[6]) {
35  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;
36  } else {
37  x6 = par[6];
38  x2 = xx - x6;
39  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 *
40  x6;
41  f = par[7] * x2 + ctp;
42  }
43  return f;
44  }
45 
50  Double_t cheby5pol1(Double_t* x, Double_t* par)
51  {
52  Double_t xx = x[0];
53  Double_t x6, x2;
54  Double_t f, ctp;
55 
56  if (xx < par[6]) {
57  f = ROOT::Math::Chebyshev5(xx, par[0], par[1], par[2], par[3], par[4], par[5]);
58  } else {
59  x6 = par[6];
60  x2 = xx - x6;
61  ctp = ROOT::Math::Chebyshev5(x6, par[0], par[1], par[2], par[3], par[4], par[5]);
62  f = par[7] * x2 + ctp;
63  }
64  return f;
65  }
66 
71  class XTFunction {
72  public:
76  explicit XTFunction(TProfile* h1)
77  {
78  m_h1 = (TProfile*)h1->Clone();
79  m_h1->SetDirectory(0);
80  if (m_mode == c_Chebyshev) {
81  m_fitFunc = new TF1("xtCheb5", cheby5pol1, 0.0, 700, 8);
82  } else {
83  m_fitFunc = new TF1("xtpol5", pol5pol1, 0.0, 700, 8);
84  }
85  }
89  XTFunction(TProfile* h1, int mode)
90  {
91  m_h1 = (TProfile*)h1->Clone();
92  m_h1->SetDirectory(0);
93  m_mode = mode;
94  if (m_mode == c_Chebyshev) {
95  m_fitFunc = new TF1("xtCheb5", cheby5pol1, 0.0, 700, 8);
96  } else {
97  m_fitFunc = new TF1("xtpol5", pol5pol1, 0.0, 700, 8);
98  }
99 
100  }
101 
105  XTFunction(TH1F* h1, int mode)
106  {
107  m_h1 = (TProfile*)h1->Clone();
108  m_h1->SetDirectory(0);
109  m_mode = mode;
110  if (m_mode == c_Chebyshev) {
111  m_fitFunc = new TF1("xtCheb5", cheby5pol1, 0.0, 700, 8);
112  } else {
113  m_fitFunc = new TF1("xtpol5", pol5pol1, 0.0, 700, 8);
114  }
115  }
116 
120  XTFunction(const XTFunction& x) :
121  m_h1(x.m_h1),
122  m_mode(x.m_mode),
123  m_debug(x.m_debug),
124  m_draw(x.m_draw),
125  m_bField(x.m_bField),
127  m_fitflag(x.m_fitflag),
129  m_tmin(x.m_tmin),
130  m_tmax(x.m_tmax)
131  {
132  m_fitFunc = (TF1*) x.m_fitFunc->Clone();
133  for (int i = 0; i < 8; ++i) {
134  m_XTParam[i] = x.m_XTParam[i];
135  m_FittedXTParams[i] = x.m_XTParam[i];
136  }
137  }
138 
142  XTFunction& operator=(const XTFunction& x)
143  {
144  if (this != &x) {
145  m_h1 = x.m_h1;
146  m_fitFunc = x.m_fitFunc;
147  m_mode = x.m_mode;
148  m_debug = x.m_debug;
149  m_draw = x.m_draw;
150  m_bField = x.m_bField;
151  m_minRequiredEntry = x.m_minRequiredEntry;
152  m_fitflag = x.m_fitflag;
153  m_Prob = x.m_Prob;
154  m_tmin = x.m_tmin;
155  m_tmax = x.m_tmax;
156 
157  for (int i = 0; i < 8; ++i) {
158  m_XTParam[i] = x.m_XTParam[i];
159  m_FittedXTParams[i] = x.m_XTParam[i];
160  }
161  }
162  return *this;
163  }
164 
165 
169  void setP6(double p6)
170  {
171  m_XTParam[6] = p6;
172  }
173 
177  bool isValid()
178  {
179  if (m_fitFunc->IsValid() == true) {
180  return true;
181  } else {
182  return false;
183  }
184  }
185 
191  void setMode(int mode)
192  {
193  m_mode = mode;
194  }
198  void setBField(bool bfield) {m_bField = bfield;}
199 
203  void setXTParams(double p[8])
204  {
205  for (int i = 0; i < 8; ++i) {m_XTParam[i] = p[i];}
206  m_tmax = p[6] + 50;
207  }
211  void setXTParams(double p0, double p1, double p2, double p3,
212  double p4, double p5, double p6, double p7)
213  {
214  m_XTParam[0] = p0; m_XTParam[1] = p1; m_XTParam[2] = p2;
215  m_XTParam[3] = p3; m_XTParam[4] = p4; m_XTParam[5] = p5;
216  m_XTParam[6] = p6; m_XTParam[7] = p7;
217  m_tmax = p6 + 50;
218  }
219 
223  void setFitRange(double tmin, double tmax)
224  {
225  m_tmin = tmin;
226  m_tmax = tmax;
227  }
232  {
233  m_minRequiredEntry = min;
234  }
238  void setDebug(bool debug)
239  {
240  m_debug = debug;
241  }
245  int getFitStatus()
246  {
247  return m_fitflag;
248  }
249 
253  double getProb()
254  {
255  return m_Prob;
256  }
260  void getFittedXTParams(double pa[8])
261  {
262  for (int i = 0; i < 8; ++i) {pa[i] = m_FittedXTParams[i];}
263  }
267  TF1* getXTFunction()
268  {
269  return m_fitFunc;
270  }
274  TProfile* getFittedHisto() {return m_h1;}
275 
279  void fitXT()
280  {
281  if (m_mode == c_Polynomial) {
283  } else if (m_mode == c_Chebyshev) {
284  FitChebyshev();
285  } else {
286  B2ERROR("Undefined fitting function");
287  }
288  }
289 
293  void FitPol5();
297  void FitChebyshev();
302  bool validate();
303  private:
304 
305  TProfile* m_h1;
306  TF1* m_fitFunc;
308  // TF1* xtpol5 = new TF1("xtpol5", pol5pol1, 0.0, 700, 8); /**< 5th order polynomial function*/
309  // TF1* xtCheb5 = new TF1("xtCheb5", Cheb5pol1, 0.0, 700, 8); /**< 5th order Cheb. polynomial function*/
310 
311  int m_mode = c_Chebyshev;
312  bool m_debug = true;
313  bool m_draw = false;
314  bool m_bField = true;
315  int m_minRequiredEntry = 800;
316  double m_XTParam[8] = {};
317  double m_FittedXTParams[8] = {};
327  int m_fitflag = 0;
328  double m_Prob = 0;
329  double m_tmin = 20;
330  double m_tmax = m_XTParam[6] + 50;
331  };
332 
333  void XTFunction::FitPol5()
334  {
335  if (m_mode != c_Polynomial) {
336  B2ERROR("Fitting function is wrong");
337  }
338  double max_dif = 0.12;
339  double max_dif2 = 0.05;
340  double par[8];
341  m_h1->Fit("pol1", "MQ", "", m_tmin, 50);
342  TF1* f1 = (TF1*)m_h1->GetFunction("pol1");
343  double p0 = f1->GetParameter(0);
344  double p1 = f1->GetParameter(1);
345  double f10 = f1->Eval(10);
346  /****************************/
347  int in = 0; /*how many time inner part change fit limit*/
348  int out = 0; /*how many time outer part change fit limit*/
349  m_fitFunc->SetParameters(p0, p1, 0, 0, 0, 0, m_XTParam[6], 0);
350  double p6default = m_XTParam[6];
351  if (m_bField) {
352  m_fitFunc->SetParLimits(7, 0.0, 0.001);
353  m_fitFunc->SetParLimits(1, 0.0, 0.01);
354  } else {
355  m_fitFunc->SetParLimits(7, 0.0, 0.01);
356  m_fitFunc->SetParLimits(1, 0.0, 0.01);
357  }
358 
359  for (int i = 0; i < 10; ++i) {
360  m_fitflag = 1;
361  std::cout << "Fitting" << std::endl;
362  double stat = m_h1->Fit(m_fitFunc, "M", "0", m_tmin, m_tmax);
363  m_fitFunc->GetParameters(par);
364 
365  /*Eval outer region,*/
366  double fp6 = m_fitFunc->Eval(par[6]);
367  double fbehindp6 = m_fitFunc->Eval(par[6] - 12) - 0.01;
368  if (fp6 < fbehindp6 || fp6 > 1) { /*may be change to good value*/
369  m_fitflag = 2;
370  out += 1;
371  m_fitFunc->SetParameters(p0, p1, 0, 0, 0, 0, p6default, 0);
372  m_fitFunc->SetParLimits(6, p6default - 10, p6default + 10 - out / 2);
373  m_tmax += 2;
374  if (m_tmax < p6default + 30) {
375  m_tmax = p6default + 30;
376  }
377  }
378  /* EVal inner region*/
379  /* compare p0 of linear fit vs p0 of xtfun fit and eval point t=10)*/
380  if (fabs(par[0] - p0) > max_dif || fabs(f10 - m_fitFunc->Eval(10)) > max_dif2) {
381  m_fitflag = 3;
382  if (i == 9) std::cout << "ERROR XT FIT inner part" << std::endl;
383  in += 1;
384  m_fitFunc->SetParameters(p0, p1, 0, 0, 0, 0, p6default, 0);
385  m_fitFunc->SetParLimits(1, 0, 0.08);
386  m_tmin -= 0.5;
387 
388  if (m_tmin < 14) {
389  m_tmin = 14; std::cout << "ERROR: tmin so small, fit iter: " << std::endl;
390  }
391  }
392 
393  if (stat == 0) {
394  m_fitflag = 0;
395  if (m_debug) {std::cout << "FIT Failure; Layer: " << std::endl;}
396  }
397 
398  if (m_debug) {printf("P6default= %3.2f, p6 fit = %3.2f", p6default, par[6]);}
399 
400  if (m_fitflag == 1) {
401  if (m_debug) std::cout << "Fit success" << std::endl;
402  m_fitFunc->GetParameters(m_FittedXTParams);
403  m_Prob = m_fitFunc->GetProb();
404  break;
405  }
406 
407  } //end loop of fitting
408  if (m_draw) {
409  TString hname = m_h1->GetName();
410  TString name = hname + ".pdf";
411  TCanvas* c1 = new TCanvas("c1", "", 800, 600);
412  m_h1->Draw();
413  m_fitFunc->SetLineColor(kBlack);
414  m_fitFunc->DrawF1(0, 400, "same");
415  c1->SaveAs(name);
416  }
417  }
418 
419  bool XTFunction::validate()
420  {
421 
422  const double p6 = m_fitFunc->GetParameter(6);
423  if (fabs(m_fitFunc->Eval(0)) > 0.2) {
424  B2WARNING("Bad xt function");
425  m_fitflag = 0;
426  return false;
427  } else if (p6 < 100.0) {
428  B2WARNING("Unrealistic p6");
429  m_fitflag = 0;
430  return false;
431  } else {
432  return true;
433  }
434  }
435 
437  {
438  if (m_mode != c_Chebyshev) {
439  B2ERROR("Fitting function is wrong");
440  }
441 
442  if (m_h1->GetEntries() < m_minRequiredEntry) {
443  m_fitflag = -1;
444  return;
445  }
446  // m_tmax = m_XTParam[6] + 100;
447  //xtCheb5->SetParameters(0.0, 0.005, 0., 0., 0., 0., m_XTParam[6], 0.001);
448  double par[8];
449  m_fitFunc->SetParLimits(7, 0., 0.001);
450  int fitresult = m_h1->Fit("chebyshev5", "QME", "", m_tmin, m_XTParam[6]);
451  if (fitresult >= 0) {
452  m_h1->GetFunction("chebyshev5")->GetParameters(par);
453  m_fitFunc->SetParameters(par[0], par[1], par[2], par[3], par[4], par[5], m_XTParam[6], 0.000);
454  }
455  double stat;
456  for (int i = 0; i < 10; ++i) {
457  stat = m_h1->Fit(m_fitFunc, "MQ", "0", m_tmin, m_tmax);
458  if (stat == 0) {
459  m_fitFunc->SetParameters(par[0], par[1], par[2], par[3], par[4], par[5], m_XTParam[6] - 20, 0.000);
460  m_tmax -= 10;
461  continue;
462  }
463  m_fitFunc->GetParameters(par);
464  if (par[1] < 0) { // negative c1
465  // std::cout << " neg c1 converted" << std::endl;
466  par[1] *= -1.0;
467  m_fitFunc->SetParLimits(1, 0., 0.01);
468  m_tmin += 10.0;
469  continue;
470  }
471 
472  /*Eval outer region,*/
473  double fp6 = m_fitFunc->Eval(par[6]);
474  double fbehindp6 = m_fitFunc->Eval(par[6] - 10) - 0.005;
475  if (fp6 < fbehindp6 || fp6 > 1) { /*may be change to good value*/
476  m_fitflag = 2;
477  // out += 1;
478  m_fitFunc->SetParameters(par[0], par[1], par[2], par[3], par[4], par[5], par[6] - 20, 0.000);
479  m_fitFunc->SetParLimits(6, par[6] - 50, par[6] - 10);
480  m_tmax -= 10;
481  // if (m_tmax < p6default + 30) {
482  // m_tmax = p6default + 30;
483  // }
484  } else {
485  break;
486  }
487  // m_tmax +=10;
488  //if (stat != 0) break;
489  }
490  if (par[1] > 0) {
491  m_fitFunc->GetParameters(m_FittedXTParams);
492  m_Prob = m_fitFunc->GetProb();
493  }
494  if (stat == 0)
495  m_fitflag = 0;
496  else
497  m_fitflag = 1;
498  // m_fitFunc->DrawF1(0,400,"same");
499  if (m_draw) {
500  TString hname = m_h1->GetName();
501  TString name = hname + ".pdf";
502  TCanvas* c1 = new TCanvas("c1", "", 800, 600);
503  m_h1->Draw();
504  m_fitFunc->SetLineColor(kBlack);
505  m_fitFunc->DrawF1(0, 400, "same");
506  c1->SaveAs(name);
507  }
508  }
509  }
511 }
Belle2::CDC::XTFunction::m_FittedXTParams
double m_FittedXTParams[8]
Fitted parameters.
Definition: XTFunction.h:325
Belle2::CDC::XTFunction::setXTParams
void setXTParams(double p[8])
Set Paramerters for fit.
Definition: XTFunction.h:211
Belle2::CDC::XTFunction::isValid
bool isValid()
Is valid.
Definition: XTFunction.h:185
Belle2::CDC::XTFunction::getFitStatus
int getFitStatus()
get fitted flag.
Definition: XTFunction.h:253
Belle2::CDC::XTFunction
Class to perform fitting for each xt function.
Definition: XTFunction.h:79
Belle2::CDC::XTFunction::m_debug
bool m_debug
Print debug durring fitting or not.
Definition: XTFunction.h:320
Belle2::CDC::XTFunction::setMode
void setMode(int mode)
Set XT mode.
Definition: XTFunction.h:199
Belle2::CDC::XTFunction::getFittedXTParams
void getFittedXTParams(double pa[8])
get fit parameters.
Definition: XTFunction.h:268
Belle2::CDC::XTFunction::setSmallestEntryRequired
void setSmallestEntryRequired(int min)
Set minimum number of entry required for fit.
Definition: XTFunction.h:239
Belle2::CDC::XTFunction::m_draw
bool m_draw
Draw and store png plot of each histo or not.
Definition: XTFunction.h:321
Belle2::CDC::XTFunction::setFitRange
void setFitRange(double tmin, double tmax)
Set Fit range.
Definition: XTFunction.h:231
Belle2::CDC::XTFunction::m_h1
TProfile * m_h1
Histogram of xt relation.
Definition: XTFunction.h:313
Belle2::CDC::cheby5pol1
Double_t cheby5pol1(Double_t *x, Double_t *par)
helper function to initialize xt function with 5th order Chebshev Polynomial + linear.
Definition: XTFunction.h:58
Belle2::CDC::XTFunction::getXTFunction
TF1 * getXTFunction()
Get XT function.
Definition: XTFunction.h:275
Belle2::CDC::XTFunction::XTFunction
XTFunction(TProfile *h1)
Initialized with TProfile histogram.
Definition: XTFunction.h:84
Belle2::CDC::XTFunction::setDebug
void setDebug(bool debug)
Set Debug.
Definition: XTFunction.h:246
Belle2::CDC::XTFunction::m_tmin
double m_tmin
lower boundary of fit range
Definition: XTFunction.h:337
Belle2::CDC::XTFunction::FitChebyshev
void FitChebyshev()
Fit xt histogram incase 5th order Chebeshev polynomial is used.
Definition: XTFunction.h:444
Belle2::CDC::XTFunction::m_fitFunc
TF1 * m_fitFunc
Fit function.
Definition: XTFunction.h:314
Belle2::CDC::XTFunction::operator=
XTFunction & operator=(const XTFunction &x)
Assignment operator.
Definition: XTFunction.h:150
Belle2::CDC::XTFunction::getProb
double getProb()
Get the chi2 probability.
Definition: XTFunction.h:261
Belle2::CDC::XTFunction::m_Prob
double m_Prob
Chi2 prob of fitting.
Definition: XTFunction.h:336
Belle2::CDC::XTFunction::m_mode
int m_mode
XT mode, 0 is for 5th order polynomial, 1 is Chebshev polynomial.
Definition: XTFunction.h:319
Belle2::CDC::XTFunction::m_minRequiredEntry
int m_minRequiredEntry
Minimum entry required for each histo.
Definition: XTFunction.h:323
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDC::XTFunction::validate
bool validate()
Validate the xt has proper shape.
Definition: XTFunction.h:427
Belle2::CDC::XTFunction::m_bField
bool m_bField
With magnetic field or not.
Definition: XTFunction.h:322
Belle2::CDC::XTFunction::m_XTParam
double m_XTParam[8]
Parameter fo xt.
Definition: XTFunction.h:324
Belle2::CDC::XTFunction::FitPol5
void FitPol5()
Fit xt histogram incase 5th order polynomial is used.
Definition: XTFunction.h:341
Belle2::CDC::XTFunction::setBField
void setBField(bool bfield)
set to use BField
Definition: XTFunction.h:206
Belle2::CDC::XTFunction::m_fitflag
int m_fitflag
Fit Flag =-1: low statitic =1: good =0: Fit failure =2: Error Outer =3: Error Inner part;.
Definition: XTFunction.h:335
Belle2::CDC::pol5pol1
Double_t pol5pol1(Double_t *x, Double_t *par)
helper function to initialize xt function with 5th order polynomial + linear.
Definition: XTFunction.h:37
Belle2::CDC::XTFunction::m_tmax
double m_tmax
upper boundary of fit range
Definition: XTFunction.h:338
Belle2::CDC::XTFunction::fitXT
void fitXT()
Do fitting.
Definition: XTFunction.h:287
Belle2::CDC::XTFunction::setP6
void setP6(double p6)
Set Parameter 6 for polynomia fit.
Definition: XTFunction.h:177
Belle2::CDC::XTFunction::getFittedHisto
TProfile * getFittedHisto()
Get histogram.
Definition: XTFunction.h:282