ATLAS Offline Software
Loading...
Searching...
No Matches
CscBipolarStripFitter.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3*/
4
5// CscBipolarStripFitter.cxx
6#include "CscBipolarStripFitter.h"
7
8#include "CscCalibTools/ICscCalibTool.h"
9#include "MuonIdHelpers/CscIdHelper.h"
10#include "MuonPrepRawData/CscStripPrepData.h"
11#include "MuonReadoutGeometry/MuonDetectorManager.h"
12
14// K. Nikolopoulos
15// V. Chernyatin and M. Schernau
16// Oct 2007
17//
18// CSC Strip fitter using the bipolar waveform of the pulse
19//
20//
21
22using Muon::CscStripPrepData;
23
24using Result = ICscStripFitter::Result;
25using ChargeList = ICscStripFitter::ChargeList;
26
27//**********************************************************************
28
29CscBipolarStripFitter::CscBipolarStripFitter(const std::string &type, const std::string &aname, const IInterface *parent) :
30 AthAlgTool(type, aname, parent), m_phelper(nullptr), m_n(0), m_n2(0), m_zmax(0), m_bipolarNormalization(0), m_tsampling(0) {
31 declareInterface<ICscStripFitter>(this);
32 declareProperty("chargeError", m_qerr = 5500.0);
33 declareProperty("timeError", m_terr = 5.0);
34 declareProperty("failChargeError", m_qerr_fail = 5000.0);
35 declareProperty("failTimeError", m_terr_fail = 50.0);
36 declareProperty("chargeCalibrationError", m_qerrprop = 0.002);
37}
38
39//**********************************************************************
40
41StatusCode CscBipolarStripFitter::initialize() {
42 ATH_MSG_DEBUG("Initializing " << name());
43
44 ATH_MSG_DEBUG("Properties for " << name() << ":");
45 ATH_MSG_DEBUG(" Charge error: " << m_qerr);
46 ATH_MSG_DEBUG(" Time error: " << m_terr);
47 ATH_MSG_DEBUG(" Fail charge error: " << m_qerr_fail);
48 ATH_MSG_DEBUG(" Fail time error: " << m_terr_fail);
49 ATH_MSG_DEBUG(" Charge calibration error: " << m_qerrprop);
50
52 ATH_CHECK_RECOVERABLE(m_cscCalibTool.retrieve());
53
54 const MuonGM::MuonDetectorManager *muDetMgr = nullptr;
55 ATH_CHECK_RECOVERABLE(detStore()->retrieve(muDetMgr));
56 ATH_MSG_DEBUG("Retrieved geometry.");
57
58 m_phelper = muDetMgr->cscIdHelper();
59
60 m_n = m_cscCalibTool->getNumberOfIntegration(); // 12.;
61 m_n2 = m_cscCalibTool->getNumberOfIntegration2(); // 11.66;
62 m_zmax = m_n + m_n2 + 2. - sqrt((m_n + 2. + m_n2) * (m_n + 2. + m_n2) - 4. * m_n * (m_n2 + 1));
63 m_zmax *= 0.5;
64 // m_n+1 - sqrt(m_n+1.0);
65 m_bipolarNormalization = FindPow(m_zmax) * (1 - m_zmax / (m_n2 + 1)) * std::exp(-m_zmax);
66 // FindPow(m_zmax)*(1-m_zmax/(m_n+1))*std::exp(-m_zmax);
67 m_tsampling = m_cscCalibTool->getSamplingTime() / 2.; // 50/2 = 25.;
68
69 return StatusCode::SUCCESS;
70}
71
72//**********************************************************************
73
74Result CscBipolarStripFitter::fit(const ChargeList &chgs, double period, Identifier &stripId) const {
75 ATH_MSG_DEBUG("Fitting with sample time " << period);
76 ATH_MSG_DEBUG(" Number of charges =" << chgs.size() << " Charges: " << chgs);
77 Result res;
78
79 IdentifierHash stripHash;
80 if (m_phelper->get_channel_hash(stripId, stripHash)) {
81 ATH_MSG_WARNING("Unable to get CSC striphash id "
82 << " the identifier is ");
83 stripId.show();
84 }
85 ATH_MSG_DEBUG("CalibCscStripFitter:: " << stripId << " " << (unsigned int)stripHash);
86
91 double initValues[3];
92 for (int i = 0; i < 3; ++i) { initValues[i] = 0.; }
93 // int imax = -1;
94 initValues[2] = 8.; // m_width
96
97 double noise = m_cscCalibTool->stripNoise(stripHash);
98 res.status = 0;
99 if (m_qerr > 0)
100 res.dcharge = m_qerr;
101 else
102 res.dcharge = noise;
103
104 res.dtime = m_terr;
105 res.charge = initValues[0];
106 res.time = initValues[1];
107
108 initValues[1] = res.time - 2.5;
109 // int imeas =4;
110 // int meas[4] = {true,true,true,true};
111 // int ipar = 3;
112 // int par[3] = {true,true,false};
113 // double m_chi2;
114 //double result[3];
116
117 // Add an error proportional to the charge.
118 double dqprop = m_qerrprop * res.charge;
119 res.dcharge = sqrt(res.dcharge * res.dcharge + dqprop * dqprop);
120 // Display result.
121 ATH_MSG_DEBUG(" Status: " << res.status);
122 ATH_MSG_DEBUG(" Charge: " << res.charge);
123 ATH_MSG_DEBUG(" Charge err: " << res.dcharge);
124 ATH_MSG_DEBUG(" Time: " << res.time);
125 //ATH_MSG_DEBUG(" fit : " << result[0] << " " << result[1] << " " << result[2]);
126 return res;
127}
128
129double CscBipolarStripFitter::FindInitValues(double *x, double *initValues, int *maxsample) {
130 // find maximum sample imax:
131 double peakingTime = -99.; // interpolated peaking time in samples
132 double amplitude = -99.; // interpolated amplitude
133 double amin, amax;
134 int imax = 0;
135 const int numSample = 4;
136 amax = x[0];
137 amin = x[0];
138 for (int j = 1; j < numSample; j++) {
139 if (amax < x[j]) {
140 amax = x[j];
141 imax = j;
142 }
143 if (amin > x[j]) amin = x[j];
144 }
145
146 // calculate peak and amplitude:
147 (*maxsample) = imax;
148 if ((imax == 0) || (imax == numSample - 1)) // no interpolation possible
149 {
150 peakingTime = imax;
151 amplitude = amax;
152 } else // do parabolic interpolation
153 {
154 double a, b, c; // coeffients of parabola y=a*x*x+b*x+c
155 a = 0.5 * (x[imax + 1] + x[imax - 1] - 2.0 * x[imax]);
156 b = 0.5 * (x[imax + 1] - x[imax - 1]);
157 c = x[imax];
158
159 peakingTime = -0.5 * b / a;
160 amplitude = a * peakingTime * peakingTime + b * peakingTime + c;
161 peakingTime += imax; // it was relative to imax
162 }
163
164 initValues[0] = amplitude;
165 initValues[1] = peakingTime;
166 // - m_zmax*initValues[2]/m_tsampling;
167 return x[imax];
168}
169
170double CscBipolarStripFitter::FindPow(double z) const {
171 double zpower = std::exp(m_n * std::log(z));
172 return zpower;
173}
174
175void CscBipolarStripFitter::InvertMatrix(double matrix[][3], const int dim, const int *correspdim) {
176 // invert 2x2 or 3x3 symmetric matrix
177 if (dim == 1) {
178 int ii = correspdim[0];
179 matrix[ii][ii] = 1. / matrix[ii][ii];
180
181 } else if (dim == 2) {
182 int ii = correspdim[0];
183 int jj = correspdim[1];
184 double determinant = -matrix[jj][ii] * matrix[ii][jj] + matrix[ii][ii] * matrix[jj][jj];
185 // if(std::abs(determinant)<1.e-13)
186 // std::cout<<" zero determinant "<<std::endl;
187 double i00 = matrix[ii][ii];
188 matrix[ii][ii] = matrix[jj][jj] / determinant;
189 matrix[jj][jj] = i00 / determinant;
190 matrix[ii][jj] = -matrix[ii][jj] / determinant;
191 matrix[jj][ii] = matrix[ii][jj];
192
193 } else if (dim == 3) {
194 double sm12 = matrix[0][1] * matrix[0][1];
195 double sm13 = matrix[0][2] * matrix[0][2];
196 double sm23 = matrix[1][2] * matrix[1][2];
197 double determinant = matrix[0][0] * matrix[1][1] * matrix[2][2] - matrix[0][0] * sm23 - sm12 * matrix[2][2] +
198 2. * matrix[0][1] * matrix[0][2] * matrix[1][2] - matrix[1][1] * sm13;
199
200 // if(std::abs(determinant)<1.e-13)
201 // std::cout << "zero determinant"<<std::endl;
202 double i00 = matrix[1][1] * matrix[2][2] - sm23;
203 double i11 = matrix[0][0] * matrix[2][2] - sm13;
204 double i22 = matrix[0][0] * matrix[1][1] - sm12;
205
206 matrix[1][0] = (matrix[0][2] * matrix[1][2] - matrix[2][2] * matrix[0][1]) / determinant;
207 matrix[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1]) / determinant;
208 matrix[2][1] = (matrix[0][1] * matrix[0][2] - matrix[0][0] * matrix[1][2]) / determinant;
209 matrix[0][1] = matrix[1][0];
210 matrix[0][2] = matrix[2][0];
211 matrix[1][2] = matrix[2][1];
212 matrix[0][0] = i00 / determinant;
213 matrix[1][1] = i11 / determinant;
214 matrix[2][2] = i22 / determinant;
215
216 } else {
217 // int ierr;
218 // matrix.invert(ierr);
219 printf("this is not a 1x1 or 2x2 or 3x3 matrix\n");
220 }
221}
222
223void CscBipolarStripFitter::InvertSymmetric4x4(double W[][4]) {
224 double Determinant =
225 W[0][3] * W[0][3] * W[1][2] * W[1][2] - 2. * W[0][2] * W[0][3] * W[1][2] * W[1][3] + W[0][2] * W[0][2] * W[1][3] * W[1][3] -
226 W[0][3] * W[0][3] * W[1][1] * W[2][2] + 2. * W[0][1] * W[0][3] * W[1][3] * W[2][2] - W[0][0] * W[1][3] * W[1][3] * W[2][2] +
227 2. * W[0][2] * W[0][3] * W[1][1] * W[2][3] - 2. * W[0][1] * W[0][3] * W[1][2] * W[2][3] -
228 2. * W[0][1] * W[0][2] * W[1][3] * W[2][3] + 2. * W[0][0] * W[1][2] * W[1][3] * W[2][3] + W[0][1] * W[0][1] * W[2][3] * W[2][3] -
229 W[0][0] * W[1][1] * W[2][3] * W[2][3] - W[0][2] * W[0][2] * W[1][1] * W[3][3] + 2. * W[0][1] * W[0][2] * W[1][2] * W[3][3] -
230 W[0][0] * W[1][2] * W[1][2] * W[3][3] - W[0][1] * W[0][1] * W[2][2] * W[3][3] + W[0][0] * W[1][1] * W[2][2] * W[3][3];
231
232 W[0][1] = W[3][0] * W[3][1] * W[2][2] - W[3][0] * W[2][1] * W[3][2] - W[2][0] * W[3][1] * W[3][2] + W[1][0] * W[3][2] * W[3][2] +
233 W[2][0] * W[2][1] * W[3][3] - W[1][0] * W[2][2] * W[3][3];
234 W[0][2] = -W[3][0] * W[2][1] * W[3][1] + W[2][0] * W[3][1] * W[3][1] + W[3][0] * W[1][1] * W[3][2] - W[1][0] * W[3][1] * W[3][2] -
235 W[2][0] * W[1][1] * W[3][3] + W[1][0] * W[2][1] * W[3][3];
236 W[0][3] = W[3][0] * W[2][1] * W[2][1] - W[2][0] * W[2][1] * W[3][1] - W[3][0] * W[1][1] * W[2][2] + W[1][0] * W[3][1] * W[2][2] +
237 W[2][0] * W[1][1] * W[3][2] - W[1][0] * W[2][1] * W[3][2];
238 W[1][2] = W[3][0] * W[3][0] * W[2][1] - W[2][0] * W[3][0] * W[3][1] - W[1][0] * W[3][0] * W[3][2] + W[0][0] * W[3][1] * W[3][2] +
239 W[1][0] * W[2][0] * W[3][3] - W[0][0] * W[2][1] * W[3][3];
240
241 W[1][3] = -W[2][0] * W[3][0] * W[2][1] + W[2][0] * W[2][0] * W[3][1] + W[1][0] * W[3][0] * W[2][2] - W[0][0] * W[3][1] * W[2][2] -
242 W[1][0] * W[2][0] * W[3][2] + W[0][0] * W[2][1] * W[3][2];
243 W[2][3] = W[2][0] * W[3][0] * W[1][1] - W[1][0] * W[3][0] * W[2][1] - W[1][0] * W[2][0] * W[3][1] + W[0][0] * W[2][1] * W[3][1] +
244 W[1][0] * W[1][0] * W[3][2] - W[0][0] * W[1][1] * W[3][2];
245
246 double W00 = -W[3][1] * W[3][1] * W[2][2] + 2. * W[2][1] * W[3][1] * W[3][2] - W[1][1] * W[3][2] * W[3][2] -
247 W[2][1] * W[2][1] * W[3][3] + W[1][1] * W[2][2] * W[3][3];
248 double W11 = -W[3][0] * W[3][0] * W[2][2] + 2. * W[2][0] * W[3][0] * W[3][2] - W[0][0] * W[3][2] * W[3][2] -
249 W[2][0] * W[2][0] * W[3][3] + W[0][0] * W[2][2] * W[3][3];
250 double W22 = -W[3][0] * W[3][0] * W[1][1] + 2. * W[1][0] * W[3][0] * W[3][1] - W[0][0] * W[3][1] * W[3][1] -
251 W[1][0] * W[1][0] * W[3][3] + W[0][0] * W[1][1] * W[3][3];
252 double W33 = -W[2][0] * W[2][0] * W[1][1] + 2. * W[1][0] * W[2][0] * W[2][1] - W[0][0] * W[2][1] * W[2][1] -
253 W[1][0] * W[1][0] * W[2][2] + W[0][0] * W[1][1] * W[2][2];
254
255 for (int i = 0; i < 3; i++) {
256 for (int j = 1; j < 4; j++) {
257 if (i >= j) continue;
258 W[i][j] = W[i][j] / Determinant;
259 W[j][i] = W[i][j];
260 }
261 }
262 W[0][0] = W00 / Determinant;
263 W[1][1] = W11 / Determinant;
264 W[2][2] = W22 / Determinant;
265 W[3][3] = W33 / Determinant;
266}
267void CscBipolarStripFitter::Derivative(double A[][3], double fp[][1], double p0[][1], int imeas, const int *meas) const {
268 // calculate the derivatives and the 0th order approximation
269 // around the ADC samplings
270 double norm = p0[0][0];
271 // double parm[3]={1.,p0[1][0],p0[2][0]};
272 for (int i = 0; i < imeas; i++) {
273 int ii = meas[i];
274 double z = (ii - p0[1][0]) * m_tsampling / p0[2][0];
275 double repquant = 0.;
276 double dFdzNormalized = 0.;
277 if (z > 0.) {
278 repquant = FindPow(z) * std::exp(-z) / m_bipolarNormalization;
279 dFdzNormalized = repquant * (m_n / z + z / (m_n2 + 1) - (m_n + 1.) / (m_n2 + 1.) - 1.); // repquant*(m_n/z+z/(m_n+1)-2.);
280 }
281
282 A[ii][0] = repquant * (1. - z / (m_n2 + 1.)); // repquant*(1.-z/(m_n+1.));
283 // A[ii][0] = bipolar(&z,parm);
284 fp[ii][0] = norm * A[ii][0];
285
286 // double normOverZmax = norm/m_bipolarNormalization;
287 double commonpart = norm * dFdzNormalized; //(z,parm);
288 A[ii][1] = commonpart * (-m_tsampling / p0[2][0]);
289 A[ii][2] = commonpart * (-z / p0[2][0]);
290 }
291 // end of derivative/zeroth order calculations
292}
293
294int CscBipolarStripFitter::TheFitter(double *x, const double ex, const double *initValues, int imeas, int *meas, int ipar, int *par, double *chi2,
295 double *result) const {
296 // maximum iterations
297 const int maxIter = 7;
298 // tolerances
299 double fitTolerance0 = 0.1;
300 double fitTolerance1 = 0.01;
301 // CLHEP::HepMatrix p0(3,1,0); // the matrix of the initial fit parameters
302 double p0[3][1];
303 for (int j = 0; j < 3; j++) p0[j][0] = initValues[j];
304
305 // CLHEP::HepMatrix m(4,1,0); // the matrix of ADC measurements (samples: 0,1,2,3)
306 double m[4][1];
307 // CLHEP::HepMatrix W(4,4,0); // the error matrix of the ADC measurements
308 double W[4][4];
309 for (int i = 0; i < 4; i++) {
310 m[i][0] = x[i];
311 W[i][i] = ex * ex;
312 }
313 // covariances
314 W[0][1] = 0.03 * ex * ex;
315 W[0][2] = -0.411 * ex * ex;
316 W[0][3] = -0.188 * ex * ex;
317 W[1][2] = 0.0275 * ex * ex;
318 W[1][3] = -0.4303 * ex * ex;
319 W[2][3] = 0. * ex * ex;
320 W[1][0] = W[0][1];
321 W[2][0] = W[0][2];
322 W[3][0] = W[0][3];
323 W[2][1] = W[1][2];
324 W[3][1] = W[1][3];
325 W[3][2] = W[2][3];
326
327 // WW.invert(ierr);
328 InvertSymmetric4x4(W);
329
330 // Taylor std::expansion of the bipolar pulse model around the
331 // samplings : F(x) = F(p0) + A *(p-p0) + higher.order
332 // CLHEP::HepMatrix fp(4,1,0); // the matrix of 0th order approximation
333 double fp[4][1];
334 // CLHEP::HepMatrix A(4,3,0); // the matrix of derivatives
335 double A[4][3];
336 for (int i = 0; i < 4; ++i) {
337 fp[i][0] = 0.;
338 for (int j = 0; j < 3; ++j) { A[i][j] = 0.; }
339 }
340 // remarks :
341 // if the pulse peaks in the last sampling fit with a constant shaping time
342 // if the pulse peaks in the first sampling fit without using the last sampling
343 // (too large contribution to the chi2
344 int counter = 0;
345 bool converged = false;
346 double amplitudeChangeOld = 0.;
347 bool diverganceCandidate = false;
348 // CLHEP::HepMatrix weight(3,3,1); // weight matrix allocated once
349 // the non-fitted parts are taken care appropriately
350 // at least if the fitting parameters or measurements
351 // don't change during the fitting procedure
352 double weight[3][3];
353
354 // CLHEP::HepMatrix residFactor(3,4,0); // residFactor allocated once
355 double residFactor[3][4];
356
357 for (int i = 0; i < 3; ++i) {
358 for (int j = 0; j < 4; ++j) {
359 if (j < 3) { weight[i][j] = 0.; }
360 residFactor[i][j] = 0.;
361 }
362 }
363 weight[0][0] = 1.;
364 weight[1][1] = 1.;
365 weight[2][2] = 1.;
366
367 while (!converged && counter < maxIter) // fit loop
368 {
369 Derivative(A, fp, p0, imeas, meas); // calculate the matrix of derivatives and 0th order approximation
370 // matrix multiplication
371 // the weight matrix is symmetric
372 // weight= A.T()*W*A;//.assign(A.T()*W*A);
373
374 double helpmatrix[4][3];
375 for (int i = 0; i < 4; ++i) {
376 for (int j = 0; j < 3; ++j) { helpmatrix[i][j] = 0.; }
377 }
378
379 for (int i = 0; i < imeas; i++) {
380 int ii = meas[i];
381 for (int j = 0; j < ipar; j++) {
382 int jj = par[j];
383 for (int k = 0; k < imeas; k++) {
384 int kk = meas[k];
385 helpmatrix[ii][jj] += W[ii][kk] * A[kk][jj];
386 }
387 }
388 }
389 for (int i = 0; i < ipar; i++) {
390 int ii = par[i];
391 for (int j = i; j < ipar; j++) {
392 int jj = par[j];
393 weight[ii][jj] = 0.;
394 for (int k = 0; k < imeas; k++) {
395 int kk = meas[k];
396 weight[ii][jj] += A[kk][ii] * helpmatrix[kk][jj]; // A[kk][ii]*A[kk][jj];
397 }
398 // weight[ii][jj]*=W[0][0];
399 weight[jj][ii] = weight[ii][jj];
400 }
401 }
402 // weight.invert(ierr); // inversion of weight matrix
403 // hand-made inversion of 2x2 or 3x3 symmetric matrix
404 InvertMatrix(weight, ipar, par);
405
406 // calculate W*(A.T()*W)
407 // residFactor = weight*(A.T()*W);
408 double helpmatrix2[3][4];
409 for (int i = 0; i < 3; ++i) {
410 for (int j = 0; j < 4; ++j) { helpmatrix2[i][j] = 0.; }
411 }
412
413 for (int i = 0; i < ipar; i++) {
414 int ii = par[i];
415 for (int j = 0; j < imeas; j++) {
416 int jj = meas[j];
417 for (int k = 0; k < imeas; k++) {
418 int kk = meas[k];
419 helpmatrix2[ii][jj] += A[kk][ii] * W[kk][jj];
420 }
421 }
422 }
423
424 for (int i = 0; i < ipar; i++) {
425 int ii = par[i];
426 for (int j = 0; j < imeas; j++) {
427 int jj = meas[j];
428 residFactor[ii][jj] = 0.;
429 for (int k = 0; k < ipar; k++) {
430 int kk = par[k];
431 residFactor[ii][jj] += weight[ii][kk] * helpmatrix2[kk][jj];
432 }
433 // residFactor[ii][jj]*=W[0][0];
434 }
435 }
436
437 double paramDiff[3];
438 for (int i = 0; i < 3; ++i) { paramDiff[i] = 0.; }
439
440 for (int i = 0; i < ipar; i++) {
441 int ii = par[i];
442 // estimation of new parameters
443 // paramDiff[i][0] += (weight*(A.T()*W)*(m-fp))[i][0];
444 for (int j = 0; j < imeas; j++) {
445 int jj = meas[j];
446 paramDiff[ii] += residFactor[ii][jj] * (m[jj][0] - fp[jj][0]);
447 }
448 p0[ii][0] += paramDiff[ii];
449 }
450 std::cout << "##### " << p0[0][0] << " " << p0[1][0] << " " << p0[2][0] << std::endl;
451 // if the parameters are not physical, keep them sensible
452 // if peaking time less than -0.5
453 double peakingTime = p0[1][0] + m_zmax * p0[2][0] / m_tsampling;
454 if (peakingTime < -0.5 || peakingTime > 3.) p0[1][0] = initValues[1];
455
456 if (p0[0][0] < 0.) p0[0][0] = initValues[0];
457 double amplitudeChangeNew = std::abs(paramDiff[0]);
458 if (std::abs(paramDiff[0]) < fitTolerance0 && std::abs(paramDiff[1]) < fitTolerance1) {
459 converged = true;
460 // calculate chi2
461 // (m-fp).T()*W*(m-fp)
462 double residual[4];
463 for (int i = 0; i < 4; ++i) { residual[i] = 0.; }
464
465 for (int i = 0; i < imeas; i++) {
466 int ii = meas[i];
467 residual[i] = m[ii][0] - fp[ii][0];
468 }
469
470 double tmpChi2 = 0.;
471 double helpmatrixchi2[4][1];
472 for (int i = 0; i < 4; ++i) { helpmatrixchi2[i][0] = 0.; }
473 for (int i = 0; i < imeas; i++) {
474 int ii = meas[i];
475 for (int k = 0; k < imeas; k++) {
476 int kk = meas[k];
477 helpmatrixchi2[ii][0] += W[ii][kk] * residual[kk];
478 }
479 }
480 for (int k = 0; k < imeas; k++) {
481 int kk = meas[k];
482 tmpChi2 += residual[kk] * helpmatrixchi2[kk][0];
483 // std::cout << residual[kk] << " ";//*W[kk][kk]*residual[kk]<< " ";
484 }
485 // std::cout<<std::endl;
486 (*chi2) = tmpChi2;
487 } else if (counter > 4 && (amplitudeChangeNew > 4. * amplitudeChangeOld)) {
488 if (diverganceCandidate) {
489 // diverging fit
490 // return parabola interpolation
491 printf("%3.2f %3.2f %3.2f %3.2f\n ", x[0], x[1], x[2], x[3]);
492 printf("Diverging fit\n");
493 return 4;
494 } else
495 diverganceCandidate = true;
496 }
497 if (p0[0][0] < 0.) {
498 // negative amplitude
499 // fit diverged
500 // return parabola
501 return 4;
502 }
503 // if after a couple of iterations the amplitude is low
504 // reduce the tolerances (or the maximum iterations)
505 // low amplitude pulses tend to oscillate and exhaust all iterations
506 if (p0[0][0] < 20.) {
507 fitTolerance0 = 0.1;
508 fitTolerance1 = 0.05;
509 }
510 amplitudeChangeOld = amplitudeChangeNew;
511 counter++;
512 }
513 /*
514 std::cout << " Error matrix "<<std::endl;
515 for(int i=0;i<3;i++)
516 {
517 for(int j=0;j<3;j++)
518 std::cout << weight[i][j]<<"\t";
519 std::cout<<std::endl;
520 }
521 */
522
523 result[0] = p0[0][0];
524 result[1] = m_zmax * p0[2][0] / m_tsampling + p0[1][0];
525 result[2] = p0[2][0];
526
527 if (counter == maxIter) return 3;
528 return 0;
529}
530
531//**********************************************************************
ATH_CHECK_RECOVERABLE
#define ATH_CHECK_RECOVERABLE
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:48
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
ChargeList
ICscStripFitter::ChargeList ChargeList
Definition CalibCscStripFitter.cxx:14
res
std::pair< std::vector< unsigned int >, bool > res
Definition JetGroupProductTest.cxx:11
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
imax
int imax(int i, int j)
Definition TileLaserTimingTool.cxx:33
x
#define x
z
#define z
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
Definition AthCommonDataStore.h:95
CscBipolarStripFitter::CscBipolarStripFitter
CscBipolarStripFitter(const std::string &, const std::string &, const IInterface *)
Definition CscBipolarStripFitter.cxx:29
CscBipolarStripFitter::fit
Result fit(const ChargeList &charges, double samplingTime, Identifier &stripId) const
Definition CscBipolarStripFitter.cxx:74
CscBipolarStripFitter::Derivative
void Derivative(double A[][3], double fp[][1], double p0[][1], int imeas, const int *meas) const
Definition CscBipolarStripFitter.cxx:267
CscBipolarStripFitter::InvertMatrix
static void InvertMatrix(double matrix[][3], const int dim, const int *)
Definition CscBipolarStripFitter.cxx:175
CscBipolarStripFitter::FindPow
double FindPow(double z) const
Definition CscBipolarStripFitter.cxx:170
CscBipolarStripFitter::TheFitter
int TheFitter(double *x, const double ex, const double *initValues, int imeas, int *meas, int ipar, int *par, double *chi2, double *result) const
Definition CscBipolarStripFitter.cxx:294
CscBipolarStripFitter::m_cscCalibTool
ToolHandle< ICscCalibTool > m_cscCalibTool
Definition CscBipolarStripFitter.h:51
CscBipolarStripFitter::m_phelper
const CscIdHelper * m_phelper
Definition CscBipolarStripFitter.h:48
CscBipolarStripFitter::initialize
StatusCode initialize() override
Definition CscBipolarStripFitter.cxx:41
CscBipolarStripFitter::InvertSymmetric4x4
static void InvertSymmetric4x4(double W[][4])
Definition CscBipolarStripFitter.cxx:223
CscBipolarStripFitter::FindInitValues
static double FindInitValues(double *x, double *initValues, int *maxsample)
Definition CscBipolarStripFitter.cxx:129
ICscStripFitter::ChargeList
std::vector< float > ChargeList
Definition ICscStripFitter.h:55
IdentifierHash
This is a "hash" representation of an Identifier.
Definition IdentifierHash.h:25
Identifier::show
void show() const
Print out in hex form.
Definition Identifier.cxx:30
MuonGM::MuonDetectorManager
The MuonDetectorManager stores the transient representation of the Muon Spectrometer geometry and pro...
Definition MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonDetectorManager.h:51
Muon::CscStripPrepData
Class representing the raw data of one CSC strip (for clusters look at Muon::CscPrepData).
Definition CscStripPrepData.h:40
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
A
hold the test vectors and ease the comparison
type
result
Generated on for ATLAS Offline Software by doxygen 1.14.0