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

Last change on this file since 14472 was 14472, checked in by Jens Buss, 12 years ago
class for pulse fit modell
  • Property svn:executable set to *
File size: 7.8 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 mHisto = NULL;
13 mOptions = "";
14 mBsl = 0;
15 mHeight = 0;
16 mStart = 0;
17 mRising = 0;
18 mTauRising = 0;
19 mTauFalling = 0;
20 mIntegral = 0;
21 mAmplitude = 0;
22 mPhE = 0;
23}
24
25Pulse::Pulse(TH1* histo)
26{
27 mHisto = histo;
28 mOptions = "";
29 mBsl = 0;
30 mHeight = 0;
31 mStart = 0;
32 mRising = 0;
33 mTauRising = 0;
34 mTauFalling = 0;
35 mIntegral = 0;
36 mAmplitude = 0;
37 mPhE = 0;
38}
39
40Pulse::Pulse(TH1* histo, TString options)
41{
42 mHisto = histo;
43 mOptions = options;
44 mBsl = 0;
45 mHeight = 0;
46 mStart = 0;
47 mRising = 0;
48 mTauRising = 0;
49 mTauFalling = 0;
50 mIntegral = 0;
51 mAmplitude = 0;
52 mPhE = 0;
53}
54
55Pulse::~Pulse()
56{
57
58}
59
60void
61Pulse::Fit(
62 TString fitName,
63 TString fitOptions
64 )
65{
66 int fitMin = mHisto->GetXaxis()->GetFirst();
67 int fitMax = mHisto->GetXaxis()->GetLast();
68
69 Fit(
70 fitName,
71 fitOptions,
72 fitMin,
73 fitMax
74 );
75}
76
77void
78Pulse::Fit(
79 TString fitName,
80 TString fitOptions,
81 int fitMin,
82 int fitMax
83 )
84{
85 TF1 fit(fitName, shapeFunc, fitMin, fitMax, 6 );
86
87 // ======================================================================
88 // Calculate startvaues for fit
89 double tau = 30.0;
90 double bsl = 0; //MW der ersten zehn slices
91
92 int first_bin = mHisto->GetXaxis()->GetFirst();
93
94 for (int i = 0; i < 10; i++){
95 bsl += mHisto->GetBinContent(first_bin+0);
96 }
97 bsl = bsl/10;
98
99 double stop = mHisto->GetMaximumBin(); //pos of max
100 double height = mHisto->GetBinContent(stop);
101 double start = stop-10; //pos 10 slices before maximum
102 double rising = (stop - start)/tau;
103 // ======================================================================
104
105 double para[] = {bsl, height, start, stop, tau, tau};
106
107
108
109 fit.SetParameters(para);
110 fit.SetParNames("BSL", "Height", "Start", "Rising", "Tau1", "TauFalling");
111 fit.SetLineColor(kRed);
112 mHisto->Fit(&fit, fitOptions);
113
114 mBsl = fit.GetParameter(0);
115 mHeight = fit.GetParameter(1);
116 mStart = fit.GetParameter(2);
117 mRising = fit.GetParameter(3);
118 mTauRising = fit.GetParameter(4);
119 mTauFalling = fit.GetParameter(5);
120}
121
122void
123Pulse::Fit2(
124 TString fitName,
125 TString fitOptions
126 )
127{
128 int fitMin = mHisto->GetXaxis()->GetFirst();
129 int fitMax = mHisto->GetXaxis()->GetLast();
130
131 Fit2(
132 fitName,
133 fitOptions,
134 fitMin,
135 fitMax
136 );
137}
138
139void
140Pulse::Fit2(
141 TString fitName,
142 TString fitOptions,
143 int fitMin,
144 int fitMax
145 )
146{
147 TF1 fit(fitName, shapeFunc2, fitMin, fitMax, 5 );
148
149 // ======================================================================
150 // Calculate startvaues for fit
151 double tau = 30.0;
152 double bsl = 0; //MW der ersten zehn slices
153
154 int first_bin = mHisto->GetXaxis()->GetFirst();
155
156 for (int i = 0; i < 10; i++){
157 bsl += mHisto->GetBinContent(first_bin+0);
158 }
159 bsl = bsl/10;
160
161 double stop = mHisto->GetMaximumBin(); //pos of max
162 double amplitude= mHisto->GetBinContent(stop);
163 double start = stop-10; //pos 10 slices before maximum
164 // ======================================================================
165
166
167 double para[] = {bsl, amplitude, start, tau, tau};
168 fit.SetParameters(para);
169 fit.SetParNames("BSL", "Amplitude", "Start", "Tau1", "TauFalling");
170 fit.SetLineColor(kBlue);
171 mHisto->Fit(&fit, fitOptions);
172
173 mBsl = fit.GetParameter(0);
174 mAmplitude = fit.GetParameter(1);
175 mStart = fit.GetParameter(2);
176 mTauRising = fit.GetParameter(3);
177 mTauFalling = fit.GetParameter(4);
178}
179
180double Pulse::GetBsl(){ return mBsl;}
181double Pulse::GetHeight(){ return mHeight;}
182double Pulse::GetAvalancheStart(){ return mStart;}
183double Pulse::GetAvalancheEnd(){ return mRising;}
184double Pulse::GetTimeConstRising(){ return mTauRising;}
185double Pulse::GetTimeConstFalling(){ return mTauFalling;}
186double Pulse::GetIntegral(){ return mIntegral;}
187double Pulse::GetAmplitude(){ return mAmplitude;}
188int Pulse::GetPE(){ return mPhE;}
189
190//double
191int Heaviside(double val)
192{
193 if( val < 0) return 0;
194 if( val >= 0) return 1;
195}
196
197//Pulse::
198double shapeFunc(
199 double* t,
200 double* par
201 )
202{
203 double returnval= 0.0;
204
205 // name parameters
206 double bsl = par[0];
207 double height = par[1];
208 double start = par[2];
209// double rising = par[3];
210 double tau1 = par[4];
211 double tau2 = par[5];
212 double stop = par[3];
213
214 // helper variables
215// double max = start + rising * tau1;
216 double e1 = height * (1 - Exp(-(t[0]-start)/tau1 ) );
217// double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
218 double e2 = (-1) * height * (1 - Exp(-(t[0]-stop)/tau2 ) );
219 // calculate return value
220 returnval += bsl;
221
222 if (t[0] > start)
223 returnval += e1;
224 if (t[0] > stop)
225 returnval += e2;
226
227 return returnval;
228}
229
230/*
231[18:42:13] Dominik: können wir nachher skypen? haste n headset?
232[18:47:34] Dominik: hallo?
233[18:53:38] Dominik: so also ich hab das file glaubich gefunden
234[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
235[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
236[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
237[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
238[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
239[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
240[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.
241[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
242[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.
243Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.
244
245p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters
246
247bsl = p[0] # just the baseline
248 height = p[1] # this is not the pulse height, but the maximum of the discharging edge
249 # it *would* be the pulse height only in case the
250 # recharging process would set in later...
251
252 start = p[2] # start of avalanche aka start of discharge of diode capacitance
253 stop = p[3] # stop of avalanche, or start of recharging of diode
254 tau1 = p[4] # time constant of discharging process
255 tau2 = p[5] # time constant of recharging process
256
257*/
258
259
260double shapeFunc2(
261 double* t,
262 double* par
263 )
264{
265 double returnval= 0.0;
266
267 // name parameters
268 double bsl = par[0];
269 double gain = par[1];
270 double t_0 = par[2];
271 double tau1 = par[3];
272 double tau2 = par[4];
273
274 // helper variables
275 double e1 = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
276 double e2 = Exp(-(t[0]-t_0)/tau2 ) ;
277
278 // calculate return value
279 returnval += bsl;
280 returnval += gain*e1*e2*Heaviside(t[0]-t_0);
281
282 return returnval;
283}
284
Note: See TracBrowser for help on using the repository browser.