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

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