Project
Loading...
Searching...
No Matches
CaloRawFitterGS.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15#include <gsl/span>
16
19
20using namespace o2::phos;
27
29{
31 // prepare fitting arrays, once per lifetime
33 ma0[0] = 0.;
34 ma1[0] = 0.;
35 ma2[0] = 0.;
36 ma3[0] = 0.;
37 ma4[0] = 0.;
38 for (int i = 1; i < NMAXSAMPLES; i++) {
39 double xi = k * i;
40 mexp[i] = exp(-xi);
41 double s = mexp[i] * mexp[i];
42 ma0[i] = ma0[i - 1] + s;
43 ma1[i] = ma1[i - 1] + s * xi;
44 ma2[i] = ma2[i - 1] + s * xi * xi;
45 ma3[i] = ma3[i - 1] + s * xi * xi * xi;
46 ma4[i] = ma4[i - 1] + s * xi * xi * xi * xi;
47 }
48}
49
50CaloRawFitterGS::FitStatus CaloRawFitterGS::evaluate(gsl::span<short unsigned int> signal)
51{
52
53 // Pedestal analysis mode
54 if (mPedestalRun) {
55 int nPed = signal.size();
56 mAmp = 0.;
57 mTime = 0.;
58 for (auto a : signal) {
59 mAmp += a;
60 mTime += a * a;
61 }
62 if (nPed > 0) {
63 mAmp /= nPed;
64 mTime = mTime / nPed - mAmp * mAmp;
65 if (mTime > 0.) {
66 mTime = sqrt(mTime);
67 }
68 }
69 mOverflow = false;
70 mStatus = kOK;
71 return kOK;
72 }
73
75 // Extract amplitude and time
76 mStatus = evalFit(signal);
77 return mStatus;
78}
79
80CaloRawFitterGS::FitStatus CaloRawFitterGS::evalFit(gsl::span<short unsigned int> signal)
81{
82 // Calculate signal parameters (energy, time, quality) from array of samples
83 // Fit with semi-gaus function with free parameters time and amplitude
84 // Signal overflows if there are at least 3 samples of the same amplitude above 900
85
86 // Calculate signal parameters (energy, time, quality) from array of samples
87 // Energy is a maximum sample minus pedestal 9
88 // Time is the first time bin
89 // Signal overflows is there are at least 3 samples of the same amplitude above 900
90
91 int nSamples = signal.size();
92 if (nSamples == 0) {
93 mAmp = 0;
94 mTime = 0.;
95 mChi2 = 0.;
96 return kEmptyBunch;
97 }
98 if (nSamples == 1) {
99 mAmp = signal[0];
100 mTime = 0.;
101 mChi2 = 1.;
102 return kOK;
103 }
104
105 mOverflow = false;
106
107 // if pedestal should be subtracted first evaluate it
108 float pedMean = 0;
109 int nPed = 0;
110 if (mPedSubtract) {
111 // remember inverse time order
112 for (auto it = signal.rbegin(); (nPed < mPreSamples) && it != signal.rend(); ++it) {
113 nPed++;
114 pedMean += *it;
115 }
116 if (nPed > 0) {
117 pedMean /= nPed;
118 }
119 nSamples -= mPreSamples;
120 if (nSamples <= 0) { // empty bunch left
121 mAmp = 0;
122 mTime = 0.;
123 mChi2 = 0.;
124 return kEmptyBunch;
125 }
126 }
127
128 float maxSample = 0.; // maximal sample value
129 int nMax = 0; // number of consequitive maximal samples
130 bool spike = false; // spike in previoud signal bin?
131 short ap = -1, app = -1; // remember previous values to evaluate spikes
132 double b0 = 0., b1 = 0., b2 = 0., y2 = 0.; // fit coeficients
133 double sa0 = 0., sa1 = 0., sa2 = 0., sa3 = 0., sa4 = 0.; // corrections in case of overflow
134
135 int firstS = nSamples;
136 int j = TMath::Min(nSamples, NMAXSAMPLES - 1);
137 for (int i = 1; i <= j; i++) {
138 short a = signal[firstS - i] - pedMean; // remember inverse order of samples
139 float xi = i * mDecTime;
140 if (a > maxSample) {
141 maxSample = a;
142 nMax = 1;
143 } else {
144 if (a == maxSample) {
145 nMax++;
146 }
147 }
148 // check if there was a spike in previous step?
149 if (app > 0 && ap > 0) {
150 spike = (ap - a > mSpikeThreshold) && (ap - app > mSpikeThreshold);
151 }
152 if (spike) { // Try to recover: subtract last point contribution and replace by average of "app" and "a"
153 float atmp = 0.5 * (app + a);
154 float xiprev = xi - mDecTime;
155 float st = (atmp - ap) * mexp[i - 1];
156 b0 += st;
157 b1 += st * xiprev;
158 b2 += st * xiprev * xiprev;
159 y2 += atmp * atmp - ap * ap;
160 ap = a;
161 } else {
162 app = ap;
163 ap = a;
164 }
165 // Check if in saturation
166 if (maxSample > 900 && nMax >= 3) {
167 // Remove overflow points from the fit
168 if (!mOverflow) { // first time in this sample: remove two previous points
169 sa0 = ma0[i] - ma0[i - 2]; // can not appear at i<2
170 sa1 = ma1[i] - ma1[i - 2];
171 sa2 = ma2[i] - ma2[i - 2];
172 sa3 = ma3[i] - ma3[i - 2];
173 sa4 = ma4[i] - ma4[i - 2];
174 float st = ap * mexp[i - 1];
175 float xiprev = xi - mDecTime;
176 b0 -= st;
177 b1 -= st * xiprev;
178 b2 -= st * xiprev * xiprev;
179 y2 -= ap * ap;
180 st = app * mexp[i - 2];
181 xiprev -= mDecTime;
182 b0 -= st;
183 b1 -= st * xiprev;
184 b2 -= st * xiprev * xiprev;
185 y2 -= ap * ap;
186 }
187 mOverflow = true;
188 }
189 if (!mOverflow) {
190 // to calculate time
191 float st = a * mexp[i];
192 b0 += st;
193 b1 += st * xi;
194 b2 += st * xi * xi;
195 y2 += a * a;
196 } else { // do not add current point and subtract contributions to amx[] arrays
197 sa0 = ma0[i] - ma0[i - 1]; // can not appear at i<2
198 sa1 = ma1[i] - ma1[i - 1];
199 sa2 = ma2[i] - ma2[i - 1];
200 sa3 = ma3[i] - ma3[i - 1];
201 sa4 = ma4[i] - ma4[i - 1];
202 }
203 } // Scanned full
204
205 // too small amplitude, assing max to max Amp and time to zero and do not calculate height
206 if (maxSample < mMinTimeCalc) {
207 mAmp = maxSample;
208 mTime = 0.;
209 mChi2 = 0.;
210 return kOK;
211 }
212
213 if (mOverflow && b0 == 0) { // strong overflow, no reasonable counts, can not extract anything
214 mAmp = 0.;
215 mTime = 0.;
216 mChi2 = 900.;
217 return kOverflow;
218 }
219
220 // calculate time, amp and chi2
221 double a, b, c, d, e; // Polinomial coefficients
222 if (!mOverflow) {
223 a = ma1[j] * b0 - ma0[j] * b1;
224 b = ma0[j] * b2 + 2. * ma1[j] * b1 - 3. * ma2[j] * b0;
225 c = 3. * (ma3[j] * b0 - ma1[j] * b2);
226 d = 3. * ma2[j] * b2 - ma4[j] * b0 - 2. * ma3[j] * b1;
227 e = ma4[j] * b1 - ma3[j] * b2;
228 } else { // account removed points in overflow
229 a = (ma1[j] - sa1) * b0 - (ma0[j] - sa0) * b1;
230 b = (ma0[j] - sa0) * b2 + 2. * (ma1[j] - sa1) * b1 - 3. * (ma2[j] - sa2) * b0;
231 c = 3. * ((ma3[j] - sa3) * b0 - (ma1[j] - sa1) * b2);
232 d = 3. * (ma2[j] - sa2) * b2 - (ma4[j] - sa4) * b0 - 2. * (ma3[j] - sa4) * b1;
233 e = (ma4[j] - sa4) * b1 - (ma3[j] - sa3) * b2;
234 }
235
236 // Find zero of 4-order polinomial
237 // first use linear extrapolation to reach correct root of four
238 double z = -1.;
239 if (ma0[j] * b1 - ma1[j] * b0 != 0) {
240 z = (ma1[j] * b1 - ma2[j] * b0) / (ma0[j] * b1 - ma1[j] * b0) - 1.; // linear fit + offset
241 }
242 double q = 0., dq = 0., ddq = 0., lq = 0., dz = 0.1;
243 double z2 = z * z;
244 double z3 = z2 * z;
245 double z4 = z2 * z2;
246 q = a * z4 + b * z3 + c * z2 + d * z + e; // polinomial
247 dq = 4. * a * z3 + 3. * b * z2 + 2. * c * z + d; // Derivative
248 ddq = 12. * a * z2 + 6. * b * z + 2. * c; // Second derivative
249 if (dq != 0.) {
250 lq = q * ddq / (dq * dq);
251 }
252 // dz = -q/dq ; // Newton ~7 terations
253 // dz =-(1+0.5*lq)*q/dq ; // Chebyshev ~3 iterations to reach |q|<1.e-11
254 double ttt = dq * (1. - 0.5 * lq); // Halley’s method ~3 iterations, a bit more precise
255 if (ttt != 0) {
256 dz = -q / ttt;
257 } else {
258 dz = 0.1; // step off saddle point
259 }
260 int it = 0;
261 while (TMath::Abs(q) > 0.0001 && (++it < 15)) {
262 z += dz;
263 z2 = z * z;
264 z3 = z2 * z;
265 z4 = z2 * z2;
266 q = a * z4 + b * z3 + c * z2 + d * z + e;
267 dq = 4. * a * z3 + 3. * b * z2 + 2. * c * z + d;
268 ddq = 12. * a * z2 + 6. * b * z + 2. * c;
269 if (dq != 0) {
270 lq = q * ddq / (dq * dq);
271 ttt = dq * (1. - 0.5 * lq);
272 // dz = -q/dq ; //Newton
273 // dz =-(1+0.5*lq)*q/dq ; //Chebyshev
274 if (ttt != 0) {
275 dz = -q / ttt; // Halley’s
276 } else {
277 dz = -q / dq;
278 }
279 } else {
280 dz = 0.5 * dz; // step off saddle point
281 }
282 }
283
284 // check that result is reasonable
285 double denom = ma4[j] - 4. * ma3[j] * z + 6. * ma2[j] * z * z - 4. * ma1[j] * z * z * z + ma0[j] * z * z * z * z;
286 if (denom != 0.) {
287 mAmp = 4. * exp(-2 - z) * (b2 - 2. * b1 * z + b0 * z * z) / denom;
288 } else {
289 mAmp = 0.;
290 }
291
292 if ((TMath::Abs(q) < mQAccuracy) && (mAmp < 1.2 * maxSample)) { // converged and estimated amplitude is not mush larger than Max
293 mTime = z / mDecTime;
294 mChi2 = (y2 - 0.25 * exp(2. + z) * mAmp * (b2 - 2 * b1 * z + b0 * z2)) / nSamples;
295 return kOK;
296 } else { // too big difference, fit failed
297 mAmp = maxSample;
298 mTime = 0; // First count in sample
299 mChi2 = 999.;
300 return kFitFailed;
301 }
302}
uint64_t exp(uint64_t base, uint8_t exp) noexcept
int32_t i
const GPUTPCGMMerger::trackCluster & b1
uint32_t j
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
benchmark::State & st
FitStatus evaluate(gsl::span< short unsigned int > signal) final
Evaluation Amplitude and TOF.
FitStatus evalFit(gsl::span< short unsigned int > signal)
static constexpr int NMAXSAMPLES
float mAmp
amplitude of last processed sample
bool mOverflow
is last sample saturated
bool mPedSubtract
should one evaluate and subtract pedestals
float mTime
time of last processed sample
float mChi2
chi2 calculated in last fit
bool mPedestalRun
analyze as pedestal run
FitStatus mStatus
status of last evaluated sample: -1: not yet evaluated; 0: OK; 1: overflow; 2: too large RMS; 3: sing...
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
float mSampleDecayTime
Time parameter in Gamma2 function (1/tau, 100.e-9/2.1e-6)