Belle II Software  release-05-01-25
ECLDspFitterPure.cc
1 #include <ecl/digitization/ECLDspFitterPure.h>
2 #include <TVectorD.h>
3 #include <TMatrixD.h>
4 #include <TDecompLU.h>
5 #include <numeric>
6 #include <algorithm>
7 #include <iostream>
8 
9 using namespace Belle2;
10 using namespace Belle2::ECL;
11 using namespace std;
12 
13 double Belle2::ECL::func(int i, int ifine,
14  EclConfigurationPure::signalsamplepure_t& signal)
15 {
16  if (i == 0) return 0;
17  if (ifine < 0) {
18  i--;
19  if (i == 0) return 0;
20  ifine = EclConfigurationPure::m_ndtPure + ifine;
21  }
22  return signal.m_ft[(i - 1) * EclConfigurationPure::m_ndtPure + ifine ];
23 }
24 
25 double Belle2::ECL::func1(int i, int ifine,
26  EclConfigurationPure::signalsamplepure_t& signal)
27 {
28  if (i == 0) return 0;
29  if (ifine < 0) {
30  i--;
31  if (i == 0) return 0;
32  ifine = EclConfigurationPure::m_ndtPure + ifine;
33  }
34  return signal.m_ft1[(i - 1) * EclConfigurationPure::m_ndtPure + ifine ];
35 }
36 
37 void Belle2::ECL::initParams(EclConfigurationPure::fitparamspure_t& params,
38  EclConfigurationPure::signalsamplepure_t& signal)
39 {
40 
41  for (int i = 0; i < 16; i++) {
42  for (int k = -EclConfigurationPure::m_ndtPure; k < EclConfigurationPure::m_ndtPure; k++) {
43  // cout << i << ":" << k << " " << func( i, k, signal) << " " << func1(i, k, signal) << endl;
44  if (k < 0) params.f[i][ EclConfigurationPure::m_ndtPure - 1 - k] = func(i, k, signal);
45  else params.f[i][k] = func(i, k, signal);
46  }
47  }
48 
49  double kdt = 0.001 / (EclConfigurationPure::m_tickPure / EclConfigurationPure::m_ndtPure) ; // df / dtau = df / dt * dt / dtau
50 
51  params.c002 = 0;
52  for (int h = 0; h < 16; h++) {
53  params.c001[h] = 0;
54  for (int j = 0; j < 16; j++) {
55  params.c001[h] -= params.invC[h][j];
56  params.c002 += params.invC[h][j];
57  }
58  params.c001[h] *= 2;
59  }
60 
61  for (int k = 0; k < EclConfigurationPure::m_ndtPure; k++) {
62  params.c110[k] = 0;
63  params.c200[k] = 0;
64  params.c020[k] = 0;
65  params.c101[k] = 0;
66  params.c011[k] = 0;
67  for (int i = 0; i < 16; i++) {
68  for (int j = 0; j < 16; j++) {
69  params.c110[k] += kdt * func1(i, k, signal) * params.invC[i][j] * func(j, k, signal);
70  params.c200[k] += func(i, k, signal) * params.invC[i][j] * func(j, k, signal);
71  params.c020[k] += kdt * func1(i, k, signal) * params.invC[i][j] * func1(j, k, signal) * kdt;
72  params.c101[k] += func(i, k, signal) * params.invC[i][j];
73  params.c011[k] += kdt * func1(i, k, signal) * params.invC[i][j];
74  }
75  }
76  params.c110[k] *= 2;
77  params.c101[k] *= 2;
78  params.c011[k] *= 2;
79  }
80 
81  int off = EclConfigurationPure::m_ndtPure - 1;
82  for (int k = 1; k < EclConfigurationPure::m_ndtPure; k++) {
83  params.c110[k + off] = 0;
84  params.c200[k + off] = 0;
85  params.c020[k + off] = 0;
86  params.c101[k + off] = 0;
87  params.c011[k + off] = 0;
88  for (int i = 0; i < 16; i++) {
89  for (int j = 0; j < 16; j++) {
90  params.c110[k + off] += kdt * func1(i, -k, signal) * params.invC[i][j] * func(j, -k, signal);
91  params.c200[k + off] += func(i, -k, signal) * params.invC[i][j] * func(j, -k, signal);
92  params.c020[k + off] += kdt * func1(i, -k, signal) * params.invC[i][j] * func1(j, -k, signal) * kdt;
93  params.c101[k + off] += func(i, -k, signal) * params.invC[i][j];
94  params.c011[k + off] += kdt * func1(i, -k, signal) * params.invC[i][j];
95  }
96  }
97  params.c110[k + off] *= 2;
98  params.c101[k + off] *= 2;
99  params.c011[k + off] *= 2;
100  }
101 
102 
103  for (int h = 0; h < 16; h++) {
104 
105  for (int k = 0; k < EclConfigurationPure::m_ndtPure; k++) {
106  params.c100[h][k] = 0;
107  params.c010[h][k] = 0;
108  for (int i = 0; i < 16; i++) {
109  params.c100[h][k] -= func(i, k, signal) * params.invC[i][h];
110  params.c010[h][k] -= kdt * func1(i, k, signal) * params.invC[i][h];
111  }
112  params.c100[h][k] *= 2;
113  params.c010[h][k] *= 2;
114  }
115  for (int k = 1; k < EclConfigurationPure::m_ndtPure; k++) {
116  params.c100[h][k + off] = 0;
117  params.c010[h][k + off] = 0;
118  for (int i = 0; i < 16; i++) {
119  params.c100[h][k + off] -= func(i, -k, signal) * params.invC[i][h];
120  params.c010[h][k + off] -= kdt * func1(i, -k, signal) * params.invC[i][h];
121  }
122  params.c100[h][k + off] *= 2;
123  params.c010[h][k + off] *= 2;
124  }
125 
126  }
127 }
128 
129 void Belle2::ECL::DSPFitterPure(const EclConfigurationPure::fitparamspure_t& f , const int* FitA, const int ttrig, int& amp,
130  double& time, double& chi2, int& iter)
131 {
132  int baseline = accumulate(FitA, FitA + 16, 0) / 16;
133  const int* imax = max_element(FitA, FitA + EclConfigurationPure::m_nsmp);
134 
135  amp = *imax - baseline;
136  //time = imax - FitA;
137 
138 
139  int dt = ttrig;
140  double y[16];
141  y[0] = baseline;
142  for (int i = 16; i < EclConfigurationPure::m_nsmp; i++) y[i - 15] = FitA[i];
143 
144  //cout << "ttrig = " << dt << endl;
145  // for (int i=0; i<16; i++ ) cout << "y_" << i << " = " << y[i] << endl;
146 
147  double k100 { 0 }, k010{ 0 }, k001{ 0 };
148  int offset = EclConfigurationPure::m_ndtPure;
149  if (dt < 0) dt = offset - 1 - dt;
150  for (int i = 0; i < 16; i++) {
151  k100 += f.c100[i][dt] * y[i];
152  k010 += f.c010[i][dt] * y[i];
153  k001 += f.c001[i] * y[i];
154  }
155 
156  double a[] = { 2 * f.c200[dt], f.c110[dt], f.c101[dt] ,
157  f.c110[dt], 2 * f.c020[dt], f.c011[dt],
158  f.c101[dt], f.c011[dt], 2 * f.c002
159  };
160 
161  TMatrixD M1(3, 3);
162  M1.SetMatrixArray(a);
163 
164  TVectorD C(3);
165  C[0] = - k100;
166  C[1] = - k010;
167  C[2] = - k001;
168  /*
169  cout << "(" << -k100 << "==" << C[0] << ","
170  << "(" << -k010 << "==" << C[1] << ","
171  << "(" << -k001 << "==" << C[2] << ","
172  << endl;
173 
174  */
175  // C.Print();
176  TDecompLU solver(M1);
177  bool ok;
178  C = solver.Solve(C, ok);
179  /*
180  if ( !ok ) {
181  cout << "Fit failed" << endl;
182  cout << "dt = " << dt << endl;
183  M1.Print();
184  cout << "solution -->" << endl;
185  C.Print();
186  }
187  */
188  amp = C[0];
189 
190  /*
191  cout << "t0 = " << ttrig << endl;
192  cout << "P = " << C[2] << endl;
193  cout << "A = " << C[0] << endl;
194  cout << "B = A * t = " << C[1] << endl;
195  cout << "dt = B / A = " << C[1] / C[0] << endl;
196  */
197 
198  double A = C[0];
199  double B = C[1];
200  double P = C[2];
201 
202  /*
203  for (int i=1; i<16; i++)
204  cout << "A= " << C[0] << " f = " << f.f[i][dt] << " A*f = " << C[0]*f.f[i][dt] << " y = " << y[i] << endl;
205  */
206 
207  chi2 = f.c200[dt] * A * A + f.c020[dt] * B * B + f.c002 * P * P +
208  f.c110[dt] * A * B + f.c101[dt] * A * P + f.c011[dt] * B * P +
209  k100 * A + k010 * B + k001 * P ;
210  for (int i = 0; i < 16; i++)
211  for (int j = 0; j < 16; j++)
212  chi2 += y[i] * f.invC[i][j] * y[j];
213 
214  // cout << "chi2 = " << chi2 << endl;
215  double deltaT = C[1] / C[0] ;
216  int intTime = ttrig + round(deltaT) ;
217  // time = -ttrig - C[1] / C[0] ;
218  time = -(ttrig + deltaT);
219  // cout << "time : " << time << endl;
220  bool intime = abs(intTime) < EclConfigurationPure::m_ndtPure - 1;
221  if (iter < 3 && intime) {
222  iter++;
223  DSPFitterPure(f, FitA, intTime, amp, time, chi2, iter);
224  }
225 
226  if (intime && abs(deltaT) > 0.5 && iter++ < 10)
227  DSPFitterPure(f, FitA, intTime, amp, time, chi2, iter);
228 }
Belle2::ECL::EclConfigurationPure::fitparamspure_t
a struct for the fit parameters for the pure CsI calorimeter
Definition: EclConfigurationPure.h:57
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::EclConfigurationPure::m_ndtPure
static constexpr int m_ndtPure
number of points per ADC tick where signal fit procedure parameters are evaluated
Definition: EclConfigurationPure.h:36
Belle2::ECL::EclConfigurationPure::m_nsmp
static constexpr int m_nsmp
number of ADC measurements for signal fitting
Definition: EclConfigurationPure.h:29
Belle2::ECL::EclConfigurationPure::signalsamplepure_t
a struct for a signal sample for the pure CsI calorimeter
Definition: EclConfigurationPure.h:40