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

Last change on this file since 14818 was 14814, checked in by Jens Buss, 12 years ago
new construktor is pulse order, new calculation of phe
  • Property svn:executable set to *
File size: 14.7 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
57Pulse::Pulse(TString name, TH1* histo, TString options, int type, int order)
58{
59 InitMembers();
60 mName = name;
61 mHisto = histo;
62 mOptions += options;
63 mType = type;
64 mOrder = order;
65 Fit(mName, mOptions, mType);
66}
67
68
69Pulse::~Pulse()
70{
71 delete mModel;
72}
73
74// ===========================================================================
75// PUBLIC OPERATIONS
76// ===========================================================================
77
78void
79Pulse::Fit(
80 TString fitName,
81 TString fitOptions,
82 int type
83 )
84{
85 int fitMin = mHisto->GetXaxis()->GetBinLowEdge( mHisto->GetXaxis()->GetFirst() );
86 int fitMax = mHisto->GetXaxis()->GetBinUpEdge( mHisto->GetXaxis()->GetLast() );
87
88 Fit(
89 fitName,
90 fitOptions,
91 type,
92 fitMin,
93 fitMax
94 );
95}
96
97void
98Pulse::Fit(
99 TString fitName,
100 TString fitOptions,
101 int type,
102 int fitMin,
103 int fitMax
104 )
105{
106// mHisto->Sumw2();
107
108 fitOptions += "R";
109 if (type == 0) //use sectionwise function
110 FitSectionWise(
111 fitName,
112 fitOptions,
113 fitMin,
114 fitMax
115 );
116 else if (type == 1) //use continous function
117 FitContious(
118 fitName,
119 fitOptions,
120 fitMin,
121 fitMax
122 );
123 else{
124 cout << "not a correct type number --> fitting skipped!" << endl;
125 return;
126 }
127
128 CalculateParameters();
129 return;
130}
131
132// ===========================================================================
133// PRIVATE OPERATIONS
134// ===========================================================================
135
136void
137Pulse::InitMembers()
138{
139 mName = "pulse";
140 mHisto = NULL;
141 mModel = NULL;
142 mOptions = "S";
143 mBsl.first = 0;
144 mHeight.first = 0;
145 mT0.first = 0;
146 mT1.first = 0;
147 mTau1.first = 0;
148 mTau2.first = 0;
149 mIntegral.first = 0;
150 mAmplitude.first = 0;
151 mPhE.first = 0;
152 mBsl.second = 0;
153 mHeight.second = 0;
154 mT0.second = 0;
155 mT1.second = 0;
156 mTau1.second = 0;
157 mTau2.second = 0;
158 mIntegral.second = 0;
159 mAmplitude.second = 0;
160 mPhE.second = 0;
161 mRisetime.first = 0;
162 mRisetime.second = 0;
163 mOrder = 0;
164 mType = 0;
165 mFitProb = 0;
166 mFitNCalls = 0;
167 mFitNdf = 0;
168 mChi2 = 0;
169 mTimeResolution = 0.5;
170}
171
172// Fit pulseshape with a section wise defined funciton
173// expressed by a sum of section wise defined exp-functions
174//controled by heaviside function
175void
176Pulse::FitSectionWise(
177 TString fitName,
178 TString fitOptions,
179 int fitMin,
180 int fitMax
181 )
182{
183 mFitMin = fitMin;
184 mFitMax = fitMax;
185 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 7 );
186
187 // ======================================================================
188 // Calculate startvaues for fit
189 double tau = 30.0;
190 double bsl = 0.8; //shall be around 0
191 double gain = 10; //startvalue for gain
192
193 int first_bin = mHisto->GetXaxis()->GetFirst();
194
195 for (int i = 0; i < 30; i++){
196 bsl += mHisto->GetBinContent(first_bin+0);
197 }
198 bsl = bsl/30;
199
200 double stop = mHisto->GetMaximumBin(); //pos of max
201 double amplitude = mHisto->GetBinContent(stop);
202 double start = stop-10; //pos 10 slices before maximum
203 int phe = 0;
204 if (mOrder < 0)
205 {
206 phe = (amplitude-bsl)/gain;
207 }
208 else
209 {
210 phe = mOrder+1;
211 }
212 // ======================================================================
213
214 double para[] = {bsl, amplitude, start, stop, tau, tau, phe};
215
216
217
218 mModel->SetParameters(para);
219 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2", "PhE");
220 mModel->SetLineColor(kRed);
221
222 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
223 mBsl.first = mModel->GetParameter(0);
224 mBsl.second = mModel->GetParError(0);
225 mHeight.first = mModel->GetParameter(1);
226 mHeight.second = mModel->GetParError(1);
227 mT0.first = mModel->GetParameter(2);
228 mT0.second = mModel->GetParError(2);
229 mT1.first = mModel->GetParameter(3);
230 mT1.second = mModel->GetParameter(3);
231 mTau1.first = mModel->GetParameter(4);
232 mTau1.second = mModel->GetParError(4);
233 mTau2.first = mModel->GetParameter(5);
234 mTau2.second = mModel->GetParError(5);
235 mPhE.first = mModel->GetParameter(6);
236 mPhE.second = mModel->GetParError(6);
237}
238
239// Fit pulseshape with a continous funciton expressed by a product of two exp-functions
240void
241Pulse::FitContious(
242 TString fitName,
243 TString fitOptions,
244 int fitMin,
245 int fitMax
246 )
247{
248 mFitMin = fitMin;
249 mFitMax = fitMax;
250 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 6 );
251
252 // ======================================================================
253 // Calculate startvaues for fit
254 double tau = 30.0;
255 double bsl = 0; //MW der ersten zehn slices
256 double gain = 10; //startvalue for gain
257 int first_bin = mHisto->GetXaxis()->GetFirst();
258
259 for (int i = 0; i < 30; i++){
260 bsl += mHisto->GetBinContent(first_bin+0);
261 }
262 bsl = bsl/30;
263
264 double stop = mHisto->GetMaximumBin(); //pos of max
265 double amplitude= mHisto->GetBinContent(stop);
266 double start = stop-10; //pos 10 slices before maximum
267 int phe = 0;
268 if (mOrder < 0)
269 {
270 phe = (amplitude-bsl)/gain;
271 }
272 else
273 {
274 phe = mOrder+1;
275 }
276 // ======================================================================
277
278
279 double para[] = {bsl, amplitude, start, tau, tau, phe};
280 mModel->SetParameters(para);
281 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2", "PhE");
282 mModel->SetLineColor(kBlue);
283
284 mFitResultPtr = mHisto->Fit(mModel, fitOptions);
285 mBsl.first = mModel->GetParameter(0);
286 mBsl.second = mModel->GetParError(0);
287 mHeight.first = mModel->GetParameter(1);
288 mHeight.second = mModel->GetParError(1);
289 mT0.first = mModel->GetParameter(2);
290 mT0.second = mModel->GetParError(2);
291 mTau1.first = mModel->GetParameter(3);
292 mTau1.second = mModel->GetParError(3);
293 mTau2.first = mModel->GetParameter(4);
294 mTau2.second = mModel->GetParError(4);
295 mPhE.first = mModel->GetParameter(5);
296 mPhE.second = mModel->GetParError(5);
297
298// cout << "test" << mBsl.first << "\t +/-" << mBsl.second << endl;
299}
300
301void
302Pulse::CalculateParameters()
303{
304 mIntegral.first = mModel->Integral(mFitMin, mFitMax);
305 //TODO setze konstantes glied auf 0 dann integrieren und fehler bestimmen
306 mIntegral.second = mModel->IntegralError(mFitMin, mFitMax);
307 mAmplitude.first = mModel->GetMaximum() - mBsl.first;
308// mAmplitude.second = mModel->GetMaximum() - mBsl.first;
309 mRisetime.first = mTimeResolution*CalculateHistEdgeRisetime();
310 mFitProb = mFitResultPtr->Prob();
311 mFitNCalls = mFitResultPtr->NCalls();
312 mFitNdf = mFitResultPtr->Ndf();
313 mChi2 = mFitResultPtr->Chi2();
314}
315
316double
317Pulse::CalculateHistEdgeRisetime() // in bins
318{
319 //rise-tim for pulse for pulse height from 10% to 90%
320 double riseTime = 0;
321
322 // Get Pulse Maximum value and bin
323 int maxBin = mHisto->GetMaximumBin();
324 double max = mHisto->GetBinContent(maxBin);
325
326 // values of bin with 10% and 90%
327 int tenPctBin = 0;
328 int ninetyPctBin = 0;
329
330 // flag points
331 bool finishedLow = false;
332 bool finishedHigh = false;
333
334 //2nd runvariable
335 int j = maxBin;
336
337 //loop over from 20 bins bevore max to maximum
338 for (int i = maxBin - 20; i<maxBin; i++){
339 j--;
340 if (j<0) break;
341
342 if (mHisto->GetBinContent(i) >= 0.1*max && !finishedLow)
343 {
344 tenPctBin = i;
345 finishedLow = true;
346 }
347
348
349 if (mHisto->GetBinContent(j) <= 0.9*max && !finishedHigh)
350 {
351 ninetyPctBin = j;
352 finishedHigh = true;
353 }
354
355 if (finishedLow && finishedHigh)
356 {
357 break;
358 }
359 }
360
361 riseTime = mHisto->GetBinLowEdge(ninetyPctBin) - mHisto->GetBinLowEdge(tenPctBin);
362 return riseTime;
363}
364
365
366
367
368// ===========================================================================
369// ACCESS
370// ===========================================================================
371
372TString Pulse::GetName(){ return mName;}
373double Pulse::GetBsl(){ return mBsl.first;}
374double Pulse::GetBslErr(){ return mBsl.second;}
375double Pulse::GetHeight(){ return mHeight.first;}
376double Pulse::GetHeightErr(){ return mHeight.second;}
377double Pulse::GetT0(){ return mT0.first;}
378double Pulse::GetT0Err(){ return mT0.second;}
379double Pulse::GetT1(){ return mT1.first;}
380double Pulse::GetT1Err(){ return mT1.second;}
381double Pulse::GetTau1(){ return mTau1.first;}
382double Pulse::GetTau1Err(){ return mTau1.second;}
383double Pulse::GetTau2(){ return mTau2.first;}
384double Pulse::GetTau2Err(){ return mTau2.second;}
385double Pulse::GetIntegral(){ return mIntegral.first;}
386double Pulse::GetIntegralErr(){ return mIntegral.second;}
387double Pulse::GetAmplitude(){ return mAmplitude.first;}
388double Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
389int Pulse::GetPhE(){ return mPhE.first;}
390int Pulse::GetPhEErr(){ return mPhE.second;}
391double Pulse::GetRiseTime(){ return mRisetime.first; }
392double Pulse::GetRiseTimeErr(){ return mRisetime.second;}
393int Pulse::GetOrder(){ return mOrder;}
394int Pulse::GetType(){return mType;}
395int Pulse::GetFitMin(){return mFitMin;}
396int Pulse::GetFitMax(){return mFitMax;}
397double Pulse::GetFitProb(){return mFitProb;}
398double Pulse::GetFitNCalls(){return mFitNCalls;}
399double Pulse::GetFitNdf(){return mFitNdf;}
400double Pulse::GetChi2(){return mChi2;}
401double Pulse::GetTimeResolution(){return mTimeResolution;}
402void Pulse::SetTimeResolution(double timeResolution){mTimeResolution = timeResolution;}
403
404
405// ===========================================================================
406// NON CLASS MATH FUNCTIONS
407// ===========================================================================
408
409//double
410int Heaviside(double val)
411{
412 if( val >= 0) return 1;
413 else return 0;
414}
415
416//Pulse::
417double shapeFunc(
418 double* t,
419 double* par
420 )
421{
422 double returnval= 0.0;
423
424 // name parameters
425 double bsl = par[0];
426 double height = par[1];
427 double start = par[2];
428// double rising = par[3];
429 double stop = par[3];
430 double tau1 = par[4];
431 double tau2 = par[5];
432 int phe = par[6];
433
434 // helper variables
435// double max = start + rising * tau1;
436 double e1 = (1 - Exp(-(t[0]-start)/tau1 ) );
437// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
438 double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
439 // calculate return value
440 returnval += Heaviside(t[0]-start)*e1;
441 returnval += Heaviside(t[0]-stop)*e2;
442 returnval *= phe*height;
443 returnval += bsl;
444// if (t[0] > start)
445// returnval += e1;
446// if (t[0] > stop)
447// returnval += e2;
448
449 return returnval;
450
451
452}
453
454
455
456/*
457[18:42:13] Dominik: können wir nachher skypen? haste n headset?
458[18:47:34] Dominik: hallo?
459[18:53:38] Dominik: so also ich hab das file glaubich gefunden
460[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
461[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
462[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
463[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
464[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
465[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
466[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.
467[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
468[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.
469Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
470
471p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
472
473bsl = p[0] # just the baseline
474 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
475 # it *would* be the pulse height only in case the
476 # recharging process would set in later...
477
478 start = p[2] # start of avalanche aka start of discharge of diode capacitance
479 stop = p[3] # stop of avalanche, or start of recharging of diode
480 tau1 = p[4] # time constant of discharging process
481 tau2 = p[5] # time constant of recharging process
482
483*/
484
485
486double shapeFunc2(
487 double* t,
488 double* par
489 )
490{
491 double returnval= 0.0;
492
493 // name parameters
494 double bsl = par[0];
495 double gain = par[1];
496 double t_0 = par[2];
497 double tau1 = par[3];
498 double tau2 = par[4];
499 int phe = par[5];
500
501 // helper variables
502 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
503 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
504
505 // calculate return value
506 returnval += bsl;
507 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0);
508
509 return returnval;
510}
Note: See TracBrowser for help on using the repository browser.