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

Last change on this file since 14476 was 14476, checked in by Jens Buss, 12 years ago
calculated model attributes
  • Property svn:executable set to *
File size: 9.2 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
63void
64Pulse::InitMembers()
65{
66 mName = "pulse";
67 mHisto = NULL;
68 mModel = NULL;
69 mOptions = "S";
70 mBsl = 0;
71 mHeight = 0;
72 mT0 = 0;
73 mT1 = 0;
74 mTau1 = 0;
75 mTau2 = 0;
76 mIntegral = 0;
77 mAmplitude = 0;
78 mPhE = 0;
79 mType = 0;
80 mFitProb = 0;
81 mFitNCalls = 0;
82 mFitNdf = 0;
83 mChi2 = 0;
84}
85
86void
87Pulse::Fit(
88 TString fitName,
89 TString fitOptions,
90 int type
91 )
92{
93 int fitMin = mHisto->GetXaxis()->GetFirst();
94 int fitMax = mHisto->GetXaxis()->GetLast();
95
96 Fit(
97 fitName,
98 fitOptions,
99 type,
100 fitMin,
101 fitMax
102 );
103}
104
105void
106Pulse::Fit(
107 TString fitName,
108 TString fitOptions,
109 int type,
110 int fitMin,
111 int fitMax
112 )
113{
114 fitOptions += "R";
115 if (type == 0)
116 FitSectionWise(
117 "Model0",
118 fitOptions,
119 fitMin,
120 fitMax
121 );
122 else if (type == 1)
123 FitContious(
124 "Model1",
125 fitOptions,
126 fitMin,
127 fitMax
128 );
129 else{
130 cout << "not a correct type number --> fitting skipped!" << endl;
131 return;
132 }
133
134 CalculateParameters();
135 return;
136}
137
138void
139Pulse::FitSectionWise(
140 TString fitName,
141 TString fitOptions,
142 int fitMin,
143 int fitMax
144 )
145{
146 mFitMin = fitMin;
147 mFitMax = fitMax;
148 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 6 );
149
150 // ======================================================================
151 // Calculate startvaues for fit
152 double tau = 30.0;
153 double bsl = 0; //MW der ersten zehn slices
154
155 int first_bin = mHisto->GetXaxis()->GetFirst();
156
157 for (int i = 0; i < 30; i++){
158 bsl += mHisto->GetBinContent(first_bin+0);
159 }
160 bsl = bsl/30;
161
162 double stop = mHisto->GetMaximumBin(); //pos of max
163 double height = mHisto->GetBinContent(stop);
164 double start = stop-10; //pos 10 slices before maximum
165 // ======================================================================
166
167 double para[] = {bsl, height, start, stop, tau, tau};
168
169
170
171 mModel->SetParameters(para);
172 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2");
173 mModel->SetLineColor(kRed);
174
175 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
176
177 mBsl = mModel->GetParameter(0);
178 mHeight = mModel->GetParameter(1);
179 mT0 = mModel->GetParameter(2);
180 mT1 = mModel->GetParameter(3);
181 mTau1 = mModel->GetParameter(4);
182 mTau2 = mModel->GetParameter(5);
183}
184
185void
186Pulse::FitContious(
187 TString fitName,
188 TString fitOptions,
189 int fitMin,
190 int fitMax
191 )
192{
193 mFitMin = fitMin;
194 mFitMax = fitMax;
195 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 5 );
196
197 // ======================================================================
198 // Calculate startvaues for fit
199 double tau = 30.0;
200 double bsl = 0; //MW der ersten zehn slices
201
202 int first_bin = mHisto->GetXaxis()->GetFirst();
203
204 for (int i = 0; i < 30; i++){
205 bsl += mHisto->GetBinContent(first_bin+0);
206 }
207 bsl = bsl/30;
208
209 double stop = mHisto->GetMaximumBin(); //pos of max
210 double amplitude= mHisto->GetBinContent(stop);
211 double start = stop-10; //pos 10 slices before maximum
212 // ======================================================================
213
214
215 double para[] = {bsl, amplitude, start, tau, tau};
216 mModel->SetParameters(para);
217 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2");
218 mModel->SetLineColor(kBlue);
219
220 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
221
222 mBsl = mModel->GetParameter(0);
223 mHeight = mModel->GetParameter(1);
224 mT0 = mModel->GetParameter(2);
225 mTau1 = mModel->GetParameter(3);
226 mTau2 = mModel->GetParameter(4);
227}
228
229double Pulse::GetBsl(){ return mBsl;}
230double Pulse::GetHeight(){ return mHeight;}
231double Pulse::GetAvalancheStart(){ return mT0;}
232double Pulse::GetAvalancheEnd(){ return mT1;}
233double Pulse::GetTimeConstRising(){ return mTau1;}
234double Pulse::GetTimeConstFalling(){ return mTau2;}
235double Pulse::GetIntegral(){ return mIntegral;}
236double Pulse::GetAmplitude(){ return mAmplitude;}
237int Pulse::GetPE(){ return mPhE;}
238
239void
240Pulse::CalculateParameters()
241{
242 mIntegral = mModel->Integral(mFitMin, mFitMax);
243 mAmplitude = mModel->GetMaximum() - mBsl;
244 mFitProb = mFitResultPtr->Prob();
245 mFitNCalls = mFitResultPtr->NCalls();
246 mFitNdf = mFitResultPtr->Ndf();
247 mChi2 = mFitResultPtr->Chi2();
248}
249
250
251//double
252int Heaviside(double val)
253{
254 if( val >= 0) return 1;
255 else return 0;
256}
257
258//Pulse::
259double shapeFunc(
260 double* t,
261 double* par
262 )
263{
264 double returnval= 0.0;
265
266 // name parameters
267 double bsl = par[0];
268 double height = par[1];
269 double start = par[2];
270// double rising = par[3];
271 double tau1 = par[4];
272 double tau2 = par[5];
273 double stop = par[3];
274
275 // helper variables
276// double max = start + rising * tau1;
277 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
278// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
279 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
280 // calculate return value
281 returnval += Heaviside(t[0]-start)*e1;
282 returnval += Heaviside(t[0]-stop)*e2;
283 returnval *= height;
284 returnval += bsl;
285// if (t[0] > start)
286// returnval += e1;
287// if (t[0] > stop)
288// returnval += e2;
289
290 return returnval;
291
292
293}
294
295
296
297/*
298[18:42:13] Dominik: können wir nachher skypen? haste n headset?
299[18:47:34] Dominik: hallo?
300[18:53:38] Dominik: so also ich hab das file glaubich gefunden
301[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
302[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
303[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
304[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
305[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
306[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
307[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.
308[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
309[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.
310Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
311
312p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
313
314bsl = p[0] # just the baseline
315 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
316 # it *would* be the pulse height only in case the
317 # recharging process would set in later...
318
319 start = p[2] # start of avalanche aka start of discharge of diode capacitance
320 stop = p[3] # stop of avalanche, or start of recharging of diode
321 tau1 = p[4] # time constant of discharging process
322 tau2 = p[5] # time constant of recharging process
323
324*/
325
326
327double shapeFunc2(
328 double* t,
329 double* par
330 )
331{
332 double returnval= 0.0;
333
334 // name parameters
335 double bsl = par[0];
336 double gain = par[1];
337 double t_0 = par[2];
338 double tau1 = par[3];
339 double tau2 = par[4];
340
341 // helper variables
342 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
343 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
344
345 // calculate return value
346 returnval += bsl;
347 returnval += gain*e1*e2*Heaviside(t[0]-t_0);
348
349 return returnval;
350}
Note: See TracBrowser for help on using the repository browser.