Project
Loading...
Searching...
No Matches
CaloRawFitter.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#include <cmath>
17
20
21using namespace o2::phos;
22
30
31CaloRawFitter::FitStatus CaloRawFitter::evaluate(gsl::span<short unsigned int> signal)
32{
33
34 //Pedestal analysis mode
35 if (mPedestalRun) {
36 int nPed = signal.size();
37 float mean = 0.;
38 float rms = 0.;
39 for (auto a : signal) {
40 mean += a;
41 rms += a * a;
42 }
43 if (nPed > 0) {
44 mean /= nPed;
45 rms = rms / nPed - mean * mean;
46 if (rms > 0.) {
47 rms = sqrt(rms);
48 }
49 }
50 mAmp = mean;
51 mTime = rms; // only in Pedestal mode!
52 mOverflow = false;
53 mStatus = kOK;
54 return kOK;
55 }
56
57 // Extract amplitude and time using maximum and k-level methods
58 return evalKLevel(signal);
59}
60
61CaloRawFitter::FitStatus CaloRawFitter::evalKLevel(gsl::span<short unsigned int> signal) //const ushort *signal, int sigStart, int sigLength)
62{
63 // Calculate signal parameters (energy, time, quality) from array of samples
64 // Energy is a maximum sample minus pedestal 9
65 // Time is the first time bin
66 // Signal overflows is there are at least 3 samples of the same amplitude above 900
67
68 int sigLength = signal.size();
69 if (sigLength == 0) {
70 return kEmptyBunch;
71 }
72
73 float pedMean = 0;
74 int nPed = 0;
75 mMaxSample = 0;
76 int nMax = 0; //number of consequitive maximal samples
77 bool spike = false;
78 mOverflow = false;
79
80 int ap = -1, app = -1; //remember previous values to evaluate spikes
81 for (auto it = signal.rbegin(); it != signal.rend(); ++it) {
82 uint16_t a = *it;
83 if (mPedSubtract) {
84 if (nPed < mPreSamples) { //inverse signal time order
85 nPed++;
86 pedMean += a;
87 }
88 }
89 if (a > mMaxSample) {
90 mMaxSample = a;
91 nMax = 0;
92 }
93 if (a == mMaxSample) {
94 nMax++;
95 }
96 //check if there is a spike
97 if (app >= 0 && ap >= 0) {
98 spike |= (2 * ap - (a + app) > 2 * mSpikeThreshold);
99 }
100 app = ap;
101 ap = a;
102 }
103 mAmp = (float)mMaxSample;
104
105 if (spike) {
106 mTime = -2;
107 mOverflow = false;
108 return kSpike;
109 }
110
111 if (mMaxSample > 900 && nMax > 2) {
112 mOverflow = true;
113 }
114
115 float pedestal = 0;
116 if (mPedSubtract) {
117 if (nPed > 0) {
118 pedMean /= nPed;
119 } else {
120 mAmp = 0.;
121 mTime = -2;
122 mOverflow = false;
123 return kBadPedestal;
124 }
125 } else {
126 pedMean = 0.;
127 }
128
129 mAmp -= pedMean;
130 if (mAmp < mBaseLine) {
131 mAmp = 0;
132 }
133
134 //Evaluate time
135 mTime = -2;
136 const int nLine = 6; //Parameters of fitting
137 const float eMinTOF = 10.; //Choosed from beam-test and cosmic analyis
138 const float kAmp = 0.35; //Result slightly depends on them, so no getters
139 // Avoid too low peak:
140 if (mAmp < eMinTOF) {
141 return kOK; //use estimated time
142 }
143
144 // Find index posK (kLevel is a level of "timestamp" point Tk):
145 int posK = sigLength - 1; //last point before crossing k-level
146 float levelK = pedestal + kAmp * mAmp;
147 while (posK >= 0 && signal[posK] <= levelK) {
148 posK--;
149 }
150 posK++;
151
152 if (posK == 0 || posK == sigLength - 1) {
153 return kNoTime; //
154 }
155
156 // Find crossing point by solving linear equation (least squares method)
157 int np = 0;
158 int iup = posK - 1;
159 int idn = posK;
160 double sx = 0., sy = 0., sxx = 0., sxy = 0.;
161 double x, y;
162
163 while (np < nLine) {
164 //point above crossing point
165 if (iup >= 0) {
166 x = sigLength - iup - 1;
167 y = signal[iup];
168 sx += x;
169 sy += y;
170 sxx += (x * x);
171 sxy += (x * y);
172 np++;
173 iup--;
174 }
175 //Point below crossing point
176 if (idn < sigLength) {
177 if (signal[idn] < pedestal) {
178 idn = sigLength - 1; //do not scan further
179 idn++;
180 continue;
181 }
182 x = sigLength - idn - 1;
183 y = signal[idn];
184 sx += x;
185 sy += y;
186 sxx += (x * x);
187 sxy += (x * y);
188 np++;
189 idn++;
190 }
191 if (idn >= sigLength && iup < 0) {
192 break; //can not fit futher
193 }
194 }
195
196 double det = np * sxx - sx * sx;
197 if (det == 0) {
198 return kNoTime;
199 }
200 if (np == 0) {
201 return kEmptyBunch;
202 }
203 double c1 = (np * sxy - sx * sy) / det; //slope
204 double c0 = (sy - c1 * sx) / np; //offset
205 if (c1 == 0) {
206 return kNoTime;
207 }
208
209 // Find where the line cross kLevel:
210 mTime += (levelK - c0) / c1 - 5.; //5: mean offset between k-Level and start times
211
212 if (mOverflow) {
213 return kOverflow;
214 } else {
215 return kOK;
216 }
217}
bool const GPUTPCGMMerger::trackCluster * c1
FitStatus evalKLevel(gsl::span< short unsigned int > signal)
float mAmp
amplitude of last processed sample
bool mOverflow
is last sample saturated
short mMaxSample
maximal sample
bool mPedSubtract
should one evaluate and subtract pedestals
float mTime
time of last processed sample
virtual FitStatus evaluate(gsl::span< short unsigned int > signal)
Evaluation Amplitude and TOF return status -1: not evaluated/empty bunch; 0: OK; 1: overflow; 4: sing...
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...
GLint GLenum GLint x
Definition glcorearb.h:403
GLint y
Definition glcorearb.h:270
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
short mSpikeThreshold
Single spike >100 ADC channels.
short mPreSamples
number of pre-samples readout before sample (if no pedestal subtrauction)