pxar
 All Classes Namespaces Functions Variables Typedefs Friends
PixInitFunc.cc
1 #include "PixInitFunc.hh"
2 
3 #include "TROOT.h"
4 #include "TMath.h"
5 #include "TCanvas.h"
6 
7 #include <iostream>
8 #include <iomanip>
9 
10 
11 ClassImp(PixInitFunc)
12 
13 using namespace std;
14 
15 namespace {
16 
17  // ----------------------------------------------------------------------
18  // par[0]: "step"
19  // par[1]: "slope" the smaller the steeper
20  // par[2]: "floor" 1 -> floor is at 0
21  // par[3]: "plateau" y -> plateau is at 2*y
22 
23  double PIF_err(double *x, double *par) {
24  return par[3]*(TMath::Erf((x[0]-par[0])/par[1])+par[2]);
25  }
26 
27  double PIF_weibullCdf(double *x, double *par) {
28  if (par[1] < 0) return 0.;
29  if (par[2] < 0) return 0.;
30  return par[0]*(1 - TMath::Exp(-TMath::Power(x[0]/par[1], par[2])));
31  }
32 
33 
34 
35  const double xCut = TMath::Pi()/2. - 0.0005;
36  const double tanXCut = TMath::Tan(xCut);
37 
38  double PIF_gpTanPol( Double_t *x, Double_t *par) {
39  if (par[0]*x[0] - par[4] > xCut) return tanXCut + (x[0] - (xCut + par[4])/par[0])* 1e8;
40  double result = TMath::Tan(par[0]*x[0] - par[4]) + par[1]*x[0]*x[0]*x[0] + par[5]*x[0]*x[0] + par[2]*x[0] + par[3];
41  // cout << result << endl;
42  return result;
43  }
44 
45 
46  double PIF_gpTanH( Double_t *x, Double_t *par) {
47  return par[3] + par[2] * TMath::TanH(par[0]*x[0] - par[1]);
48  }
49 
50 }
51 
52 
53 // ----------------------------------------------------------------------
54 PixInitFunc::PixInitFunc() {
55 
56 }
57 
58 // ----------------------------------------------------------------------
59 PixInitFunc::~PixInitFunc() {
60 
61 }
62 
63 // ----------------------------------------------------------------------
64 void PixInitFunc::resetLimits() {
65  for (int i = 0; i < 20; ++i) {
66  fLimit[i] = false;
67  fLimitLo[i] = 0.;
68  fLimitHi[i] = 0.;
69  }
70 }
71 
72 // ----------------------------------------------------------------------
73 void PixInitFunc::limitPar(int ipar, double lo, double hi) {
74  fLimit[ipar] = true;
75  fLimitLo[ipar] = lo;
76  fLimitHi[ipar] = hi;
77 }
78 
79 
80 // ----------------------------------------------------------------------
81 TF1* PixInitFunc::gpTanPol(TH1 *h) {
82 
83  double lo = h->GetBinLowEdge(1);
84  double hi = h->FindLastBinAbove(0.9*h->GetMaximum());
85  hi = h->GetBinLowEdge(h->GetNbinsX()+1);
86  // -- setup function
87  TF1* f = (TF1*)gROOT->FindObject("PIF_gpTanPol");
88  if (0 == f) {
89  f = new TF1("PIF_gpTanPol", PIF_gpTanPol, lo, hi, 6);
90  f->SetNpx(1000);
91  // f->SetParNames("norm", "scale", "shape");
92  f->SetRange(lo, hi);
93  } else {
94  f->ReleaseParameter(0);
95  f->ReleaseParameter(1);
96  f->ReleaseParameter(2);
97  f->ReleaseParameter(3);
98  f->ReleaseParameter(4);
99  f->ReleaseParameter(5);
100  f->SetRange(lo, hi);
101  }
102 
103  int ilo = h->FindFirstBinAbove(0);
104  double xlo = h->GetBinCenter(ilo);
105  double ylo = h->GetBinContent(ilo);
106 
107  int iup = h->FindFirstBinAbove(0.8*h->GetMaximum());
108  double xup = h->GetBinCenter(iup);
109  double yup = h->GetBinContent(iup);
110 
111  double slope(0.5);
112  if (ilo != iup && (TMath::Abs(xup - xlo) > 1e-5)) slope = (yup-ylo)/(xup-xlo);
113 
114  f->SetParameter(2, slope);
115  f->SetParameter(3, yup - slope*xup);
116 
117  double par0 = (TMath::Pi()/2. - 1.4) / h->GetBinCenter(h->FindLastBinAbove(0.));
118  f->SetParameter(0, par0);
119  f->FixParameter(1, 0.);
120  f->SetParameter(4, -1.4);
121 
122  if (xup > 0.) f->SetParameter(5, (yup - (TMath::Tan(f->GetParameter(0)*xup - f->GetParameter(4))
123  + f->GetParameter(1)*xup*xup*xup + slope*xup + f->GetParameter(3)))/(xup*xup));
124  else f->SetParameter(5, 0.);
125 
126 // cout << f->GetParameter(0) << endl;
127 // cout << f->GetParameter(1) << endl;
128 // cout << f->GetParameter(2) << endl;
129 // cout << f->GetParameter(3) << endl;
130 // cout << f->GetParameter(4) << endl;
131 // cout << f->GetParameter(5) << endl;
132 
133  return f;
134 }
135 
136 
137 // ----------------------------------------------------------------------
138 TF1* PixInitFunc::gpTanH(TH1 *h) {
139 
140  double lo = 0.; //h->GetBinLowEdge(1);
141  double hi = 1700.; //h->FindLastBinAbove(0.9*h->GetMaximum());
142 
143  // -- setup function
144  TF1* f = (TF1*)gROOT->FindObject("PIF_gpTanH");
145  if (0 == f) {
146  f = new TF1("PIF_gpTanH", PIF_gpTanH, lo, hi, 4);
147  f->SetNpx(1000);
148  f->SetRange(lo, hi);
149  } else {
150  f->ReleaseParameter(0);
151  f->ReleaseParameter(1);
152  f->ReleaseParameter(2);
153  f->ReleaseParameter(3);
154  f->SetRange(lo, hi);
155  }
156 
157  f->SetParameter(0, 1.4e-3);
158  f->SetParLimits(0, 1.e-3, 2.e-3);
159 
160  f->SetParameter(1, 0.8);
161  f->SetParLimits(1, 0., 20.);
162 
163  double middle = h->GetMaximum() - h->GetMinimum();
164  double step = h->GetMaximum() - middle;
165 
166  f->SetParameter(2, middle);
167  f->SetParLimits(2, 0., 2*middle);
168 
169  f->SetParameter(3, step);
170 
171  return f;
172 }
173 
174 
175 // ----------------------------------------------------------------------
176 TF1* PixInitFunc::weibullCdf(TH1 *h) {
177  fDoNotFit = false;
178 
179  double hi = h->FindLastBinAbove(0.9*h->GetMaximum());
180  double lo = h->GetBinLowEdge(1);
181  // require 3 consecutive bins at zero
182  for (int i = 3; i < h->GetNbinsX(); ++i) {
183  if (h->GetBinContent(i-2) < 1 && h->GetBinContent(i-1) < 1 && h->GetBinContent(i) < 1) {
184  lo = h->GetBinLowEdge(i-2);
185  break;
186  }
187  }
188 
189  // -- setup function
190  TF1* f = (TF1*)gROOT->FindObject("PIF_weibullCdf");
191  if (0 == f) {
192  f = new TF1("PIF_weibullCdf", PIF_weibullCdf, h->GetBinLowEdge(1), h->GetBinLowEdge(h->GetNbinsX()+1), 3);
193  f->SetParNames("norm", "scale", "shape");
194  f->SetNpx(1000);
195  f->SetRange(lo, hi);
196  } else {
197  f->ReleaseParameter(0);
198  f->ReleaseParameter(1);
199  f->ReleaseParameter(2);
200  f->ReleaseParameter(3);
201  f->SetRange(lo, hi);
202  }
203 
204  f->SetParameter(0, h->GetMaximum());
205  f->SetParameter(1, h->GetBinCenter(h->FindFirstBinAbove(0.5*h->GetMaximum())));
206  f->SetParameter(2, 1.5);
207 
208  return f;
209 }
210 
211 
212 // ----------------------------------------------------------------------
213 TF1* PixInitFunc::errScurve(TH1 *h) {
214 
215  fDoNotFit = false;
216 
217  // -- determine step function and start of function range (to exclude spurious low-threshold readouts)
218  int STARTBIN(2);
219  int ibin(-1), jbin(-1);
220  double hmax(h->GetMaximum());
221  for (int i = STARTBIN; i <= h->GetNbinsX(); ++i) {
222  if (h->GetBinContent(i) > 0) {
223  ibin = i; // first bin above zero
224  break;
225  }
226  }
227 
228  // require 2 consecutive bins on plateau
229  for (int i = STARTBIN; i < h->GetNbinsX(); ++i) {
230  if (h->GetBinContent(i) > 0.9*hmax && h->GetBinContent(i+1) > 0.9*hmax) {
231  jbin = i; // first bin "really" on plateau
232  break;
233  }
234  }
235 
236  int plateauWidth = h->FindLastBinAbove(0.9*hmax) - jbin;
237  double hi = h->GetBinLowEdge(jbin + 0.7*plateauWidth);
238 
239  double lo = h->GetBinLowEdge(1);
240  // require 3 consecutive bins at zero
241  int foundLo(-1);
242  for (int i = 3; i < h->GetNbinsX(); ++i) {
243  if (h->GetBinLowEdge(i-2) > hi) break;
244  if (h->GetBinContent(i-2) < 1 && h->GetBinContent(i-1) < 1 && h->GetBinContent(i) < 1) {
245  lo = h->GetBinLowEdge(i-2);
246  foundLo = i-2;
247  break;
248  }
249  }
250  if (foundLo > -1) {
251  lo = lo + 0.6*(h->GetBinLowEdge(ibin) - lo);
252  }
253 
254  // -- setup function
255  TF1* f = (TF1*)gROOT->FindObject("PIF_err");
256  if (0 == f) {
257  f = new TF1("PIF_err", PIF_err, h->GetBinLowEdge(1), h->GetBinLowEdge(h->GetNbinsX()+1), 4);
258  f->SetParNames("step", "slope", "floor", "plateau");
259  f->SetNpx(1000);
260  f->SetRange(lo, hi);
261  } else {
262  f->ReleaseParameter(0);
263  f->ReleaseParameter(1);
264  f->ReleaseParameter(2);
265  f->ReleaseParameter(3);
266  f->SetRange(lo, hi);
267  }
268 
269  f->SetParameter(0, h->GetBinCenter(0.5*(ibin+jbin)));
270  f->SetParameter(1, 2.0/(h->GetBinLowEdge(jbin) - h->GetBinLowEdge(ibin)));
271  // cout << "init parameters 0: " << f->GetParameter(0) << " 1: " << f->GetParameter(1)
272  // << " ibin: " << ibin << " jbin: " << jbin
273  // << " lo: " << lo << " hi: " << hi
274  // << " foundLo: " << foundLo
275  // << endl;
276  if (jbin == ibin) {
277  // cout << "XXXXXXXXXXX PixInitFunc: STEP FUNCTION " << h->GetTitle() << " ibin = " << ibin << " jbin = " << jbin << endl;
278  f->FixParameter(0, h->GetBinCenter(jbin));
279  f->FixParameter(1, 1.e2);
280  fDoNotFit = true;
281  }
282  f->FixParameter(2, 1.);
283  f->FixParameter(3, 0.5*h->GetMaximum());
284  return f;
285 }
286 
287 // ----------------------------------------------------------------------
288 TF1* PixInitFunc::xrayScan(TH1 *h) {
289 
290  // -- determine step function and start of function range
291  double lo = h->GetBinCenter(h->FindFirstBinAbove(0.));
292  double hi = h->GetBinCenter(h->FindLastBinAbove(0.9*h->GetMaximum()));
293  double hmax = h->GetMaximum();
294 
295  double pol0(0.);
296  double threshold(-1.), plateau(-1), width(-1.);
297  for (int i = h->FindBin(lo); i < h->GetNbinsX(); i += 5) {
298  if (h->GetBinCenter(i) > hi - 10.) break;
299  h->Fit("pol0", "qr", "", h->GetBinCenter(i), hi);
300  pol0 = h->GetFunction("pol0")->GetParameter(0);
301  // cout << h->GetBinCenter(i) << " pol0 = " << pol0 << endl;
302 
303  double a = pol0/hmax;
304  if (a > 0.5 && threshold < 0) {
305  threshold = h->GetBinCenter(i);
306  }
307  plateau = pol0;
308  }
309 
310  width = 0.5*(threshold - lo);
311 
312  // cout << "lo = " << lo << endl;
313  // cout << "hi = " << hi << endl;
314  // cout << "threshold = " << threshold << endl;
315  // cout << "plateau = " << plateau << endl;
316 
317  // -- setup function
318  TF1* f = (TF1*)gROOT->FindObject("PIF_err");
319  if (0 == f) {
320  f = new TF1("PIF_err", PIF_err, h->GetBinLowEdge(1), h->GetBinLowEdge(h->GetNbinsX()+1), 4);
321  f->SetParNames("step", "slope", "floor", "plateau");
322  f->SetNpx(1000);
323  f->SetRange(lo, hi);
324  } else {
325  f->ReleaseParameter(0);
326  f->ReleaseParameter(1);
327  f->ReleaseParameter(2);
328  f->ReleaseParameter(3);
329  f->SetRange(lo, hi);
330  }
331 
332  f->SetParameter(0, threshold);
333  f->SetParameter(1, width);
334  f->SetParameter(2, 1.);
335  f->SetParameter(3, 0.5*plateau);
336  return f;
337 }
338 
339 
340 
341 
342 // ----------------------------------------------------------------------
343 void PixInitFunc::initExpo(double &p0, double &p1, TH1 *h) {
344 
345  int EDG(4), NB(EDG+1);
346  int lbin(1), hbin(h->GetNbinsX()+1);
347  if (fLo < fHi) {
348  lbin = h->FindBin(fLo);
349  hbin = h->FindBin(fHi);
350  }
351 
352  double dx = h->GetBinLowEdge(hbin) - h->GetBinLowEdge(lbin);
353  double ylo = h->Integral(lbin, lbin+EDG)/NB;
354  double yhi = h->Integral(hbin-EDG, hbin)/NB;
355 
356  if (ylo > 0 && yhi > 0) {
357  p1 = (TMath::Log(yhi) - TMath::Log(ylo))/dx;
358  p0 = ylo/TMath::Exp(p1*fLo);
359  } else {
360  if (yhi > ylo) {
361  p1 = 1.;
362  } else {
363  p1 = -1.;
364  }
365  p0 = 50.;
366  }
367 
368 // cout << "fLo: " << fLo << " fHi: " << fHi << endl;
369 // cout << "ylo: " << ylo << " yhi: " << yhi << endl;
370 // cout << "p0: " << p0 << " p1: " << p1 << endl;
371 
372 }
373 
374 // ----------------------------------------------------------------------
375 void PixInitFunc::initPol1(double &p0, double &p1, TH1 *h) {
376 
377  int EDG(4), NB(EDG+1);
378  int lbin(1), hbin(h->GetNbinsX()+1);
379  if (fLo < fHi) {
380  lbin = h->FindBin(fLo);
381  hbin = h->FindBin(fHi);
382  }
383 
384  double dx = h->GetBinLowEdge(hbin) - h->GetBinLowEdge(lbin);
385  double ylo = h->Integral(lbin, lbin+EDG)/NB;
386  double yhi = h->Integral(hbin-EDG, hbin)/NB;
387 
388  // cout << "ylo: " << ylo << " yhi: " << yhi << endl;
389 
390  p1 = (yhi-ylo)/dx;
391  p0 = ylo - p1*fLo;
392 }
393 
Generated on Tue Mar 3 2015 10:03:36 for pxar by   doxygen 1.8.6