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

Last change on this file since 14749 was 14748, checked in by Jens Buss, 13 years ago
comments and changed start values for fit
  • Property svn:executable set to *
File size: 12.3 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()->GetFirst();
75 int fitMax = 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 mAmplitude.first = mModel->GetMaximum() - mBsl.first;
275 mFitProb = mFitResultPtr->Prob();
276 mFitNCalls = mFitResultPtr->NCalls();
277 mFitNdf = mFitResultPtr->Ndf();
278 mChi2 = mFitResultPtr->Chi2();
279}
280
281// ===========================================================================
282// ACCESS
283// ===========================================================================
284
285TString Pulse::GetName(){ return mName;}
286double Pulse::GetBsl(){ return mBsl.first;}
287double Pulse::GetHeight(){ return mHeight.first;}
288double Pulse::GetT0(){ return mT0.first;}
289double Pulse::GetT1(){ return mT1.first;}
290double Pulse::GetTau1(){ return mTau1.first;}
291double Pulse::GetTau2(){ return mTau2.first;}
292double Pulse::GetIntegral(){ return mIntegral.first;}
293double Pulse::GetAmplitude(){ return mAmplitude.first;}
294int Pulse::GetPhE(){ return mPhE.first;}
295double Pulse::GetBslErr(){ return mBsl.second;}
296double Pulse::GetHeightErr(){ return mHeight.second;}
297double Pulse::GetT0Err(){ return mT0.second;}
298double Pulse::GetT1Err(){ return mT1.second;}
299double Pulse::GetTau1Err(){ return mTau1.second;}
300double Pulse::GetTau2Err(){ return mTau2.second;}
301double Pulse::GetIntegralErr(){ return mIntegral.second;}
302double Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
303int Pulse::GetPhEErr(){ return mPhE.second;}
304int Pulse::GetOrder(){ return mOrder;}
305int Pulse::GetType(){return mType;}
306int Pulse::GetFitMin(){return mFitMin;}
307int Pulse::GetFitMax(){return mFitMax;}
308double Pulse::GetFitProb(){return mFitProb;}
309double Pulse::GetFitNCalls(){return mFitNCalls;}
310double Pulse::GetFitNdf(){return mFitNdf;}
311double Pulse::GetChi2(){return mChi2;}
312
313
314// ===========================================================================
315// NON CLASS MATH FUNCTIONS
316// ===========================================================================
317
318//double
319int Heaviside(double val)
320{
321 if( val >= 0) return 1;
322 else return 0;
323}
324
325//Pulse::
326double shapeFunc(
327 double* t,
328 double* par
329 )
330{
331 double returnval= 0.0;
332
333 // name parameters
334 double bsl = par[0];
335 double height = par[1];
336 double start = par[2];
337// double rising = par[3];
338 double stop = par[3];
339 double tau1 = par[4];
340 double tau2 = par[5];
341 int phe = par[6];
342
343 // helper variables
344// double max = start + rising * tau1;
345 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
346// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
347 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
348 // calculate return value
349 returnval += Heaviside(t[0]-start)*e1;
350 returnval += Heaviside(t[0]-stop)*e2;
351 returnval *= phe*height;
352 returnval += bsl;
353// if (t[0] > start)
354// returnval += e1;
355// if (t[0] > stop)
356// returnval += e2;
357
358 return returnval;
359
360
361}
362
363
364
365/*
366[18:42:13] Dominik: können wir nachher skypen? haste n headset?
367[18:47:34] Dominik: hallo?
368[18:53:38] Dominik: so also ich hab das file glaubich gefunden
369[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
370[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
371[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
372[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
373[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
374[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
375[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.
376[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
377[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.
378Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
379
380p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
381
382bsl = p[0] # just the baseline
383 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
384 # it *would* be the pulse height only in case the
385 # recharging process would set in later...
386
387 start = p[2] # start of avalanche aka start of discharge of diode capacitance
388 stop = p[3] # stop of avalanche, or start of recharging of diode
389 tau1 = p[4] # time constant of discharging process
390 tau2 = p[5] # time constant of recharging process
391
392*/
393
394
395double shapeFunc2(
396 double* t,
397 double* par
398 )
399{
400 double returnval= 0.0;
401
402 // name parameters
403 double bsl = par[0];
404 double gain = par[1];
405 double t_0 = par[2];
406 double tau1 = par[3];
407 double tau2 = par[4];
408 int phe = par[5];
409
410 // helper variables
411 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
412 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
413
414 // calculate return value
415 returnval += bsl;
416 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0);
417
418 return returnval;
419}
Note: See TracBrowser for help on using the repository browser.