ATLAS Offline Software
Loading...
Searching...
No Matches
CscBipolarStripFitter.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 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 \n" << stripId);
83 }
84 ATH_MSG_DEBUG("CalibCscStripFitter:: " << stripId << " " << (unsigned int)stripHash);
85
90 double initValues[3];
91 for (int i = 0; i < 3; ++i) { initValues[i] = 0.; }
92 // int imax = -1;
93 initValues[2] = 8.; // m_width
95
96 double noise = m_cscCalibTool->stripNoise(stripHash);
97 res.status = 0;
98 if (m_qerr > 0)
99 res.dcharge = m_qerr;
100 else
101 res.dcharge = noise;
102
103 res.dtime = m_terr;
104 res.charge = initValues[0];
105 res.time = initValues[1];
106
107 initValues[1] = res.time - 2.5;
108 // int imeas =4;
109 // int meas[4] = {true,true,true,true};
110 // int ipar = 3;
111 // int par[3] = {true,true,false};
112 // double m_chi2;
113 //double result[3];
115
116 // Add an error proportional to the charge.
117 double dqprop = m_qerrprop * res.charge;
118 res.dcharge = sqrt(res.dcharge * res.dcharge + dqprop * dqprop);
119 // Display result.
120 ATH_MSG_DEBUG(" Status: " << res.status);
121 ATH_MSG_DEBUG(" Charge: " << res.charge);
122 ATH_MSG_DEBUG(" Charge err: " << res.dcharge);
123 ATH_MSG_DEBUG(" Time: " << res.time);
124 //ATH_MSG_DEBUG(" fit : " << result[0] << " " << result[1] << " " << result[2]);
125 return res;
126}
127
128double CscBipolarStripFitter::FindInitValues(double *x, double *initValues, int *maxsample) {
129 // find maximum sample imax:
130 double peakingTime = -99.; // interpolated peaking time in samples
131 double amplitude = -99.; // interpolated amplitude
132 double amin, amax;
133 int imax = 0;
134 const int numSample = 4;
135 amax = x[0];
136 amin = x[0];
137 for (int j = 1; j < numSample; j++) {
138 if (amax < x[j]) {
139 amax = x[j];
140 imax = j;
141 }
142 if (amin > x[j]) amin = x[j];
143 }
144
145 // calculate peak and amplitude:
146 (*maxsample) = imax;
147 if ((imax == 0) || (imax == numSample - 1)) // no interpolation possible
148 {
149 peakingTime = imax;
150 amplitude = amax;
151 } else // do parabolic interpolation
152 {
153 double a, b, c; // coeffients of parabola y=a*x*x+b*x+c
154 a = 0.5 * (x[imax + 1] + x[imax - 1] - 2.0 * x[imax]);
155 b = 0.5 * (x[imax + 1] - x[imax - 1]);
156 c = x[imax];
157
158 peakingTime = -0.5 * b / a;
159 amplitude = a * peakingTime * peakingTime + b * peakingTime + c;
160 peakingTime += imax; // it was relative to imax
161 }
162
163 initValues[0] = amplitude;
164 initValues[1] = peakingTime;
165 // - m_zmax*initValues[2]/m_tsampling;
166 return x[imax];
167}
168
169double CscBipolarStripFitter::FindPow(double z) const {
170 double zpower = std::exp(m_n * std::log(z));
171 return zpower;
172}
173
174void CscBipolarStripFitter::InvertMatrix(double matrix[][3], const int dim, const int *correspdim) {
175 // invert 2x2 or 3x3 symmetric matrix
176 if (dim == 1) {
177 int ii = correspdim[0];
178 matrix[ii][ii] = 1. / matrix[ii][ii];
179
180 } else if (dim == 2) {
181 int ii = correspdim[0];
182 int jj = correspdim[1];
183 double determinant = -matrix[jj][ii] * matrix[ii][jj] + matrix[ii][ii] * matrix[jj][jj];
184 // if(std::abs(determinant)<1.e-13)
185 // std::cout<<" zero determinant "<<std::endl;
186 double i00 = matrix[ii][ii];
187 matrix[ii][ii] = matrix[jj][jj] / determinant;
188 matrix[jj][jj] = i00 / determinant;
189 matrix[ii][jj] = -matrix[ii][jj] / determinant;
190 matrix[jj][ii] = matrix[ii][jj];
191
192 } else if (dim == 3) {
193 double sm12 = matrix[0][1] * matrix[0][1];
194 double sm13 = matrix[0][2] * matrix[0][2];
195 double sm23 = matrix[1][2] * matrix[1][2];
196 double determinant = matrix[0][0] * matrix[1][1] * matrix[2][2] - matrix[0][0] * sm23 - sm12 * matrix[2][2] +
197 2. * matrix[0][1] * matrix[0][2] * matrix[1][2] - matrix[1][1] * sm13;
198
199 // if(std::abs(determinant)<1.e-13)
200 // std::cout << "zero determinant"<<std::endl;
201 double i00 = matrix[1][1] * matrix[2][2] - sm23;
202 double i11 = matrix[0][0] * matrix[2][2] - sm13;
203 double i22 = matrix[0][0] * matrix[1][1] - sm12;
204
205 matrix[1][0] = (matrix[0][2] * matrix[1][2] - matrix[2][2] * matrix[0][1]) / determinant;
206 matrix[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1]) / determinant;
207 matrix[2][1] = (matrix[0][1] * matrix[0][2] - matrix[0][0] * matrix[1][2]) / determinant;
208 matrix[0][1] = matrix[1][0];
209 matrix[0][2] = matrix[2][0];
210 matrix[1][2] = matrix[2][1];
211 matrix[0][0] = i00 / determinant;
212 matrix[1][1] = i11 / determinant;
213 matrix[2][2] = i22 / determinant;
214
215 } else {
216 // int ierr;
217 // matrix.invert(ierr);
218 printf("this is not a 1x1 or 2x2 or 3x3 matrix\n");
219 }
220}
221
222void CscBipolarStripFitter::InvertSymmetric4x4(double W[][4]) {
223 double Determinant =
224 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] -
225 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] +
226 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] -
227 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] -
228 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] -
229 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];
230
231 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] +
232 W[2][0] * W[2][1] * W[3][3] - W[1][0] * W[2][2] * W[3][3];
233 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] -
234 W[2][0] * W[1][1] * W[3][3] + W[1][0] * W[2][1] * W[3][3];
235 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] +
236 W[2][0] * W[1][1] * W[3][2] - W[1][0] * W[2][1] * W[3][2];
237 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] +
238 W[1][0] * W[2][0] * W[3][3] - W[0][0] * W[2][1] * W[3][3];
239
240 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] -
241 W[1][0] * W[2][0] * W[3][2] + W[0][0] * W[2][1] * W[3][2];
242 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] +
243 W[1][0] * W[1][0] * W[3][2] - W[0][0] * W[1][1] * W[3][2];
244
245 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] -
246 W[2][1] * W[2][1] * W[3][3] + W[1][1] * W[2][2] * W[3][3];
247 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] -
248 W[2][0] * W[2][0] * W[3][3] + W[0][0] * W[2][2] * W[3][3];
249 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] -
250 W[1][0] * W[1][0] * W[3][3] + W[0][0] * W[1][1] * W[3][3];
251 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] -
252 W[1][0] * W[1][0] * W[2][2] + W[0][0] * W[1][1] * W[2][2];
253
254 for (int i = 0; i < 3; i++) {
255 for (int j = 1; j < 4; j++) {
256 if (i >= j) continue;
257 W[i][j] = W[i][j] / Determinant;
258 W[j][i] = W[i][j];
259 }
260 }
261 W[0][0] = W00 / Determinant;
262 W[1][1] = W11 / Determinant;
263 W[2][2] = W22 / Determinant;
264 W[3][3] = W33 / Determinant;
265}
266void CscBipolarStripFitter::Derivative(double A[][3], double fp[][1], double p0[][1], int imeas, const int *meas) const {
267 // calculate the derivatives and the 0th order approximation
268 // around the ADC samplings
269 double norm = p0[0][0];
270 // double parm[3]={1.,p0[1][0],p0[2][0]};
271 for (int i = 0; i < imeas; i++) {
272 int ii = meas[i];
273 double z = (ii - p0[1][0]) * m_tsampling / p0[2][0];
274 double repquant = 0.;
275 double dFdzNormalized = 0.;
276 if (z > 0.) {
277 repquant = FindPow(z) * std::exp(-z) / m_bipolarNormalization;
278 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.);
279 }
280
281 A[ii][0] = repquant * (1. - z / (m_n2 + 1.)); // repquant*(1.-z/(m_n+1.));
282 // A[ii][0] = bipolar(&z,parm);
283 fp[ii][0] = norm * A[ii][0];
284
285 // double normOverZmax = norm/m_bipolarNormalization;
286 double commonpart = norm * dFdzNormalized; //(z,parm);
287 A[ii][1] = commonpart * (-m_tsampling / p0[2][0]);
288 A[ii][2] = commonpart * (-z / p0[2][0]);
289 }
290 // end of derivative/zeroth order calculations
291}
292
293int CscBipolarStripFitter::TheFitter(double *x, const double ex, const double *initValues, int imeas, int *meas, int ipar, int *par, double *chi2,
294 double *result) const {
295 // maximum iterations
296 const int maxIter = 7;
297 // tolerances
298 double fitTolerance0 = 0.1;
299 double fitTolerance1 = 0.01;
300 // CLHEP::HepMatrix p0(3,1,0); // the matrix of the initial fit parameters
301 double p0[3][1];
302 for (int j = 0; j < 3; j++) p0[j][0] = initValues[j];
303
304 // CLHEP::HepMatrix m(4,1,0); // the matrix of ADC measurements (samples: 0,1,2,3)
305 double m[4][1];
306 // CLHEP::HepMatrix W(4,4,0); // the error matrix of the ADC measurements
307 double W[4][4];
308 for (int i = 0; i < 4; i++) {
309 m[i][0] = x[i];
310 W[i][i] = ex * ex;
311 }
312 // covariances
313 W[0][1] = 0.03 * ex * ex;
314 W[0][2] = -0.411 * ex * ex;
315 W[0][3] = -0.188 * ex * ex;
316 W[1][2] = 0.0275 * ex * ex;
317 W[1][3] = -0.4303 * ex * ex;
318 W[2][3] = 0. * ex * ex;
319 W[1][0] = W[0][1];
320 W[2][0] = W[0][2];
321 W[3][0] = W[0][3];
322 W[2][1] = W[1][2];
323 W[3][1] = W[1][3];
324 W[3][2] = W[2][3];
325
326 // WW.invert(ierr);
327 InvertSymmetric4x4(W);
328
329 // Taylor std::expansion of the bipolar pulse model around the
330 // samplings : F(x) = F(p0) + A *(p-p0) + higher.order
331 // CLHEP::HepMatrix fp(4,1,0); // the matrix of 0th order approximation
332 double fp[4][1];
333 // CLHEP::HepMatrix A(4,3,0); // the matrix of derivatives
334 double A[4][3];
335 for (int i = 0; i < 4; ++i) {
336 fp[i][0] = 0.;
337 for (int j = 0; j < 3; ++j) { A[i][j] = 0.; }
338 }
339 // remarks :
340 // if the pulse peaks in the last sampling fit with a constant shaping time
341 // if the pulse peaks in the first sampling fit without using the last sampling
342 // (too large contribution to the chi2
343 int counter = 0;
344 bool converged = false;
345 double amplitudeChangeOld = 0.;
346 bool diverganceCandidate = false;
347 // CLHEP::HepMatrix weight(3,3,1); // weight matrix allocated once
348 // the non-fitted parts are taken care appropriately
349 // at least if the fitting parameters or measurements
350 // don't change during the fitting procedure
351 double weight[3][3];
352
353 // CLHEP::HepMatrix residFactor(3,4,0); // residFactor allocated once
354 double residFactor[3][4];
355
356 for (int i = 0; i < 3; ++i) {
357 for (int j = 0; j < 4; ++j) {
358 if (j < 3) { weight[i][j] = 0.; }
359 residFactor[i][j] = 0.;
360 }
361 }
362 weight[0][0] = 1.;
363 weight[1][1] = 1.;
364 weight[2][2] = 1.;
365
366 while (!converged && counter < maxIter) // fit loop
367 {
368 Derivative(A, fp, p0, imeas, meas); // calculate the matrix of derivatives and 0th order approximation
369 // matrix multiplication
370 // the weight matrix is symmetric
371 // weight= A.T()*W*A;//.assign(A.T()*W*A);
372
373 double helpmatrix[4][3];
374 for (int i = 0; i < 4; ++i) {
375 for (int j = 0; j < 3; ++j) { helpmatrix[i][j] = 0.; }
376 }
377
378 for (int i = 0; i < imeas; i++) {
379 int ii = meas[i];
380 for (int j = 0; j < ipar; j++) {
381 int jj = par[j];
382 for (int k = 0; k < imeas; k++) {
383 int kk = meas[k];
384 helpmatrix[ii][jj] += W[ii][kk] * A[kk][jj];
385 }
386 }
387 }
388 for (int i = 0; i < ipar; i++) {
389 int ii = par[i];
390 for (int j = i; j < ipar; j++) {
391 int jj = par[j];
392 weight[ii][jj] = 0.;
393 for (int k = 0; k < imeas; k++) {
394 int kk = meas[k];
395 weight[ii][jj] += A[kk][ii] * helpmatrix[kk][jj]; // A[kk][ii]*A[kk][jj];
396 }
397 // weight[ii][jj]*=W[0][0];
398 weight[jj][ii] = weight[ii][jj];
399 }
400 }
401 // weight.invert(ierr); // inversion of weight matrix
402 // hand-made inversion of 2x2 or 3x3 symmetric matrix
403 InvertMatrix(weight, ipar, par);
404
405 // calculate W*(A.T()*W)
406 // residFactor = weight*(A.T()*W);
407 double helpmatrix2[3][4];
408 for (int i = 0; i < 3; ++i) {
409 for (int j = 0; j < 4; ++j) { helpmatrix2[i][j] = 0.; }
410 }
411
412 for (int i = 0; i < ipar; i++) {
413 int ii = par[i];
414 for (int j = 0; j < imeas; j++) {
415 int jj = meas[j];
416 for (int k = 0; k < imeas; k++) {
417 int kk = meas[k];
418 helpmatrix2[ii][jj] += A[kk][ii] * W[kk][jj];
419 }
420 }
421 }
422
423 for (int i = 0; i < ipar; i++) {
424 int ii = par[i];
425 for (int j = 0; j < imeas; j++) {
426 int jj = meas[j];
427 residFactor[ii][jj] = 0.;
428 for (int k = 0; k < ipar; k++) {
429 int kk = par[k];
430 residFactor[ii][jj] += weight[ii][kk] * helpmatrix2[kk][jj];
431 }
432 // residFactor[ii][jj]*=W[0][0];
433 }
434 }
435
436 double paramDiff[3];
437 for (int i = 0; i < 3; ++i) { paramDiff[i] = 0.; }
438
439 for (int i = 0; i < ipar; i++) {
440 int ii = par[i];
441 // estimation of new parameters
442 // paramDiff[i][0] += (weight*(A.T()*W)*(m-fp))[i][0];
443 for (int j = 0; j < imeas; j++) {
444 int jj = meas[j];
445 paramDiff[ii] += residFactor[ii][jj] * (m[jj][0] - fp[jj][0]);
446 }
447 p0[ii][0] += paramDiff[ii];
448 }
449 std::cout << "##### " << p0[0][0] << " " << p0[1][0] << " " << p0[2][0] << std::endl;
450 // if the parameters are not physical, keep them sensible
451 // if peaking time less than -0.5
452 double peakingTime = p0[1][0] + m_zmax * p0[2][0] / m_tsampling;
453 if (peakingTime < -0.5 || peakingTime > 3.) p0[1][0] = initValues[1];
454
455 if (p0[0][0] < 0.) p0[0][0] = initValues[0];
456 double amplitudeChangeNew = std::abs(paramDiff[0]);
457 if (std::abs(paramDiff[0]) < fitTolerance0 && std::abs(paramDiff[1]) < fitTolerance1) {
458 converged = true;
459 // calculate chi2
460 // (m-fp).T()*W*(m-fp)
461 double residual[4];
462 for (int i = 0; i < 4; ++i) { residual[i] = 0.; }
463
464 for (int i = 0; i < imeas; i++) {
465 int ii = meas[i];
466 residual[i] = m[ii][0] - fp[ii][0];
467 }
468
469 double tmpChi2 = 0.;
470 double helpmatrixchi2[4][1];
471 for (int i = 0; i < 4; ++i) { helpmatrixchi2[i][0] = 0.; }
472 for (int i = 0; i < imeas; i++) {
473 int ii = meas[i];
474 for (int k = 0; k < imeas; k++) {
475 int kk = meas[k];
476 helpmatrixchi2[ii][0] += W[ii][kk] * residual[kk];
477 }
478 }
479 for (int k = 0; k < imeas; k++) {
480 int kk = meas[k];
481 tmpChi2 += residual[kk] * helpmatrixchi2[kk][0];
482 // std::cout << residual[kk] << " ";//*W[kk][kk]*residual[kk]<< " ";
483 }
484 // std::cout<<std::endl;
485 (*chi2) = tmpChi2;
486 } else if (counter > 4 && (amplitudeChangeNew > 4. * amplitudeChangeOld)) {
487 if (diverganceCandidate) {
488 // diverging fit
489 // return parabola interpolation
490 printf("%3.2f %3.2f %3.2f %3.2f\n ", x[0], x[1], x[2], x[3]);
491 printf("Diverging fit\n");
492 return 4;
493 } else
494 diverganceCandidate = true;
495 }
496 if (p0[0][0] < 0.) {
497 // negative amplitude
498 // fit diverged
499 // return parabola
500 return 4;
501 }
502 // if after a couple of iterations the amplitude is low
503 // reduce the tolerances (or the maximum iterations)
504 // low amplitude pulses tend to oscillate and exhaust all iterations
505 if (p0[0][0] < 20.) {
506 fitTolerance0 = 0.1;
507 fitTolerance1 = 0.05;
508 }
509 amplitudeChangeOld = amplitudeChangeNew;
510 counter++;
511 }
512 /*
513 std::cout << " Error matrix "<<std::endl;
514 for(int i=0;i<3;i++)
515 {
516 for(int j=0;j<3;j++)
517 std::cout << weight[i][j]<<"\t";
518 std::cout<<std::endl;
519 }
520 */
521
522 result[0] = p0[0][0];
523 result[1] = m_zmax * p0[2][0] / m_tsampling + p0[1][0];
524 result[2] = p0[2][0];
525
526 if (counter == maxIter) return 3;
527 return 0;
528}
529
530//**********************************************************************
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:266
CscBipolarStripFitter::InvertMatrix
static void InvertMatrix(double matrix[][3], const int dim, const int *)
Definition CscBipolarStripFitter.cxx:174
CscBipolarStripFitter::FindPow
double FindPow(double z) const
Definition CscBipolarStripFitter.cxx:169
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:293
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:222
CscBipolarStripFitter::FindInitValues
static double FindInitValues(double *x, double *initValues, int *maxsample)
Definition CscBipolarStripFitter.cxx:128
ICscStripFitter::ChargeList
std::vector< float > ChargeList
Definition ICscStripFitter.h:55
IdentifierHash
This is a "hash" representation of an Identifier.
Definition IdentifierHash.h:25
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
Generated on for ATLAS Offline Software by doxygen 1.14.0