19 #include <ecl/modules/eclWaveformFit/ECLWaveformFit.h>
20 #include <ecl/dbobjects/ECLCrystalCalib.h>
21 #include <ecl/digitization/EclConfiguration.h>
22 #include <ecl/digitization/shaperdsp.h>
23 #include <ecl/dbobjects/ECLDigitWaveformParameters.h>
24 #include <ecl/dbobjects/ECLDigitWaveformParametersForMC.h>
25 #include <ecl/dbobjects/ECLAutoCovariance.h>
28 #include <framework/core/Environment.h>
29 #include <framework/database/DBObjPtr.h>
34 #include <TMatrixDSym.h>
35 #include <TDecompChol.h>
54 std::vector<double> fitA(31);
62 std::vector< std::vector<double> > currentCovMat;
66 void FCN2h(
int&,
double* grad,
double& f,
double* p,
int)
69 std::vector<double> df(N);
70 std::vector<double> da(N);
71 const double Ag = p[1], B = p[0], T = p[2], Ah = p[3];
72 double chi2 = 0, gAg = 0, gB = 0, gT = 0, gAh = 0;
80 for (
int i = 0; i < N; ++i) df[i] = fitA[i] - (Ag * ADg[i].f0 + Ah * ADh[i].f0 + B);
83 for (
int i = 0; i < N; ++i) da[i] = std::inner_product(currentCovMat[i].begin(), currentCovMat[i].end(), df.begin(), 0.0);
85 for (
int i = 0; i < N; ++i) {
86 chi2 += da[i] * df[i];
88 gAg -= da[i] * ADg[i].
f0;
89 gT -= da[i] * (ADg[i].
f1 * Ag + ADh[i].
f1 * Ah);
90 gAh -= da[i] * ADh[i].
f0;
101 void FCN2h2(
int&,
double* grad,
double& f,
double* p,
int)
104 std::vector<double> df(N);
105 std::vector<double> da(N);
106 const double A2 = p[4], T2 = p[5];
107 const double Ag = p[1], B = p[0], T = p[2], Ah = p[3];
108 double chi2 = 0, gA2 = 0, gT2 = 0;
109 double gAg = 0, gB = 0, gT = 0, gAh = 0;
118 for (
int i = 0; i < N; ++i) df[i] = fitA[i] - (Ag * ADg[i].f0 + Ah * ADh[i].f0 + A2 * AD2[i].f0 + B);
121 for (
int i = 0; i < N; ++i) da[i] = std::inner_product(currentCovMat[i].begin(), currentCovMat[i].end(), df.begin(), 0.0);
123 for (
int i = 0; i < N; ++i) {
124 chi2 += da[i] * df[i];
126 gAg -= da[i] * ADg[i].
f0;
127 gAh -= da[i] * ADh[i].
f0;
128 gT -= da[i] * (ADg[i].
f1 * Ag + ADh[i].
f1 * Ah);;
130 gA2 -= da[i] * AD2[i].
f0;
131 gT2 -= da[i] * AD2[i].
f1 * A2;
144 void regularize(
double* dst,
const double* src,
const int n,
const double u0 = 13.0,
const double u1 = 0.8)
146 for (
int k = 0; k < n; k++) dst[k] = src[k] / (1 + exp((k - u0) / u1));
150 bool makecovariance(
CovariancePacked& M,
const int nnoise,
const double* acov)
154 for (
int i = 0; i < ns; i++)
155 for (
int j = 0; j < i + 1; j++)
156 if (i - j < nnoise) E(i, j) = E(j, i) = acov[i - j];
159 const bool status = dc.Invert(E);
163 for (
int i = 0; i < ns; i++)
164 for (
int j = 0; j < i + 1; j++)
165 M[count++] = E(i, j);
166 M.
sigma = sqrtf(acov[0]);
177 currentCovMat.clear();
178 currentCovMat.resize(ns);
179 for (
int i = 0; i < ns; i++) {
180 currentCovMat[i].resize(ns);
181 for (
int j = 0; j < i + 1; j++) {
182 currentCovMat[i][j] = currentCovMat[j][i] = matrixPacked[count++];
185 aNoise = matrixPacked.
sigma;
194 setDescription(
"Module to fit offline waveforms and measure hadron scintillation component light output.");
195 setPropertyFlags(c_ParallelProcessingCertified);
196 addParam(
"EnergyThreshold", m_EnergyThreshold,
"Energy threshold of online fit result for Fitting Waveforms (GeV).", 0.03);
197 addParam(
"Chi2Threshold25dof", m_chi2Threshold25dof,
"chi2 threshold (25 dof) to classify offline fit as good fit.", 57.1);
198 addParam(
"Chi2Threshold27dof", m_chi2Threshold27dof,
"chi2 threshold (27 dof) to classify offline fit as good fit.", 60.0);
199 addParam(
"CovarianceMatrix", m_CovarianceMatrix,
200 "Option to use crystal dependent covariance matrices (false uses identity matrix).",
true);
212 m_TemplatesLoaded =
true;
214 if (m_IsMCFlag == 0) {
217 std::vector<double> Ptemp(11), Htemp(11), Dtemp(11);
218 for (
int i = 0; i < 8736; i++) {
219 for (
int j = 0; j < 11; j++) {
220 Ptemp[j] = (double)WavePars->getPhotonParameters(i + 1)[j];
221 Htemp[j] = (double)WavePars->getHadronParameters(i + 1)[j];
222 Dtemp[j] = (double)WavePars->getDiodeParameters(i + 1)[j];
231 std::vector<double> Ptemp(11), Htemp(11), Dtemp(11);
232 for (
int j = 0; j < 11; j++) {
233 Ptemp[j] = (double)WaveParsMC->getPhotonParameters()[j];
234 Htemp[j] = (double)WaveParsMC->getHadronParameters()[j];
235 Dtemp[j] = (double)WaveParsMC->getDiodeParameters()[j];
247 m_TemplatesLoaded =
false;
250 m_ADCtoEnergy.resize(8736);
251 if (Ael)
for (
int i = 0; i < 8736; i++) m_ADCtoEnergy[i] = Ael->getCalibVector()[i];
252 if (Aen)
for (
int i = 0; i < 8736; i++) m_ADCtoEnergy[i] *= Aen->getCalibVector()[i];
255 if (m_CovarianceMatrix) {
257 for (
int id = 1;
id <= 8736;
id++) {
258 constexpr
int N = 31;
259 std::vector<double> buf(N);
260 std::vector<double> reg(N);
261 cov->getAutoCovariance(
id, buf.data());
264 while (!makecovariance(m_c[
id - 1], N, reg.data()))
265 regularize(buf.data(), reg.data(), N, x0 -= 1, 1);
269 const double isigma = 1 / 7.5;
270 currentCovMat.clear();
271 currentCovMat.resize(31);
272 for (
int i = 0; i < 31; ++i) {
273 currentCovMat[i].resize(31);
274 for (
int j = 0; j < 31; ++j) {
275 currentCovMat[i][j] = (i == j) * isigma * isigma;
285 m_eclDSPs.registerInDataStore();
286 m_eclDigits.registerInDataStore();
289 m_Minit2h =
new TMinuit(4);
290 m_Minit2h->SetFCN(FCN2h);
294 m_Minit2h->mnexcm(
"SET PRIntout", arglist, 1, ierflg);
295 m_Minit2h->mnexcm(
"SET NOWarnings", arglist, 0, ierflg);
297 m_Minit2h->mnexcm(
"SET ERR", arglist, 1, ierflg);
299 m_Minit2h->mnexcm(
"SET STRategy", arglist, 1, ierflg);
301 m_Minit2h->mnexcm(
"SET GRAdient", arglist, 1, ierflg);
303 m_Minit2h->mnexcm(
"SET EPSmachine", arglist, 1, ierflg);
306 m_Minit2h2 =
new TMinuit(6);
307 m_Minit2h2->SetFCN(FCN2h2);
309 m_Minit2h2->mnexcm(
"SET PRIntout", arglist, 1, ierflg);
310 m_Minit2h2->mnexcm(
"SET NOWarnings", arglist, 0, ierflg);
312 m_Minit2h2->mnexcm(
"SET ERR", arglist, 1, ierflg);
314 m_Minit2h2->mnexcm(
"SET STRategy", arglist, 1, ierflg);
316 m_Minit2h2->mnexcm(
"SET GRAdient", arglist, 1, ierflg);
318 m_Minit2h2->mnexcm(
"SET EPSmachine", arglist, 1, ierflg);
321 m_TemplatesLoaded =
false;
328 if (!m_TemplatesLoaded) {
330 if (m_eclDSPs.getEntries() > 0) loadTemplateParameterArray();
333 for (
auto& aECLDsp : m_eclDSPs) {
335 aECLDsp.setTwoComponentTotalAmp(-1);
336 aECLDsp.setTwoComponentHadronAmp(-1);
337 aECLDsp.setTwoComponentDiodeAmp(-1);
338 aECLDsp.setTwoComponentChi2(-1);
342 aECLDsp.setTwoComponentTime(-1);
343 aECLDsp.setTwoComponentBaseline(-1);
346 const int id = aECLDsp.getCellId() - 1;
349 for (
int j = 0; j < ec.
m_nsmp; j++) fitA[j] = aECLDsp.getDspA()[j];
353 for (
const auto& aECLDigit : m_eclDigits) {
354 if (aECLDigit.getCellId() - 1 ==
id) {
356 aECLDsp.addRelationTo(&aECLDigit);
360 if (d == NULL)
continue;
363 if (d->getAmp() * m_ADCtoEnergy[
id] < m_EnergyThreshold)
continue;
366 if (m_IsMCFlag == 0) {
369 g_sih = &m_si[id][1];
377 if (m_CovarianceMatrix) unpackcovariance(m_c[
id]);
380 double p2_b, p2_a, p2_t, p2_a1, p2_chi2, p_extraPhotonEnergy, p_extraPhotonTime;
383 Fit2h(p2_b, p2_a, p2_t, p2_a1, p2_chi2);
387 if (p2_chi2 >= m_chi2Threshold27dof) {
391 Fit2hExtraPhoton(p2_b, p2_a, p2_t, p2_a1, p_extraPhotonEnergy, p_extraPhotonTime, p2_chi2);
395 if (p2_chi2 >= m_chi2Threshold25dof) {
399 Fit2h(p2_b, p2_a, p2_t, p2_a1, p2_chi2);
408 aECLDsp.setTwoComponentTotalAmp(p2_a + p2_a1);
410 aECLDsp.setTwoComponentHadronAmp(0.0);
411 aECLDsp.setTwoComponentDiodeAmp(p2_a1);
413 aECLDsp.setTwoComponentHadronAmp(p2_a1);
414 aECLDsp.setTwoComponentDiodeAmp(0.0);
416 aECLDsp.setTwoComponentChi2(p2_chi2);
417 aECLDsp.setTwoComponentTime(p2_t);
418 aECLDsp.setTwoComponentBaseline(p2_b);
419 aECLDsp.setTwoComponentFitType(fitType);
421 aECLDsp.setbackgroundPhotonEnergy(p_extraPhotonEnergy);
422 aECLDsp.setbackgroundPhotonTime(p_extraPhotonTime);
448 for (
int j = 0; j < 31; j++)
if (amax < fitA[j]) { amax = fitA[j]; jmax = j;}
449 double sumB0 = 0;
int jsum = 0;
450 for (
int j = 0; j < 31; j++)
if (j < jmax - 3 || jmax + 4 < j) { sumB0 += fitA[j]; ++jsum;}
451 double B0 = sumB0 / jsum;
453 if (amax < 0) amax = 10;
454 double T0 = dt * (4.5 - jmax);
458 m_Minit2h->mnparm(0,
"B", B0, 10, B0 / 1.5, B0 * 1.5, ierflg);
459 m_Minit2h->mnparm(1,
"Ag", A0, A0 / 20, 0, 2 * A0, ierflg);
460 m_Minit2h->mnparm(2,
"T", T0, 0.5, T0 - 2.5, T0 + 2.5, ierflg);
461 m_Minit2h->mnparm(3,
"Ah", 0., A0 / 20, -A0, 2 * A0, ierflg);
466 m_Minit2h->mnexcm(
"MIGRAD", arglist, 2, ierflg);
469 int nvpar, nparx, icstat;
470 m_Minit2h->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
474 m_Minit2h->GetParameter(0, B, ep);
475 m_Minit2h->GetParameter(1, Ag, ep);
476 m_Minit2h->GetParameter(2, T, ep);
477 m_Minit2h->GetParameter(3, Ah, ep);
486 double amax = 0;
int jmax = 6;
487 for (
int j = 0; j < 31; j++)
if (amax < fitA[j]) { amax = fitA[j]; jmax = j;}
489 double amax1 = 0;
int jmax1 = 6;
490 for (
int j = 0; j < 31; j++)
491 if (j < jmax - 3 || jmax + 4 < j) {
492 if (j == 0 && amax1 < fitA[j] && fitA[j + 1] < fitA[j]) { amax1 = fitA[j]; jmax1 = j;}
493 else if (j == 30 && amax1 < fitA[j] && fitA[j - 1] < fitA[j]) { amax1 = fitA[j]; jmax1 = j;}
494 else if (amax1 < fitA[j] && fitA[j + 1] < fitA[j] && fitA[j - 1] < fitA[j]) { amax1 = fitA[j]; jmax1 = j;}
497 double sumB0 = 0;
int jsum = 0;
498 for (
int j = 0; j < 31; j++)
if ((j < jmax - 3 || jmax + 4 < j) && (j < jmax1 - 3 || jmax1 + 4 < j)) { sumB0 += fitA[j]; ++jsum;}
499 double B0 = sumB0 / jsum;
500 amax -= B0; amax = std::max(10.0, amax);
501 amax1 -= B0; amax1 = std::max(10.0, amax1);
502 double T0 = dt * (4.5 - jmax);
503 double T01 = dt * (4.5 - jmax1);
505 double A0 = amax, A01 = amax1;
506 m_Minit2h2->mnparm(0,
"B", B0, 10, B0 / 1.5, B0 * 1.5, ierflg);
507 m_Minit2h2->mnparm(1,
"Ag", A0, A0 / 20, 0, 2 * A0, ierflg);
508 m_Minit2h2->mnparm(2,
"T", T0, 0.5, T0 - 2.5, T0 + 2.5, ierflg);
509 m_Minit2h2->mnparm(3,
"Ah", 0., A0 / 20, -A0, 2 * A0, ierflg);
510 m_Minit2h2->mnparm(4,
"A2", A01, A01 / 20, 0, 2 * A01, ierflg);
511 m_Minit2h2->mnparm(5,
"T2", T01 , 0.5 , T01 - 2.5, T01 + 2.5, ierflg);
516 m_Minit2h2->mnexcm(
"MIGRAD", arglist, 2, ierflg);
519 int nvpar, nparx, icstat;
520 m_Minit2h2->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
522 m_Minit2h2->GetParameter(0, B, ep);
523 m_Minit2h2->GetParameter(1, Ag, ep);
524 m_Minit2h2->GetParameter(2, T, ep);
525 m_Minit2h2->GetParameter(3, Ah, ep);
526 m_Minit2h2->GetParameter(4, A2, ep);
527 m_Minit2h2->GetParameter(5, T2, ep);
536 std::vector<double> p(s.begin() + 1, s.end());
537 p[1] = std::max(0.0029, p[1]);
538 p[4] = std::max(0.0029, p[4]);
543 dd_t t[(c_nt + c_ntail)*c_ndt];
544 dsp.
fillarray(
sizeof(t) /
sizeof(t[0]), t);
546 for (
int i = 0; i < c_nt * c_ndt; i++) m_F[i] = t[i];
547 for (
int i = 0; i < c_ntail; i++) m_F[c_nt * c_ndt + i] = t[c_nt * c_ndt + i * c_ndt];
548 const auto& Fm = *(std::end(m_F) - 2), &F0 = *(std::end(m_F) - 1);
549 m_r0 = F0.first / Fm.first;
550 m_r1 = F0.second / Fm.second;
559 const int iend0 = c_nt * c_ndt, iend1 = c_nt * c_ndt + c_ntail;
565 if (++A >= Aend)
return;
570 double x = t0 * c_idtn, ix = floor(x), w = x - ix;
572 double w2 = w * w, hw2 = 0.5 * w2, tw3 = ((1. / 6) * w) * w2;
573 auto I = [
this, &w, &hw2, &tw3](
int j,
double idt,
double dt) {
575 f0 = m_F[j].first, f1 = m_F[j + 1].first,
576 fp0 = m_F[j].second, fp1 = m_F[j + 1].second,
577 dfdt = (f1 - f0) * idt, fp = fp1 + fp0;
581 a[2] = -((fp + fp0) - 3 * dfdt);
582 a[3] = ((fp) - 2 * dfdt);
584 double b2 = 2 * a[2], b3 = 6 * a[3];
586 y.f0 = a[0] + dt * (a[1] * w + b2 * hw2 + b3 * tw3);
587 y.f1 = a[1] + b2 * w + b3 * hw2;
588 y.f2 = (b2 + b3 * w) * idt;
595 *A = I(i, c_idtn, c_dtn);
596 if (++A >= Aend)
return;
601 x = t0 * c_idt; ix = floor(x); w = x - ix;
603 i = (j - c_nt) + iend0;
604 w2 = w * w, hw2 = 0.5 * w2, tw3 = ((1. / 6) * w) * w2;
608 while (i < iend1 - 1) {
609 *A = I(i++, c_idt, c_dt);
610 if (++A >= Aend)
return;