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

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