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

Last change on this file since 14760 was 14755, checked in by Jens Buss, 13 years ago
calc integral error, calc fit boarders correctly
  • Property svn:executable set to *
File size: 12.5 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 mOrder = 0;
151 mType = 0;
152 mFitProb = 0;
153 mFitNCalls = 0;
154 mFitNdf = 0;
155 mChi2 = 0;
156}
157
158// Fit pulseshape with a section wise defined funciton
159// expressed by a sum of section wise defined exp-functions
160//controled by heaviside function
161void
162Pulse::FitSectionWise(
163 TString fitName,
164 TString fitOptions,
165 int fitMin,
166 int fitMax
167 )
168{
169 mFitMin = fitMin;
170 mFitMax = fitMax;
171 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 7 );
172
173 // ======================================================================
174 // Calculate startvaues for fit
175 double tau = 30.0;
176 double bsl = 0.8; //shall be around 0
177
178 int first_bin = mHisto->GetXaxis()->GetFirst();
179
180 for (int i = 0; i < 30; i++){
181 bsl += mHisto->GetBinContent(first_bin+0);
182 }
183 bsl = bsl/30;
184
185 double stop = mHisto->GetMaximumBin(); //pos of max
186 double amplitude = mHisto->GetBinContent(stop);
187 double start = stop-10; //pos 10 slices before maximum
188 int phe = amplitude/10;
189 // ======================================================================
190
191 double para[] = {bsl, amplitude, start, stop, tau, tau, phe};
192
193
194
195 mModel->SetParameters(para);
196 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2", "PhE");
197 mModel->SetLineColor(kRed);
198
199 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
200 mBsl.first = mModel->GetParameter(0);
201 mBsl.second = mModel->GetParError(0);
202 mHeight.first = mModel->GetParameter(1);
203 mHeight.second = mModel->GetParError(1);
204 mT0.first = mModel->GetParameter(2);
205 mT0.second = mModel->GetParError(2);
206 mT1.first = mModel->GetParameter(3);
207 mT1.second = mModel->GetParameter(3);
208 mTau1.first = mModel->GetParameter(4);
209 mTau1.second = mModel->GetParError(4);
210 mTau2.first = mModel->GetParameter(5);
211 mTau2.second = mModel->GetParError(5);
212 mPhE.first = mModel->GetParameter(6);
213 mPhE.second = mModel->GetParError(6);
214}
215
216// Fit pulseshape with a continous funciton expressed by a product of two exp-functions
217void
218Pulse::FitContious(
219 TString fitName,
220 TString fitOptions,
221 int fitMin,
222 int fitMax
223 )
224{
225 mFitMin = fitMin;
226 mFitMax = fitMax;
227 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 6 );
228
229 // ======================================================================
230 // Calculate startvaues for fit
231 double tau = 30.0;
232 double bsl = 0; //MW der ersten zehn slices
233
234 int first_bin = mHisto->GetXaxis()->GetFirst();
235
236 for (int i = 0; i < 30; i++){
237 bsl += mHisto->GetBinContent(first_bin+0);
238 }
239 bsl = bsl/30;
240
241 double stop = mHisto->GetMaximumBin(); //pos of max
242 double amplitude= mHisto->GetBinContent(stop);
243 double start = stop-10; //pos 10 slices before maximum
244 int phe = amplitude/10;
245 // ======================================================================
246
247
248 double para[] = {bsl, amplitude, start, tau, tau, phe};
249 mModel->SetParameters(para);
250 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2", "PhE");
251 mModel->SetLineColor(kBlue);
252
253 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
254 mBsl.first = mModel->GetParameter(0);
255 mBsl.second = mModel->GetParError(0);
256 mHeight.first = mModel->GetParameter(1);
257 mHeight.second = mModel->GetParError(1);
258 mT0.first = mModel->GetParameter(2);
259 mT0.second = mModel->GetParError(2);
260 mTau1.first = mModel->GetParameter(3);
261 mTau1.second = mModel->GetParError(3);
262 mTau2.first = mModel->GetParameter(4);
263 mTau2.second = mModel->GetParError(4);
264 mPhE.first = mModel->GetParameter(5);
265 mPhE.second = mModel->GetParError(5);
266
267// cout << "test" << mBsl.first << "\t +/-" << mBsl.second << endl;
268}
269
270void
271Pulse::CalculateParameters()
272{
273 mIntegral.first = mModel->Integral(mFitMin, mFitMax);
274 mIntegral.second = mModel->IntegralError(mFitMin, mFitMax);
275 mAmplitude.first = mModel->GetMaximum() - mBsl.first;
276// mAmplitude.second = mModel->GetMaximum() - mBsl.first;
277 mFitProb = mFitResultPtr->Prob();
278 mFitNCalls = mFitResultPtr->NCalls();
279 mFitNdf = mFitResultPtr->Ndf();
280 mChi2 = mFitResultPtr->Chi2();
281}
282
283// ===========================================================================
284// ACCESS
285// ===========================================================================
286
287TString Pulse::GetName(){ return mName;}
288double Pulse::GetBsl(){ return mBsl.first;}
289double Pulse::GetHeight(){ return mHeight.first;}
290double Pulse::GetT0(){ return mT0.first;}
291double Pulse::GetT1(){ return mT1.first;}
292double Pulse::GetTau1(){ return mTau1.first;}
293double Pulse::GetTau2(){ return mTau2.first;}
294double Pulse::GetIntegral(){ return mIntegral.first;}
295double Pulse::GetAmplitude(){ return mAmplitude.first;}
296int Pulse::GetPhE(){ return mPhE.first;}
297double Pulse::GetBslErr(){ return mBsl.second;}
298double Pulse::GetHeightErr(){ return mHeight.second;}
299double Pulse::GetT0Err(){ return mT0.second;}
300double Pulse::GetT1Err(){ return mT1.second;}
301double Pulse::GetTau1Err(){ return mTau1.second;}
302double Pulse::GetTau2Err(){ return mTau2.second;}
303double Pulse::GetIntegralErr(){ return mIntegral.second;}
304double Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
305int Pulse::GetPhEErr(){ return mPhE.second;}
306int Pulse::GetOrder(){ return mOrder;}
307int Pulse::GetType(){return mType;}
308int Pulse::GetFitMin(){return mFitMin;}
309int Pulse::GetFitMax(){return mFitMax;}
310double Pulse::GetFitProb(){return mFitProb;}
311double Pulse::GetFitNCalls(){return mFitNCalls;}
312double Pulse::GetFitNdf(){return mFitNdf;}
313double Pulse::GetChi2(){return mChi2;}
314
315
316// ===========================================================================
317// NON CLASS MATH FUNCTIONS
318// ===========================================================================
319
320//double
321int Heaviside(double val)
322{
323 if( val >= 0) return 1;
324 else return 0;
325}
326
327//Pulse::
328double shapeFunc(
329 double* t,
330 double* par
331 )
332{
333 double returnval= 0.0;
334
335 // name parameters
336 double bsl = par[0];
337 double height = par[1];
338 double start = par[2];
339// double rising = par[3];
340 double stop = par[3];
341 double tau1 = par[4];
342 double tau2 = par[5];
343 int phe = par[6];
344
345 // helper variables
346// double max = start + rising * tau1;
347 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
348// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
349 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
350 // calculate return value
351 returnval += Heaviside(t[0]-start)*e1;
352 returnval += Heaviside(t[0]-stop)*e2;
353 returnval *= phe*height;
354 returnval += bsl;
355// if (t[0] > start)
356// returnval += e1;
357// if (t[0] > stop)
358// returnval += e2;
359
360 return returnval;
361
362
363}
364
365
366
367/*
368[18:42:13] Dominik: können wir nachher skypen? haste n headset?
369[18:47:34] Dominik: hallo?
370[18:53:38] Dominik: so also ich hab das file glaubich gefunden
371[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
372[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
373[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
374[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
375[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
376[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
377[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.
378[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
379[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.
380Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
381
382p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
383
384bsl = p[0] # just the baseline
385 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
386 # it *would* be the pulse height only in case the
387 # recharging process would set in later...
388
389 start = p[2] # start of avalanche aka start of discharge of diode capacitance
390 stop = p[3] # stop of avalanche, or start of recharging of diode
391 tau1 = p[4] # time constant of discharging process
392 tau2 = p[5] # time constant of recharging process
393
394*/
395
396
397double shapeFunc2(
398 double* t,
399 double* par
400 )
401{
402 double returnval= 0.0;
403
404 // name parameters
405 double bsl = par[0];
406 double gain = par[1];
407 double t_0 = par[2];
408 double tau1 = par[3];
409 double tau2 = par[4];
410 int phe = par[5];
411
412 // helper variables
413 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
414 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
415
416 // calculate return value
417 returnval += bsl;
418 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0);
419
420 return returnval;
421}
Note: See TracBrowser for help on using the repository browser.