44 double zpower =
z*
z*
z;
57 double peakingTime=-99.;
58 double amplitude=-99.;
61 const int numSample=4;
64 for(
int j=1; j<numSample; j++)
89 peakingTime = -0.5*
b/
a;
90 amplitude =
a*peakingTime*peakingTime +
b*peakingTime +
c;
94 initValues[0] = amplitude;
129 for(
int i=0;
i<imeas;
i++)
133 double repquant = 0.;
134 double dFdzNormalized = 0.;
138 dFdzNormalized= repquant*(
m_n/
z+
z/13.-2.);
141 A[ii][0] = repquant*(1.-
z/13.);
146 double commonpart =
norm* dFdzNormalized;
148 A[ii][2] = commonpart * (-
z/
p0[2][0]);
159 int ii=correspdim[0];
160 int jj=correspdim[1];
162 double i00 =
matrix[ii][ii];
164 matrix[jj][jj] = i00/determinant;
189 matrix[0][0] = i00/determinant;
190 matrix[1][1] = i11/determinant;
191 matrix[2][2] = i22/determinant;
198 printf(
"this is not a 2x2 or 3x3 matrix\n");
236 double initValues[3]={0.};
238 initValues[2]=predefinedwidth;
240 result[0] = initValues[0];
241 result[1] = initValues[1];
242 result[2] = initValues[2];
245 if(samplemax < 5.*ex)
248 bool onesaturated =
false;
264 if(
x[
imax]+pedestal>4000.)
266 if(
x[
imax-1]+pedestal>4000. ||
x[
imax+1]+pedestal>4000.)
275 bool fitpar[3] = {
true,
true,
false};
276 bool usemeas[4] = {
true,
true,
true,
true};
290 usemeas[
imax]=
false;
314 result[0] = initValues[0];
328 double Determinant =
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]
329 +
W[0][2]*
W[0][2]*
W[1][3]*
W[1][3]-
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]+2.*
W[0][2]*
W[0][3]*
W[1][1]*
W[2][3]
330 - 2.*
W[0][1]*
W[0][3]*
W[1][2]*
W[2][3] - 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]-
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]-
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];
332 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]+
W[2][0]*
W[2][1]*
W[3][3]-
W[1][0]*
W[2][2]*
W[3][3];
333 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]-
W[2][0]*
W[1][1]*
W[3][3]+
W[1][0]*
W[2][1]*
W[3][3];
334 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]+
W[2][0]*
W[1][1]*
W[3][2]-
W[1][0]*
W[2][1]*
W[3][2];
335 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]+
W[1][0]*
W[2][0]*
W[3][3]-
W[0][0]*
W[2][1]*
W[3][3];
337 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]-
W[1][0]*
W[2][0]*
W[3][2]+
W[0][0]*
W[2][1]*
W[3][2];
338 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]+
W[1][0]*
W[1][0]*
W[3][2]-
W[0][0]*
W[1][1]*
W[3][2];
340 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]-
W[2][1]*
W[2][1]*
W[3][3]+
W[1][1]*
W[2][2]*
W[3][3];
341 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]-
W[2][0]*
W[2][0]*
W[3][3]+
W[0][0]*
W[2][2]*
W[3][3];
342 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]-
W[1][0]*
W[1][0]*
W[3][3]+
W[0][0]*
W[1][1]*
W[3][3];
343 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]-
W[1][0]*
W[1][0]*
W[2][2]+
W[0][0]*
W[1][1]*
W[2][2];
351 W[
i][j] =
W[
i][j]/Determinant;
355 W[0][0] = W00/Determinant;
356 W[1][1] = W11/Determinant;
357 W[2][2] = W22/Determinant;
358 W[3][3] = W33/Determinant;
365 const int maxIter = 7;
367 double fitTolerance0 = 0.10;
368 double fitTolerance1 = 0.01;
370 double p0[3][1]={{0.}};
372 p0[j][0] = initValues[j];
375 double m[4][1]={{0.}};
377 double W[4][4]={{0.}};
384 W[0][1] = 0.03*ex*ex;
385 W[0][2] = -0.411*ex*ex;
386 W[0][3] = -0.188*ex*ex;
387 W[1][2] = 0.0275*ex*ex;
388 W[1][3] = -0.4303*ex*ex;
403 double fp[4][1]={{0.}};
405 double A[4][3]={{0.}};
411 bool converged=
false;
412 double amplitudeChangeOld = 0.;
413 bool diverganceCandidate =
false;
418 double weight[3][3]={{0.}};
424 double residFactor[3][4]={{0.}};
425 while(!converged &&
counter<maxIter)
432 double helpmatrix[4][3]={{0.}};
433 for(
int i=0;
i<imeas;
i++)
436 for(
int j=0;j<ipar;j++)
439 for(
int k=0;
k<imeas;
k++)
442 helpmatrix[ii][jj] +=
W[ii][
kk]*
A[
kk][jj];
446 for(
int i=0;
i<ipar;
i++)
449 for(
int j=
i;j<ipar;j++)
453 for(
int k=0;
k<imeas;
k++)
468 double helpmatrix2[3][4]={{0.}};
469 for(
int i=0;
i<ipar;
i++)
472 for(
int j=0;j<imeas;j++)
475 for(
int k=0;
k<imeas;
k++)
478 helpmatrix2[ii][jj] +=
A[
kk][ii] *
W[
kk][jj];
483 for(
int i=0;
i<ipar;
i++)
486 for(
int j=0;j<imeas;j++)
489 residFactor[ii][jj]=0.;
490 for(
int k=0;
k<ipar;
k++)
493 residFactor[ii][jj] +=
weight[ii][
kk]*helpmatrix2[
kk][jj];
499 double paramDiff[3] = {0.};
500 for(
int i=0;
i<ipar;
i++)
505 for(
int j=0;j<imeas;j++)
508 paramDiff[ii] += residFactor[ii][jj]*(
m[jj][0]-
fp[jj][0]);
510 p0[ii][0] += paramDiff[ii];
513 if(
p0[1][0]>1. ||
p0[1][0]<-3.)
514 p0[1][0] = initValues[1];
516 double amplitudeChangeNew = std::abs(paramDiff[0]);
517 if(std::abs(paramDiff[0])<fitTolerance0 && std::abs(paramDiff[1])<fitTolerance1)
522 for(
int i=0;
i<imeas;
i++)
529 double helpmatrixchi2[4][1]={{0.}};
530 for(
int i=0;
i<imeas;
i++)
533 for(
int k=0;
k<imeas;
k++)
539 for(
int k=0;
k<imeas;
k++)
546 else if(
counter>4 && (amplitudeChangeNew>4.*amplitudeChangeOld))
548 if(diverganceCandidate)
552 printf(
"Diverging fit\n");
556 diverganceCandidate =
true;
571 fitTolerance1 = 0.05;
573 amplitudeChangeOld = amplitudeChangeNew;