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

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