Belle II Software  release-05-01-25
ECLWaveformFit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * This module is used to compute the two component (photon+hadron) *
6  * fit to ecl waveforms stored offline. Hadron component energy *
7  * measured from fit is used to perform pulse shape discrimination *
8  * for particle id. *
9  * *
10  * Author: The Belle II Collaboration *
11  * Contributors: Savino Longo (longos@uvic.ca) *
12  * Alexei Sibidanov (sibid@uvic.ca) *
13  * *
14  * *
15  * This software is provided "as is" without any warranty. *
16  **************************************************************************/
17 
18 // ECL
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>
26 
27 //FRAMEWORK
28 #include <framework/core/Environment.h>
29 #include <framework/database/DBObjPtr.h>
30 
31 //ROOT
32 #include <TMinuit.h>
33 #include <TMatrixD.h>
34 #include <TMatrixDSym.h>
35 #include <TDecompChol.h>
36 
37 #include <numeric>
38 
39 using namespace Belle2;
40 using namespace ECL;
41 
42 //-----------------------------------------------------------------
43 // Register the Modules
44 //-----------------------------------------------------------------
45 REG_MODULE(ECLWaveformFit)
46 //-----------------------------------------------------------------
47 // Implementation
48 //-----------------------------------------------------------------
49 
50 //anonymous namespace for data objects used by both ECLWaveformFitModule class and FCN2h funciton for MINUIT minimization.
51 namespace {
52 
53  //adc data array
54  std::vector<double> fitA(31);
55 
56  //g_si: photon template signal shape
57  //g_sih: hadron template signal shape
58  const SignalInterpolation2* g_si;
59  const SignalInterpolation2* g_sih;
60 
61  // covariance matrix and noise level
62  std::vector< std::vector<double> > currentCovMat;
63  double aNoise;
64 
65  //Function to minimize in photon template + hadron template fit. (chi2)
66  void FCN2h(int&, double* grad, double& f, double* p, int)
67  {
68  constexpr int N = 31;
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;
73 
74  //getting photon and hadron component shapes for set of fit parameters
75  val_der_t ADg[N], ADh[N];
76  g_si->getshape(T, ADg);
77  g_sih->getshape(T, ADh);
78 
79  //computing difference between current fit result and adc data array
80  for (int i = 0; i < N; ++i) df[i] = fitA[i] - (Ag * ADg[i].f0 + Ah * ADh[i].f0 + B);
81 
82  //computing chi2.
83  for (int i = 0; i < N; ++i) da[i] = std::inner_product(currentCovMat[i].begin(), currentCovMat[i].end(), df.begin(), 0.0);
84 
85  for (int i = 0; i < N; ++i) {
86  chi2 += da[i] * df[i];
87  gB -= da[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;
91  }
92 
93  f = chi2;
94  grad[0] = 2 * gB;
95  grad[1] = 2 * gAg;
96  grad[2] = 2 * gT;
97  grad[3] = 2 * gAh;
98  }
99 
100  //Function to minimize in photon template + hadron template + background photon fit. (chi2)
101  void FCN2h2(int&, double* grad, double& f, double* p, int)
102  {
103  const int N = 31;
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;
110 
111  //getting photon and hadron component shapes for set of fit parameters
112  val_der_t ADg[N], AD2[N], ADh[N];
113  g_si->getshape(T, ADg);
114  g_si->getshape(T2, AD2);//background photon
115  g_sih->getshape(T, ADh);
116 
117  //computing difference between current fit result and adc data array
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);
119 
120  //computing chi2.
121  for (int i = 0; i < N; ++i) da[i] = std::inner_product(currentCovMat[i].begin(), currentCovMat[i].end(), df.begin(), 0.0);
122 
123  for (int i = 0; i < N; ++i) {
124  chi2 += da[i] * df[i];
125  gB -= da[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);;
129 
130  gA2 -= da[i] * AD2[i].f0;
131  gT2 -= da[i] * AD2[i].f1 * A2;
132  }
133  f = chi2;
134  grad[0] = 2 * gB;
135  grad[1] = 2 * gAg;
136  grad[2] = 2 * gT;
137  grad[3] = 2 * gAh;
138  grad[4] = 2 * gA2;
139  grad[5] = 2 * gT2;
140  }
141 
142  // regularize autocovariance function by multipling it by the step
143  // function so elements above u0 become 0 and below are untouched.
144  void regularize(double* dst, const double* src, const int n, const double u0 = 13.0, const double u1 = 0.8)
145  {
146  for (int k = 0; k < n; k++) dst[k] = src[k] / (1 + exp((k - u0) / u1));
147  }
148 
149  // transform autocovariance function of 31 elements to the covariance matrix
150  bool makecovariance(CovariancePacked& M, const int nnoise, const double* acov)
151  {
152  const int ns = 31;
153  TMatrixDSym E(ns);
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];
157 
158  TDecompChol dc(E);
159  const bool status = dc.Invert(E);
160 
161  if (status) {
162  int count = 0;
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]);
167  }
168  return status;
169  }
170 
171  // to save space we keep only upper triangular part of the covariance matrix in float format
172  // here we inflate it to full square form in double format
173  void unpackcovariance(const CovariancePacked& matrixPacked)
174  {
175  const int ns = 31;
176  int count = 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++];
183  }
184  }
185  aNoise = matrixPacked.sigma;
186  }
187 
188 }
189 
190 // constructor
192 {
193  // Set module properties
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);
201 }
202 
203 // destructor
205 {
206 }
207 
208 //callback for loading templates from database
210 {
211 
212  m_TemplatesLoaded = true;
213 
214  if (m_IsMCFlag == 0) {
215  //load data templates
216  DBObjPtr<ECLDigitWaveformParameters> WavePars("ECLDigitWaveformParameters");
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];
223  }
224  new(&m_si[i][0]) SignalInterpolation2(Ptemp);
225  new(&m_si[i][1]) SignalInterpolation2(Htemp);
226  new(&m_si[i][2]) SignalInterpolation2(Dtemp);
227  }
228  } else {
229  //load mc template
230  DBObjPtr<ECLDigitWaveformParametersForMC> WaveParsMC("ECLDigitWaveformParametersForMC");
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];
236  }
237  new(&m_si[0][0]) SignalInterpolation2(Ptemp);
238  new(&m_si[0][1]) SignalInterpolation2(Htemp);
239  new(&m_si[0][2]) SignalInterpolation2(Dtemp);
240  }
241 }
242 
244 {
245 
246  m_IsMCFlag = Environment::Instance().isMC();
247  m_TemplatesLoaded = false;
248 
249  DBObjPtr<ECLCrystalCalib> Ael("ECLCrystalElectronics"), Aen("ECLCrystalEnergy");
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];
253 
254  //Load covariance matricies from database;
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());
262  double x0 = N;
263  reg = buf;
264  while (!makecovariance(m_c[id - 1], N, reg.data()))
265  regularize(buf.data(), reg.data(), N, x0 -= 1, 1);
266  }
267  } else {
268  //default covariance matrix is identity for all crystals
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;
276  }
277  }
278  }
279 
280 }
281 
283 {
284  // ECL dataobjects
285  m_eclDSPs.registerInDataStore();
286  m_eclDigits.registerInDataStore();
287 
288  //initializing fit minimizer
289  m_Minit2h = new TMinuit(4);
290  m_Minit2h->SetFCN(FCN2h);
291  double arglist[10];
292  int ierflg = 0;
293  arglist[0] = -1;
294  m_Minit2h->mnexcm("SET PRIntout", arglist, 1, ierflg);
295  m_Minit2h->mnexcm("SET NOWarnings", arglist, 0, ierflg);
296  arglist[0] = 1;
297  m_Minit2h->mnexcm("SET ERR", arglist, 1, ierflg);
298  arglist[0] = 0;
299  m_Minit2h->mnexcm("SET STRategy", arglist, 1, ierflg);
300  arglist[0] = 1;
301  m_Minit2h->mnexcm("SET GRAdient", arglist, 1, ierflg);
302  arglist[0] = 1e-6;
303  m_Minit2h->mnexcm("SET EPSmachine", arglist, 1, ierflg);
304 
305  //initializing fit minimizer photon+hadron + background photon
306  m_Minit2h2 = new TMinuit(6);
307  m_Minit2h2->SetFCN(FCN2h2);
308  arglist[0] = -1;
309  m_Minit2h2->mnexcm("SET PRIntout", arglist, 1, ierflg);
310  m_Minit2h2->mnexcm("SET NOWarnings", arglist, 0, ierflg);
311  arglist[0] = 1;
312  m_Minit2h2->mnexcm("SET ERR", arglist, 1, ierflg);
313  arglist[0] = 0;
314  m_Minit2h2->mnexcm("SET STRategy", arglist, 1, ierflg);
315  arglist[0] = 1;
316  m_Minit2h2->mnexcm("SET GRAdient", arglist, 1, ierflg);
317  arglist[0] = 1e-6;
318  m_Minit2h2->mnexcm("SET EPSmachine", arglist, 1, ierflg);
319 
320  //flag for callback to load templates each run
321  m_TemplatesLoaded = false;
322 }
323 
325 {
327 
328  if (!m_TemplatesLoaded) {
329  //load templates once per run in first event that has saved waveforms.
330  if (m_eclDSPs.getEntries() > 0) loadTemplateParameterArray();
331  }
332 
333  for (auto& aECLDsp : m_eclDSPs) {
334 
335  aECLDsp.setTwoComponentTotalAmp(-1);
336  aECLDsp.setTwoComponentHadronAmp(-1);
337  aECLDsp.setTwoComponentDiodeAmp(-1);
338  aECLDsp.setTwoComponentChi2(-1);
339  aECLDsp.setTwoComponentSavedChi2(ECLDsp::photonHadron, -1);
340  aECLDsp.setTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton, -1);
341  aECLDsp.setTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing, -1);
342  aECLDsp.setTwoComponentTime(-1);
343  aECLDsp.setTwoComponentBaseline(-1);
344  aECLDsp.setTwoComponentFitType(ECLDsp::poorChi2);
345 
346  const int id = aECLDsp.getCellId() - 1;
347 
348  //Filling array with ADC values.
349  for (int j = 0; j < ec.m_nsmp; j++) fitA[j] = aECLDsp.getDspA()[j];
350 
351  //setting relation of eclDSP to aECLDigit
352  const ECLDigit* d = NULL;
353  for (const auto& aECLDigit : m_eclDigits) {
354  if (aECLDigit.getCellId() - 1 == id) {
355  d = &aECLDigit;
356  aECLDsp.addRelationTo(&aECLDigit);
357  break;
358  }
359  }
360  if (d == NULL) continue;
361 
362  //Skipping low amplitude waveforms
363  if (d->getAmp() * m_ADCtoEnergy[id] < m_EnergyThreshold) continue;
364 
365  //loading template for waveform
366  if (m_IsMCFlag == 0) {
367  //data cell id dependent
368  g_si = &m_si[id][0];
369  g_sih = &m_si[id][1];
370  } else {
371  // mc uses same waveform
372  g_si = &m_si[0][0];
373  g_sih = &m_si[0][1];
374  }
375 
376  //get covariance matrix for cell id
377  if (m_CovarianceMatrix) unpackcovariance(m_c[id]);
378 
379  //Calling optimized fit photon template + hadron template (fit type = 0)
380  double p2_b, p2_a, p2_t, p2_a1, p2_chi2, p_extraPhotonEnergy, p_extraPhotonTime;
382  p2_chi2 = -1;
383  Fit2h(p2_b, p2_a, p2_t, p2_a1, p2_chi2);
384  aECLDsp.setTwoComponentSavedChi2(ECLDsp::photonHadron, p2_chi2);
385 
386  //if hadron fit failed try hadron + background photon (fit type = 1)
387  if (p2_chi2 >= m_chi2Threshold27dof) {
388 
390  p2_chi2 = -1;
391  Fit2hExtraPhoton(p2_b, p2_a, p2_t, p2_a1, p_extraPhotonEnergy, p_extraPhotonTime, p2_chi2);
392  aECLDsp.setTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton, p2_chi2);
393 
394  //hadron + background photon fit failed try diode fit (fit type = 2)
395  if (p2_chi2 >= m_chi2Threshold25dof) {
396  g_sih = &m_si[0][2];//set second component to diode
397  fitType = ECLDsp::photonDiodeCrossing;
398  p2_chi2 = -1;
399  Fit2h(p2_b, p2_a, p2_t, p2_a1, p2_chi2);
400  aECLDsp.setTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing, p2_chi2);
401 
402  if (p2_chi2 >= m_chi2Threshold27dof) fitType = ECLDsp::poorChi2; //indicates all fits tried had bad chi2
403  }
404 
405  }
406 
407  //storing fit results
408  aECLDsp.setTwoComponentTotalAmp(p2_a + p2_a1);
409  if (fitType == ECLDsp::photonDiodeCrossing) {
410  aECLDsp.setTwoComponentHadronAmp(0.0);
411  aECLDsp.setTwoComponentDiodeAmp(p2_a1);
412  } else {
413  aECLDsp.setTwoComponentHadronAmp(p2_a1);
414  aECLDsp.setTwoComponentDiodeAmp(0.0);
415  }
416  aECLDsp.setTwoComponentChi2(p2_chi2);
417  aECLDsp.setTwoComponentTime(p2_t);
418  aECLDsp.setTwoComponentBaseline(p2_b);
419  aECLDsp.setTwoComponentFitType(fitType);
420  if (fitType == ECLDsp::photonHadronBackgroundPhoton) {
421  aECLDsp.setbackgroundPhotonEnergy(p_extraPhotonEnergy);
422  aECLDsp.setbackgroundPhotonTime(p_extraPhotonTime);
423  }
424  }
425 }
426 
427 // end run
429 {
430 }
431 
432 // terminate
434 {
435 }
436 
437 //Two component fit. Fits photon+hadron component hypothesis
438 void ECLWaveformFitModule::Fit2h(double& B, double& Ag, double& T, double& Ah, double& amin)
439 {
440  //minuit parameters
441  double arglist[10];
442  int ierflg = 0;
443 
444  // setting inital fit parameters
445  double dt = 0.5;
446  double amax = 0;
447  int jmax = 6;
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;
452  amax -= B0;
453  if (amax < 0) amax = 10;
454  double T0 = dt * (4.5 - jmax);
455  double A0 = amax;
456 
457  //initalize minimizer
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);
462 
463  //perform fit
464  arglist[0] = 50000;
465  arglist[1] = 1.;
466  m_Minit2h->mnexcm("MIGRAD", arglist, 2, ierflg);
467 
468  double edm, errdef;
469  int nvpar, nparx, icstat;
470  m_Minit2h->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
471 
472  //get fit results
473  double ep;
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);
478 }
479 
480 //Two component fit background ith extra photon. Fits photon+hadron component + extra photon hypothesis
481 void ECLWaveformFitModule::Fit2hExtraPhoton(double& B, double& Ag, double& T, double& Ah, double& A2, double& T2, double& amin)
482 {
483  double arglist[10];
484  int ierflg = 0;
485  double dt = 0.5;
486  double amax = 0; int jmax = 6;
487  for (int j = 0; j < 31; j++) if (amax < fitA[j]) { amax = fitA[j]; jmax = j;}
488 
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;}
495  }
496 
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);
504 
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);
512 
513  // Now ready for minimization step
514  arglist[0] = 50000;
515  arglist[1] = 1.;
516  m_Minit2h2->mnexcm("MIGRAD", arglist, 2, ierflg);
517 
518  double edm, errdef;
519  int nvpar, nparx, icstat;
520  m_Minit2h2->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
521  double ep;
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);
528 }
529 
530 
531 
532 //Signal interpolation code used for fast evaluation of shaperDSP templates
533 SignalInterpolation2::SignalInterpolation2(const std::vector<double>& s)
534 {
535  double T0 = -0.2;
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]);
539 
540  ShaperDSP_t dsp(p, s[0]);
541  dsp.settimestride(c_dtn);
542  dsp.settimeseed(T0);
543  dd_t t[(c_nt + c_ntail)*c_ndt];
544  dsp.fillarray(sizeof(t) / sizeof(t[0]), t);
545 
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;
551 }
552 
557 void SignalInterpolation2::getshape(double t0, val_der_t* A) const
558 {
559  const int iend0 = c_nt * c_ndt, iend1 = c_nt * c_ndt + c_ntail;
560  const val_der_t* Aend = A + 31;
561 
562  //if before pulse start time (negative times) return 0
563  while (t0 < 0) {
564  *A = {0, 0, 0};
565  if (++A >= Aend) return;
566  t0 += c_dt;
567  }
568 
569  //function below evaluates the template value and the first and second derivative values for the point.
570  double x = t0 * c_idtn, ix = floor(x), w = x - ix;
571  int i = 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) {
574  double a[4],
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;
578 
579  a[0] = f0;
580  a[1] = fp0;
581  a[2] = -((fp + fp0) - 3 * dfdt);
582  a[3] = ((fp) - 2 * dfdt);
583 
584  double b2 = 2 * a[2], b3 = 6 * a[3];
585  val_der_t y;
586  y.f0 = a[0] + dt * (a[1] * w + b2 * hw2 + b3 * tw3); //function value
587  y.f1 = a[1] + b2 * w + b3 * hw2; // first derivative of function value
588  y.f2 = (b2 + b3 * w) * idt; //second derivative of function value
589  return y;
590  };
591 
592  //signal interpolation for short time steps used for points at the beginning of the pulse where the pulse is quickly changing (eg rise)
593  //iend0 indicates first region of pulse
594  while (i < iend0) {
595  *A = I(i, c_idtn, c_dtn);
596  if (++A >= Aend) return;
597  i += c_ndt;
598  t0 += c_dt;
599  }
600 
601  x = t0 * c_idt; ix = floor(x); w = x - ix;
602  int j = ix;
603  i = (j - c_nt) + iend0;
604  w2 = w * w, hw2 = 0.5 * w2, tw3 = ((1. / 6) * w) * w2;
605 
606  //signal interpolation for long time steps used for points in the tail region of the pulse
607  //iend1 indicates end of pulse
608  while (i < iend1 - 1) {
609  *A = I(i++, c_idt, c_dt);
610  if (++A >= Aend) return;
611  }
612 
613  while (A < Aend) {
614  const val_der_t& p = *(A - 1);
615  val_der_t& y = *A++;
616  y.f0 = p.f0 * m_r0; // function value
617  y.f1 = p.f1 * m_r1; // first derivative of function value
618  y.f2 = 0; // second derivative of function value
619  }
620 }
Belle2::ECL::ShaperDSP_t::settimeseed
void settimeseed(double)
set initial time
Definition: shaperdsp.cc:363
Belle2::ECLDsp::photonHadron
@ photonHadron
photon + hadron template fit
Definition: ECLDsp.h:44
Belle2::ECLWaveformFitModule::beginRun
virtual void beginRun() override
begin run.
Definition: ECLWaveformFit.cc:243
Belle2::ECL::ShaperDSP_t
Class include function that calculate electronic response from energy deposit
Definition: shaperdsp.h:36
Belle2::ECL::EclConfiguration::m_nsmp
static constexpr int m_nsmp
number of ADC measurements for signal fitting
Definition: EclConfiguration.h:51
Belle2::ECLDsp::TwoComponentFitType
TwoComponentFitType
Offline two component fit type.
Definition: ECLDsp.h:42
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ECLWaveformFitModule::loadTemplateParameterArray
void loadTemplateParameterArray()
loads waveform templates from database.
Definition: ECLWaveformFit.cc:209
Belle2::val_der_t
Struct to return signal function information f0 is the function value f1 is the first derivative f2 i...
Definition: ECLWaveformFit.h:51
Belle2::ECL::ShaperDSP_t::settimestride
void settimestride(double)
set grid step for function calculation
Definition: shaperdsp.cc:353
Belle2::ECL::ShaperDSP_t::fillarray
void fillarray(int, double *) const
fill array for amplitude and time calculation
Definition: shaperdsp.cc:373
Belle2::Environment::isMC
bool isMC() const
Do we have generated, not real data?
Definition: Environment.cc:57
Belle2::val_der_t::f1
double f1
see struct description
Definition: ECLWaveformFit.h:55
Belle2::ECLWaveformFitModule::ECLWaveformFitModule
ECLWaveformFitModule()
Constructor.
Definition: ECLWaveformFit.cc:191
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::ECL::EclConfiguration
singleton class to hold the ECL configuration
Definition: EclConfiguration.h:31
Belle2::CovariancePacked
Struct to keep upper triangle of the covariance matrix.
Definition: ECLWaveformFit.h:35
Belle2::ECLWaveformFitModule::initialize
virtual void initialize() override
Initialize variables.
Definition: ECLWaveformFit.cc:282
Belle2::val_der_t::f0
double f0
see struct description
Definition: ECLWaveformFit.h:53
Belle2::ECLDsp::poorChi2
@ poorChi2
All offline fit attempts were greater than chi2 threshold.
Definition: ECLDsp.h:43
Belle2::ECLWaveformFitModule::Fit2hExtraPhoton
void Fit2hExtraPhoton(double &b, double &a0, double &t0, double &a1, double &A2, double &T2, double &chi2)
Optimized fit using hadron component model plus out of time background photon.
Definition: ECLWaveformFit.cc:481
Belle2::SignalInterpolation2::getshape
void getshape(double, val_der_t *) const
returns signal shape(+derivatives) in 31 equidistant time points starting from T0
Definition: ECLWaveformFit.cc:557
Belle2::ECL::EclConfiguration::get
static EclConfiguration & get()
return this instance
Definition: EclConfiguration.h:34
Belle2::ECLWaveformFitModule::~ECLWaveformFitModule
~ECLWaveformFitModule()
Destructor.
Definition: ECLWaveformFit.cc:204
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECLDsp::photonDiodeCrossing
@ photonDiodeCrossing
photon + diode template fit
Definition: ECLDsp.h:46
Belle2::ECLDigit
Class to store ECL digitized hits (output of ECLDigi) relation to ECLHit filled in ecl/modules/eclDig...
Definition: ECLDigit.h:34
Belle2::ECLWaveformFitModule::Fit2h
void Fit2h(double &b, double &a0, double &t0, double &a1, double &chi2)
Optimized fit using hadron component model.
Definition: ECLWaveformFit.cc:438
Belle2::SignalInterpolation2
Interpolate signal shape using function values and the first derivative.
Definition: ECLWaveformFit.h:62
Belle2::ECLWaveformFitModule::event
virtual void event() override
event per event.
Definition: ECLWaveformFit.cc:324
Belle2::CovariancePacked::sigma
float sigma
sigma noise
Definition: ECLWaveformFit.h:39
Belle2::ECLWaveformFitModule::terminate
virtual void terminate() override
terminate.
Definition: ECLWaveformFit.cc:433
Belle2::ECLDsp::photonHadronBackgroundPhoton
@ photonHadronBackgroundPhoton
photon + hadron template + pile-up photon fit
Definition: ECLDsp.h:45
Belle2::SignalInterpolation2::SignalInterpolation2
SignalInterpolation2()
Default constructor.
Definition: ECLWaveformFit.h:95
Belle2::Environment::Instance
static Environment & Instance()
Static method to get a reference to the Environment instance.
Definition: Environment.cc:31
Belle2::ECLWaveformFitModule::endRun
virtual void endRun() override
end run.
Definition: ECLWaveformFit.cc:428