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

Last change on this file since 14762 was 14762, checked in by Jens Buss, 12 years ago
risetime function
  • Property svn:executable set to *
File size: 14.0 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
57
58Pulse::~Pulse()
59{
60 delete mModel;
61}
62
63// ===========================================================================
64// PUBLIC OPERATIONS
65// ===========================================================================
66
67void
68Pulse::Fit(
69 TString fitName,
70 TString fitOptions,
71 int type
72 )
73{
74 int fitMin = mHisto->GetXaxis()->GetBinLowEdge( mHisto->GetXaxis()->GetFirst() );
75 int fitMax = mHisto->GetXaxis()->GetBinUpEdge( mHisto->GetXaxis()->GetLast() );
76
77 Fit(
78 fitName,
79 fitOptions,
80 type,
81 fitMin,
82 fitMax
83 );
84}
85
86void
87Pulse::Fit(
88 TString fitName,
89 TString fitOptions,
90 int type,
91 int fitMin,
92 int fitMax
93 )
94{
95// mHisto->Sumw2();
96
97 fitOptions += "R";
98 if (type == 0) //use sectionwise function
99 FitSectionWise(
100 fitName,
101 fitOptions,
102 fitMin,
103 fitMax
104 );
105 else if (type == 1) //use continous function
106 FitContious(
107 fitName,
108 fitOptions,
109 fitMin,
110 fitMax
111 );
112 else{
113 cout << "not a correct type number --> fitting skipped!" << endl;
114 return;
115 }
116
117 CalculateParameters();
118 return;
119}
120
121// ===========================================================================
122// PRIVATE OPERATIONS
123// ===========================================================================
124
125void
126Pulse::InitMembers()
127{
128 mName = "pulse";
129 mHisto = NULL;
130 mModel = NULL;
131 mOptions = "S";
132 mBsl.first = 0;
133 mHeight.first = 0;
134 mT0.first = 0;
135 mT1.first = 0;
136 mTau1.first = 0;
137 mTau2.first = 0;
138 mIntegral.first = 0;
139 mAmplitude.first = 0;
140 mPhE.first = 0;
141 mBsl.second = 0;
142 mHeight.second = 0;
143 mT0.second = 0;
144 mT1.second = 0;
145 mTau1.second = 0;
146 mTau2.second = 0;
147 mIntegral.second = 0;
148 mAmplitude.second = 0;
149 mPhE.second = 0;
150 mRisetime.first = 0;
151 mRisetime.second = 0;
152 mOrder = 0;
153 mType = 0;
154 mFitProb = 0;
155 mFitNCalls = 0;
156 mFitNdf = 0;
157 mChi2 = 0;
158 mTimeResolution = 0.5;
159}
160
161// Fit pulseshape with a section wise defined funciton
162// expressed by a sum of section wise defined exp-functions
163//controled by heaviside function
164void
165Pulse::FitSectionWise(
166 TString fitName,
167 TString fitOptions,
168 int fitMin,
169 int fitMax
170 )
171{
172 mFitMin = fitMin;
173 mFitMax = fitMax;
174 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 7 );
175
176 // ======================================================================
177 // Calculate startvaues for fit
178 double tau = 30.0;
179 double bsl = 0.8; //shall be around 0
180
181 int first_bin = mHisto->GetXaxis()->GetFirst();
182
183 for (int i = 0; i < 30; i++){
184 bsl += mHisto->GetBinContent(first_bin+0);
185 }
186 bsl = bsl/30;
187
188 double stop = mHisto->GetMaximumBin(); //pos of max
189 double amplitude = mHisto->GetBinContent(stop);
190 double start = stop-10; //pos 10 slices before maximum
191 int phe = amplitude/10;
192 // ======================================================================
193
194 double para[] = {bsl, amplitude, start, stop, tau, tau, phe};
195
196
197
198 mModel->SetParameters(para);
199 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2", "PhE");
200 mModel->SetLineColor(kRed);
201
202 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
203 mBsl.first = mModel->GetParameter(0);
204 mBsl.second = mModel->GetParError(0);
205 mHeight.first = mModel->GetParameter(1);
206 mHeight.second = mModel->GetParError(1);
207 mT0.first = mModel->GetParameter(2);
208 mT0.second = mModel->GetParError(2);
209 mT1.first = mModel->GetParameter(3);
210 mT1.second = mModel->GetParameter(3);
211 mTau1.first = mModel->GetParameter(4);
212 mTau1.second = mModel->GetParError(4);
213 mTau2.first = mModel->GetParameter(5);
214 mTau2.second = mModel->GetParError(5);
215 mPhE.first = mModel->GetParameter(6);
216 mPhE.second = mModel->GetParError(6);
217}
218
219// Fit pulseshape with a continous funciton expressed by a product of two exp-functions
220void
221Pulse::FitContious(
222 TString fitName,
223 TString fitOptions,
224 int fitMin,
225 int fitMax
226 )
227{
228 mFitMin = fitMin;
229 mFitMax = fitMax;
230 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 6 );
231
232 // ======================================================================
233 // Calculate startvaues for fit
234 double tau = 30.0;
235 double bsl = 0; //MW der ersten zehn slices
236
237 int first_bin = mHisto->GetXaxis()->GetFirst();
238
239 for (int i = 0; i < 30; i++){
240 bsl += mHisto->GetBinContent(first_bin+0);
241 }
242 bsl = bsl/30;
243
244 double stop = mHisto->GetMaximumBin(); //pos of max
245 double amplitude= mHisto->GetBinContent(stop);
246 double start = stop-10; //pos 10 slices before maximum
247 int phe = amplitude/10;
248 // ======================================================================
249
250
251 double para[] = {bsl, amplitude, start, tau, tau, phe};
252 mModel->SetParameters(para);
253 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2", "PhE");
254 mModel->SetLineColor(kBlue);
255
256 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
257 mBsl.first = mModel->GetParameter(0);
258 mBsl.second = mModel->GetParError(0);
259 mHeight.first = mModel->GetParameter(1);
260 mHeight.second = mModel->GetParError(1);
261 mT0.first = mModel->GetParameter(2);
262 mT0.second = mModel->GetParError(2);
263 mTau1.first = mModel->GetParameter(3);
264 mTau1.second = mModel->GetParError(3);
265 mTau2.first = mModel->GetParameter(4);
266 mTau2.second = mModel->GetParError(4);
267 mPhE.first = mModel->GetParameter(5);
268 mPhE.second = mModel->GetParError(5);
269
270// cout << "test" << mBsl.first << "\t +/-" << mBsl.second << endl;
271}
272
273void
274Pulse::CalculateParameters()
275{
276 mIntegral.first = mModel->Integral(mFitMin, mFitMax);
277 mIntegral.second = mModel->IntegralError(mFitMin, mFitMax);
278 mAmplitude.first = mModel->GetMaximum() - mBsl.first;
279// mAmplitude.second = mModel->GetMaximum() - mBsl.first;
280 mRisetime.first = mTimeResolution*CalculateHistEdgeRisetime();
281 mFitProb = mFitResultPtr->Prob();
282 mFitNCalls = mFitResultPtr->NCalls();
283 mFitNdf = mFitResultPtr->Ndf();
284 mChi2 = mFitResultPtr->Chi2();
285}
286
287double
288Pulse::CalculateHistEdgeRisetime() // in bins
289{
290 //rise-tim for pulse for pulse height from 10% to 90%
291 double riseTime = 0;
292
293 // Get Pulse Maximum value and bin
294 int maxBin = mHisto->GetMaximumBin();
295 double max = mHisto->GetBinContent(maxBin);
296
297 // values of bin with 10% and 90%
298 int tenPctBin = 0;
299 int ninetyPctBin = 0;
300
301 // flag points
302 bool finishedLow = false;
303 bool finishedHigh = false;
304
305 //2nd runvariable
306 int j = maxBin;
307
308 //loop over from 20 bins bevore max to maximum
309 for (int i = maxBin - 20; i<maxBin; i++){
310 j--;
311 if (j<0) break;
312
313 if (mHisto->GetBinContent(i) >= 0.1*max && !finishedLow)
314 {
315 tenPctBin = i;
316 finishedLow = true;
317 }
318
319
320 if (mHisto->GetBinContent(j) <= 0.9*max && !finishedHigh)
321 {
322 ninetyPctBin = j;
323 finishedHigh = true;
324 }
325
326 if (finishedLow && finishedHigh)
327 {
328 break;
329 }
330 }
331
332 riseTime = mHisto->GetBinLowEdge(ninetyPctBin) - mHisto->GetBinLowEdge(tenPctBin);
333 return riseTime;
334}
335
336
337
338
339// ===========================================================================
340// ACCESS
341// ===========================================================================
342
343TString Pulse::GetName(){ return mName;}
344double Pulse::GetBsl(){ return mBsl.first;}
345double Pulse::GetBslErr(){ return mBsl.second;}
346double Pulse::GetHeight(){ return mHeight.first;}
347double Pulse::GetHeightErr(){ return mHeight.second;}
348double Pulse::GetT0(){ return mT0.first;}
349double Pulse::GetT0Err(){ return mT0.second;}
350double Pulse::GetT1(){ return mT1.first;}
351double Pulse::GetT1Err(){ return mT1.second;}
352double Pulse::GetTau1(){ return mTau1.first;}
353double Pulse::GetTau1Err(){ return mTau1.second;}
354double Pulse::GetTau2(){ return mTau2.first;}
355double Pulse::GetTau2Err(){ return mTau2.second;}
356double Pulse::GetIntegral(){ return mIntegral.first;}
357double Pulse::GetIntegralErr(){ return mIntegral.second;}
358double Pulse::GetAmplitude(){ return mAmplitude.first;}
359double Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
360int Pulse::GetPhE(){ return mPhE.first;}
361int Pulse::GetPhEErr(){ return mPhE.second;}
362double Pulse::GetRiseTime(){ return mRisetime.first; }
363double Pulse::GetRiseTimeErr(){ return mRisetime.second;}
364int Pulse::GetOrder(){ return mOrder;}
365int Pulse::GetType(){return mType;}
366int Pulse::GetFitMin(){return mFitMin;}
367int Pulse::GetFitMax(){return mFitMax;}
368double Pulse::GetFitProb(){return mFitProb;}
369double Pulse::GetFitNCalls(){return mFitNCalls;}
370double Pulse::GetFitNdf(){return mFitNdf;}
371double Pulse::GetChi2(){return mChi2;}
372double Pulse::GetTimeResolution(){return mTimeResolution;}
373void Pulse::SetTimeResolution(double timeResolution){mTimeResolution = timeResolution;}
374
375
376// ===========================================================================
377// NON CLASS MATH FUNCTIONS
378// ===========================================================================
379
380//double
381int Heaviside(double val)
382{
383 if( val >= 0) return 1;
384 else return 0;
385}
386
387//Pulse::
388double shapeFunc(
389 double* t,
390 double* par
391 )
392{
393 double returnval= 0.0;
394
395 // name parameters
396 double bsl = par[0];
397 double height = par[1];
398 double start = par[2];
399// double rising = par[3];
400 double stop = par[3];
401 double tau1 = par[4];
402 double tau2 = par[5];
403 int phe = par[6];
404
405 // helper variables
406// double max = start + rising * tau1;
407 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
408// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
409 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
410 // calculate return value
411 returnval += Heaviside(t[0]-start)*e1;
412 returnval += Heaviside(t[0]-stop)*e2;
413 returnval *= phe*height;
414 returnval += bsl;
415// if (t[0] > start)
416// returnval += e1;
417// if (t[0] > stop)
418// returnval += e2;
419
420 return returnval;
421
422
423}
424
425
426
427/*
428[18:42:13] Dominik: können wir nachher skypen? haste n headset?
429[18:47:34] Dominik: hallo?
430[18:53:38] Dominik: so also ich hab das file glaubich gefunden
431[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
432[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
433[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
434[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
435[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
436[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
437[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.
438[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
439[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.
440Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
441
442p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
443
444bsl = p[0] # just the baseline
445 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
446 # it *would* be the pulse height only in case the
447 # recharging process would set in later...
448
449 start = p[2] # start of avalanche aka start of discharge of diode capacitance
450 stop = p[3] # stop of avalanche, or start of recharging of diode
451 tau1 = p[4] # time constant of discharging process
452 tau2 = p[5] # time constant of recharging process
453
454*/
455
456
457double shapeFunc2(
458 double* t,
459 double* par
460 )
461{
462 double returnval= 0.0;
463
464 // name parameters
465 double bsl = par[0];
466 double gain = par[1];
467 double t_0 = par[2];
468 double tau1 = par[3];
469 double tau2 = par[4];
470 int phe = par[5];
471
472 // helper variables
473 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
474 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
475
476 // calculate return value
477 returnval += bsl;
478 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0);
479
480 return returnval;
481}
Note: See TracBrowser for help on using the repository browser.