source: fact/tools/rootmacros/PulseTemplates/pulse.C

Last change on this file was 14962, checked in by Jens Buss, 12 years ago
bootstapping>: bin values extracted from bincenter
  • Property svn:executable set to *
File size: 16.4 KB
Line 
1#include <iostream>
2#include <fstream>
3#include <stdlib.h>
4
5#include "pulse.h"
6
7using namespace std;
8using namespace TMath;
9
10Pulse::Pulse()
11{
12 InitMembers();
13}
14
15Pulse::Pulse(TString name)
16{
17 InitMembers();
18 mName = name;
19}
20
21Pulse::Pulse(TString name, TH1* histo)
22{
23 InitMembers();
24 mName = name;
25 mHisto = histo;
26 Fit(mName, mOptions, mType);
27}
28
29Pulse::Pulse(TString name, TH1* histo, int type)
30{
31 InitMembers();
32 mName = name;
33 mHisto = histo;
34 mType = type;
35 Fit(mName, mOptions, mType);
36}
37
38Pulse::Pulse(TString name, TH1* histo, TString options)
39{
40 InitMembers();
41 mName = name;
42 mHisto = histo;
43 mOptions += options;
44 Fit(mName, mOptions, mType);
45}
46
47Pulse::Pulse(TString name, TH1* histo, TString options, int type)
48{
49 InitMembers();
50 mName = name;
51 mHisto = histo;
52 mOptions += options;
53 mType = type;
54 Fit(mName, mOptions, mType);
55}
56
57Pulse::Pulse(TString name, TH1* histo, TString options, int type, int order)
58{
59 InitMembers();
60 mName = name;
61 mHisto = histo;
62 mOptions += options;
63 mType = type;
64 mOrder = order;
65 Fit(mName, mOptions, mType);
66}
67
68
69Pulse::~Pulse()
70{
71 delete mModel;
72}
73
74// ===========================================================================
75// PUBLIC OPERATIONS
76// ===========================================================================
77
78void
79Pulse::Fit()
80{
81 Fit(mName, mOptions, mType);
82}
83
84void
85Pulse::Fit(
86 TString fitOptions,
87 int type
88 )
89{
90 Fit(mName, fitOptions, type);
91}
92
93void
94Pulse::Fit(
95 TString fitName,
96 TString fitOptions,
97 int type
98 )
99{
100 int fitMin = mHisto->GetXaxis()->GetBinLowEdge( mHisto->GetXaxis()->GetFirst() );
101 int fitMax = mHisto->GetXaxis()->GetBinUpEdge( mHisto->GetXaxis()->GetLast() );
102
103 Fit(
104 fitName,
105 fitOptions,
106 type,
107 fitMin,
108 fitMax
109 );
110}
111
112void
113Pulse::Fit(
114 TString fitName,
115 TString fitOptions,
116 int type,
117 int fitMin,
118 int fitMax
119 )
120{
121// mHisto->Sumw2();
122
123 fitOptions += "R";
124 if (type == 0) //use sectionwise function
125 FitSectionWise(
126 fitName,
127 fitOptions,
128 fitMin,
129 fitMax
130 );
131 else if (type == 1) //use continous function
132 FitContious(
133 fitName,
134 fitOptions,
135 fitMin,
136 fitMax
137 );
138 else{
139 cout << "not a correct type number --> fitting skipped!" << endl;
140 return;
141 }
142
143 CalculateParameters();
144 return;
145}
146
147// ===========================================================================
148// PRIVATE OPERATIONS
149// ===========================================================================
150
151void
152Pulse::InitMembers()
153{
154 mName = "pulse";
155 mHisto = NULL;
156 mModel = NULL;
157 mOptions = "ERS";
158 mBsl.first = 0;
159 mHeight.first = 0;
160 mT0.first = 0;
161 mT1.first = 0;
162 mTau1.first = 0;
163 mTau2.first = 0;
164 mIntegral.first = 0;
165 mAmplitude.first = 0;
166 mPhE.first = 0;
167 mBsl.second = 0;
168 mHeight.second = 0;
169 mT0.second = 0;
170 mT1.second = 0;
171 mTau1.second = 0;
172 mTau2.second = 0;
173 mIntegral.second = 0;
174 mAmplitude.second = 0;
175 mPhE.second = 0;
176 mRisetime.first = 0;
177 mRisetime.second = 0;
178 mOrder = 0;
179 mType = 0;
180 mFitProb = 0;
181 mFitNCalls = 0;
182 mFitNdf = 0;
183 mChi2 = 0;
184 mTimeResolution = 0.5;
185 mVerbosityLvl = 6;
186}
187
188// Fit pulseshape with a section wise defined funciton
189// expressed by a sum of section wise defined exp-functions
190//controled by heaviside function
191void
192Pulse::FitSectionWise(
193 TString fitName,
194 TString fitOptions,
195 int fitMin,
196 int fitMax
197 )
198{
199 int numPara = 7;
200 mFitMin = fitMin;
201 mFitMax = fitMax;
202 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, numPara );
203 mModel->SetNpx( mHisto->GetNbinsX());
204 // ======================================================================
205 // Calculate startvaues for fit
206 double tau = 30.0;
207 double bsl = 0.8; //shall be around 0
208 double gain = 10; //startvalue for gain
209
210 int first_bin = mHisto->GetXaxis()->GetFirst();
211
212 for (int i = 0; i < 30; i++){
213 bsl += mHisto->GetBinContent(first_bin+0);
214 }
215 bsl = bsl/30;
216
217 double stop = mHisto->GetMaximumBin(); //pos of max
218 double amplitude = mHisto->GetBinContent(stop);
219 double start = stop-10; //pos 10 slices before maximum
220 int phe = 0;
221 if (mOrder < 0)
222 {
223 phe = (amplitude-bsl)/gain;
224 }
225 else
226 {
227 phe = mOrder+1;
228 }
229 // ======================================================================
230
231 double para[] = {bsl, amplitude, start, stop, tau, tau, phe};
232
233
234 // cout startvalues
235 if (mVerbosityLvl > 3)
236 {
237 for (unsigned int i = 0; i < numPara; i++)
238 {
239 cout << "para[ " << i << "]: " << para[i] << endl;
240 }
241 }
242
243
244 mModel->SetParameters(para);
245 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2", "PhE");
246 mModel->SetLineColor(kRed);
247 mModel->FixParameter(0, bsl);
248
249 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
250 mBsl.first = mModel->GetParameter(0);
251 mBsl.second = mModel->GetParError(0);
252 mHeight.first = mModel->GetParameter(1);
253 mHeight.second = mModel->GetParError(1);
254 mT0.first = mModel->GetParameter(2);
255 mT0.second = mModel->GetParError(2);
256 mT1.first = mModel->GetParameter(3);
257 mT1.second = mModel->GetParameter(3);
258 mTau1.first = mModel->GetParameter(4);
259 mTau1.second = mModel->GetParError(4);
260 mTau2.first = mModel->GetParameter(5);
261 mTau2.second = mModel->GetParError(5);
262 mPhE.first = mModel->GetParameter(6);
263 mPhE.second = mModel->GetParError(6);
264}
265
266// Fit pulseshape with a continous funciton expressed by a product of two exp-functions
267void
268Pulse::FitContious(
269 TString fitName,
270 TString fitOptions,
271 int fitMin,
272 int fitMax
273 )
274{
275 int numPara = 6;
276 mFitMin = fitMin;
277 mFitMax = fitMax;
278 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, numPara );
279 mModel->SetNpx( mHisto->GetNbinsX());
280 // ======================================================================
281 // Calculate startvaues for fit
282 double tau = 30.0;
283 double bsl = 0; //MW der ersten zehn slices
284// double bsl2 = 0; //MW der ersten zehn slices
285 double gain = 10; //startvalue for gain
286 int counter = 0;
287
288 double stop = mHisto->GetMaximumBin(); //pos of max
289 double amplitude=0;
290
291 for (int i = (stop - 3) ; i < (stop + 3) && i <= fitMax && i >= fitMin ; i++)
292 {
293 counter++;
294 amplitude= mHisto->GetBinContent(i);
295 }
296 amplitude=amplitude/counter;
297
298 counter = 0;
299 for (int i = (stop - 15) ; i > fitMin+2 && i <= fitMax ; i--)
300 {
301 if (mHisto->GetBinContent(i) > 0)
302 {
303 break;
304 }
305 counter++;
306 bsl += mHisto->GetBinContent(i);
307 }
308 if (counter <= 0)
309 {
310 counter = 1;
311 }
312 bsl = bsl/counter;
313
314// counter = 0;
315// for (int i = stop + 100; i < stop + 200 && i <= fitMax && i >= fitMin ; i++)
316// {
317// if (mHisto->GetBinContent(i) > 0.5*amplitude )
318// {
319// break;
320// }
321// counter++;
322// bsl2 += mHisto->GetBinContent(i);
323// }
324// if (counter <= 0)
325// {
326// counter = 1;
327// }
328// bsl2 = bsl2/counter;
329
330 double start = mHisto->GetBinCenter(stop-10); //pos 10 slices before maximum
331 int phe = 0;
332 if (mOrder < 0)
333 {
334 phe = (amplitude-bsl)/gain;
335 }
336 else
337 {
338 phe = mOrder+1;
339 }
340 // ======================================================================
341
342
343 double para[] = {bsl, 2*amplitude+gain, start, tau*0.5, tau, phe};
344 mModel->SetParameters(para);
345 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2", "PhE");
346 mModel->SetLineColor(kBlue);
347 mModel->FixParameter(0, bsl);
348 mModel->FixParameter(1, 2*amplitude+gain);
349
350 // cout startvalues
351 if (mVerbosityLvl > 3)
352 {
353 cout << endl;
354 for (unsigned int i = 0; i < numPara; i++)
355 {
356 cout << "para[ " << i << "]: " << para[i] << endl;
357 }
358 }
359 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
360 mBsl.first = mModel->GetParameter(0);
361 mBsl.second = mModel->GetParError(0);
362 mHeight.first = mModel->GetParameter(1);
363 mHeight.second = mModel->GetParError(1);
364 mT0.first = mModel->GetParameter(2);
365 mT0.second = mModel->GetParError(2);
366 mTau1.first = mModel->GetParameter(3);
367 mTau1.second = mModel->GetParError(3);
368 mTau2.first = mModel->GetParameter(4);
369 mTau2.second = mModel->GetParError(4);
370 mPhE.first = mModel->GetParameter(5);
371 mPhE.second = mModel->GetParError(5);
372
373// cout << "test" << mBsl.first << "\t +/-" << mBsl.second << endl;
374}
375
376void
377Pulse::CalculateParameters()
378{
379 mIntegral.first = mModel->Integral(mFitMin, mFitMax);
380 //TODO setze konstantes glied auf 0 dann integrieren und fehler bestimmen
381 mIntegral.second = mModel->IntegralError(mFitMin, mFitMax);
382 mAmplitude.first = mModel->GetMaximum() - mBsl.first;
383// mAmplitude.second = mModel->GetMaximum() - mBsl.first;
384 mRisetime.first = mTimeResolution*CalculateHistEdgeRisetime();
385 mFitProb = mFitResultPtr->Prob();
386 mFitNCalls = mFitResultPtr->NCalls();
387 mFitNdf = mFitResultPtr->Ndf();
388 mChi2 = mFitResultPtr->Chi2();
389}
390
391double
392Pulse::CalculateHistEdgeRisetime() // in bins
393{
394 //rise-tim for pulse for pulse height from 10% to 90%
395 double riseTime = 0;
396
397 // Get Pulse Maximum value and bin
398 int maxBin = mHisto->GetMaximumBin();
399 double max = mHisto->GetBinContent(maxBin);
400
401 // values of bin with 10% and 90%
402 int tenPctBin = 0;
403 int ninetyPctBin = 0;
404
405 // flag points
406 bool finishedLow = false;
407 bool finishedHigh = false;
408
409 //2nd runvariable
410 int j = maxBin;
411
412 //loop over from 20 bins bevore max to maximum
413 for (int i = maxBin - 20; i<maxBin; i++){
414 j--;
415 if (j<0) break;
416
417 if (mHisto->GetBinContent(i) >= 0.1*max && !finishedLow)
418 {
419 tenPctBin = i;
420 finishedLow = true;
421 }
422
423
424 if (mHisto->GetBinContent(j) <= 0.9*max && !finishedHigh)
425 {
426 ninetyPctBin = j;
427 finishedHigh = true;
428 }
429
430 if (finishedLow && finishedHigh)
431 {
432 break;
433 }
434 }
435
436 riseTime = mHisto->GetBinLowEdge(ninetyPctBin) - mHisto->GetBinLowEdge(tenPctBin);
437 return riseTime;
438}
439
440
441
442
443// ===========================================================================
444// ACCESS
445// ===========================================================================
446
447TString Pulse::GetName(){ return mName;}
448double Pulse::GetBsl(){ return mBsl.first;}
449double Pulse::GetBslErr(){ return mBsl.second;}
450double Pulse::GetHeight(){ return mHeight.first;}
451double Pulse::GetHeightErr(){ return mHeight.second;}
452double Pulse::GetT0(){ return mT0.first;}
453double Pulse::GetT0Err(){ return mT0.second;}
454double Pulse::GetT1(){ return mT1.first;}
455double Pulse::GetT1Err(){ return mT1.second;}
456double Pulse::GetTau1(){ return mTau1.first;}
457double Pulse::GetTau1Err(){ return mTau1.second;}
458double Pulse::GetTau2(){ return mTau2.first;}
459double Pulse::GetTau2Err(){ return mTau2.second;}
460double Pulse::GetIntegral(){ return mIntegral.first;}
461double Pulse::GetIntegralErr(){ return mIntegral.second;}
462double Pulse::GetAmplitude(){ return mAmplitude.first;}
463double Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
464int Pulse::GetPhE(){ return mPhE.first;}
465int Pulse::GetPhEErr(){ return mPhE.second;}
466double Pulse::GetRiseTime(){ return mRisetime.first; }
467double Pulse::GetRiseTimeErr(){ return mRisetime.second;}
468int Pulse::GetOrder(){ return mOrder;}
469int Pulse::GetType(){return mType;}
470int Pulse::GetFitMin(){return mFitMin;}
471int Pulse::GetFitMax(){return mFitMax;}
472double Pulse::GetFitProb(){return mFitProb;}
473double Pulse::GetFitNCalls(){return mFitNCalls;}
474double Pulse::GetFitNdf(){return mFitNdf;}
475double Pulse::GetChi2(){return mChi2;}
476double Pulse::GetTimeResolution(){return mTimeResolution;}
477void Pulse::SetTimeResolution(double timeResolution){mTimeResolution = timeResolution;}
478void Pulse::SetVerbosity(int verbLvl){mVerbosityLvl = verbLvl;}
479
480
481// ===========================================================================
482// NON CLASS MATH FUNCTIONS
483// ===========================================================================
484
485//double
486int Heaviside(double val)
487{
488 if( val >= 0) return 1;
489 else return 0;
490}
491
492//Pulse::
493double shapeFunc(
494 double* t,
495 double* par
496 )
497{
498 double returnval= 0.0;
499
500 // name parameters
501 double bsl = par[0];
502 double height = par[1];
503 double start = par[2];
504// double rising = par[3];
505 double stop = par[3];
506 double tau1 = par[4];
507 double tau2 = par[5];
508 int phe = par[6];
509
510 // helper variables
511// double max = start + rising * tau1;
512 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
513// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
514 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
515 // calculate return value
516 returnval += Heaviside(t[0]-start)*e1;
517 returnval += Heaviside(t[0]-stop)*e2;
518 returnval *= phe*height;
519 returnval += bsl;
520// if (t[0] > start)
521// returnval += e1;
522// if (t[0] > stop)
523// returnval += e2;
524
525 return returnval;
526
527
528}
529
530
531
532/*
533[18:42:13] Dominik: können wir nachher skypen? haste n headset?
534[18:47:34] Dominik: hallo?
535[18:53:38] Dominik: so also ich hab das file glaubich gefunden
536[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
537[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
538[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
539[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
540[18:55:54] Dominik: und so wie sich ne normale abfallende e-fkt der null annähert, so nähert sich meine e1 dem wert 'height' an
541[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
542[18:57:20] Dominik: im Allgemeinen legt e2 schon los, bevor e1 auch nur in die Nähe von 'height' kommt und beginnt die Kurve nach unten zu drücken.
543[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
544[18:59:27] Dominik: aber ich vermute, wenn man die Funktion nimmt und wie inder 9ten klasse die ableitung macht um die Position und Höhe des Maximums zu bestimmen, dann geht 'height' als linearer Faktor in die Höhe ein und spielt für die Position des Maximums keine Rolle.
545Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
546
547p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
548
549bsl = p[0] # just the baseline
550 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
551 # it *would* be the pulse height only in case the
552 # recharging process would set in later...
553
554 start = p[2] # start of avalanche aka start of discharge of diode capacitance
555 stop = p[3] # stop of avalanche, or start of recharging of diode
556 tau1 = p[4] # time constant of discharging process
557 tau2 = p[5] # time constant of recharging process
558
559*/
560
561
562double shapeFunc2(
563 double* t,
564 double* par
565 )
566{
567 double returnval= 0.0;
568
569 // name parameters
570 double bsl = par[0];
571 double gain = par[1];
572 double t_0 = par[2];
573 double tau1 = par[3];
574 double tau2 = par[4];
575 int phe = par[5];
576// double bsl2 = par[6];
577
578 // helper variables
579 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
580 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
581
582 // calculate return value
583 returnval += bsl;
584 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0);
585
586 return returnval;
587}
Note: See TracBrowser for help on using the repository browser.