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

Last change on this file since 14807 was 14780, checked in by Jens Buss, 12 years ago
TODO COMMENT
  • Property svn:executable set to *
File size: 14.1 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 //TODO setze konstantes glied auf 0 dann integrieren und fehler bestimmen
278 mIntegral.second = mModel->IntegralError(mFitMin, mFitMax);
279 mAmplitude.first = mModel->GetMaximum() - mBsl.first;
280// mAmplitude.second = mModel->GetMaximum() - mBsl.first;
281 mRisetime.first = mTimeResolution*CalculateHistEdgeRisetime();
282 mFitProb = mFitResultPtr->Prob();
283 mFitNCalls = mFitResultPtr->NCalls();
284 mFitNdf = mFitResultPtr->Ndf();
285 mChi2 = mFitResultPtr->Chi2();
286}
287
288double
289Pulse::CalculateHistEdgeRisetime() // in bins
290{
291 //rise-tim for pulse for pulse height from 10% to 90%
292 double riseTime = 0;
293
294 // Get Pulse Maximum value and bin
295 int maxBin = mHisto->GetMaximumBin();
296 double max = mHisto->GetBinContent(maxBin);
297
298 // values of bin with 10% and 90%
299 int tenPctBin = 0;
300 int ninetyPctBin = 0;
301
302 // flag points
303 bool finishedLow = false;
304 bool finishedHigh = false;
305
306 //2nd runvariable
307 int j = maxBin;
308
309 //loop over from 20 bins bevore max to maximum
310 for (int i = maxBin - 20; i<maxBin; i++){
311 j--;
312 if (j<0) break;
313
314 if (mHisto->GetBinContent(i) >= 0.1*max && !finishedLow)
315 {
316 tenPctBin = i;
317 finishedLow = true;
318 }
319
320
321 if (mHisto->GetBinContent(j) <= 0.9*max && !finishedHigh)
322 {
323 ninetyPctBin = j;
324 finishedHigh = true;
325 }
326
327 if (finishedLow && finishedHigh)
328 {
329 break;
330 }
331 }
332
333 riseTime = mHisto->GetBinLowEdge(ninetyPctBin) - mHisto->GetBinLowEdge(tenPctBin);
334 return riseTime;
335}
336
337
338
339
340// ===========================================================================
341// ACCESS
342// ===========================================================================
343
344TString Pulse::GetName(){ return mName;}
345double Pulse::GetBsl(){ return mBsl.first;}
346double Pulse::GetBslErr(){ return mBsl.second;}
347double Pulse::GetHeight(){ return mHeight.first;}
348double Pulse::GetHeightErr(){ return mHeight.second;}
349double Pulse::GetT0(){ return mT0.first;}
350double Pulse::GetT0Err(){ return mT0.second;}
351double Pulse::GetT1(){ return mT1.first;}
352double Pulse::GetT1Err(){ return mT1.second;}
353double Pulse::GetTau1(){ return mTau1.first;}
354double Pulse::GetTau1Err(){ return mTau1.second;}
355double Pulse::GetTau2(){ return mTau2.first;}
356double Pulse::GetTau2Err(){ return mTau2.second;}
357double Pulse::GetIntegral(){ return mIntegral.first;}
358double Pulse::GetIntegralErr(){ return mIntegral.second;}
359double Pulse::GetAmplitude(){ return mAmplitude.first;}
360double Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
361int Pulse::GetPhE(){ return mPhE.first;}
362int Pulse::GetPhEErr(){ return mPhE.second;}
363double Pulse::GetRiseTime(){ return mRisetime.first; }
364double Pulse::GetRiseTimeErr(){ return mRisetime.second;}
365int Pulse::GetOrder(){ return mOrder;}
366int Pulse::GetType(){return mType;}
367int Pulse::GetFitMin(){return mFitMin;}
368int Pulse::GetFitMax(){return mFitMax;}
369double Pulse::GetFitProb(){return mFitProb;}
370double Pulse::GetFitNCalls(){return mFitNCalls;}
371double Pulse::GetFitNdf(){return mFitNdf;}
372double Pulse::GetChi2(){return mChi2;}
373double Pulse::GetTimeResolution(){return mTimeResolution;}
374void Pulse::SetTimeResolution(double timeResolution){mTimeResolution = timeResolution;}
375
376
377// ===========================================================================
378// NON CLASS MATH FUNCTIONS
379// ===========================================================================
380
381//double
382int Heaviside(double val)
383{
384 if( val >= 0) return 1;
385 else return 0;
386}
387
388//Pulse::
389double shapeFunc(
390 double* t,
391 double* par
392 )
393{
394 double returnval= 0.0;
395
396 // name parameters
397 double bsl = par[0];
398 double height = par[1];
399 double start = par[2];
400// double rising = par[3];
401 double stop = par[3];
402 double tau1 = par[4];
403 double tau2 = par[5];
404 int phe = par[6];
405
406 // helper variables
407// double max = start + rising * tau1;
408 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
409// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
410 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
411 // calculate return value
412 returnval += Heaviside(t[0]-start)*e1;
413 returnval += Heaviside(t[0]-stop)*e2;
414 returnval *= phe*height;
415 returnval += bsl;
416// if (t[0] > start)
417// returnval += e1;
418// if (t[0] > stop)
419// returnval += e2;
420
421 return returnval;
422
423
424}
425
426
427
428/*
429[18:42:13] Dominik: können wir nachher skypen? haste n headset?
430[18:47:34] Dominik: hallo?
431[18:53:38] Dominik: so also ich hab das file glaubich gefunden
432[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
433[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
434[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
435[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
436[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
437[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
438[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.
439[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
440[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.
441Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
442
443p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
444
445bsl = p[0] # just the baseline
446 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
447 # it *would* be the pulse height only in case the
448 # recharging process would set in later...
449
450 start = p[2] # start of avalanche aka start of discharge of diode capacitance
451 stop = p[3] # stop of avalanche, or start of recharging of diode
452 tau1 = p[4] # time constant of discharging process
453 tau2 = p[5] # time constant of recharging process
454
455*/
456
457
458double shapeFunc2(
459 double* t,
460 double* par
461 )
462{
463 double returnval= 0.0;
464
465 // name parameters
466 double bsl = par[0];
467 double gain = par[1];
468 double t_0 = par[2];
469 double tau1 = par[3];
470 double tau2 = par[4];
471 int phe = par[5];
472
473 // helper variables
474 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
475 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
476
477 // calculate return value
478 returnval += bsl;
479 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0);
480
481 return returnval;
482}
Note: See TracBrowser for help on using the repository browser.