Belle II Software  release-08-00-04
YScanner.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <top/reconstruction_cpp/YScanner.h>
10 #include <top/reconstruction_cpp/func.h>
11 #include <top/geometry/TOPGeometryPar.h>
12 #include <framework/logging/Logger.h>
13 #include <cmath>
14 
15 namespace Belle2 {
20  namespace TOP {
21 
23 
25  const InverseRaytracer::Solution& sol_dx,
26  const InverseRaytracer::Solution& sol_de,
27  const InverseRaytracer::Solution& sol_dL)
28  {
29  if (sol_dx.step == 0) B2ERROR("TOP::YScanner::Derivatives: step (dx) is zero");
30  if (sol_de.step == 0) B2ERROR("TOP::YScanner::Derivatives: step (de) is zero");
31  if (sol_dL.step == 0) B2ERROR("TOP::YScanner::Derivatives: step (dL) is zero");
32 
33  dLen_dx = dLen_d(sol, sol_dx);
34  dLen_de = dLen_d(sol, sol_de);
35  dLen_dL = dLen_d(sol, sol_dL);
36 
37  dyB_dx = dyB_d(sol, sol_dx);
38  dyB_de = dyB_d(sol, sol_de);
39  dyB_dL = dyB_d(sol, sol_dL);
40 
41  dFic_dx = dFic_d(sol, sol_dx);
42  dFic_de = dFic_d(sol, sol_de);
43  dFic_dL = dFic_d(sol, sol_dL);
44  }
45 
46 
47  YScanner::YScanner(int moduleID, unsigned N): RaytracerBase(moduleID, c_Unified, c_SemiLinear),
49  m_pixelMasks(PixelMasks(moduleID)),
51  {
52  if (N < 2) {
53  B2FATAL("TOP::YScanner: N must be > 1");
54  return;
55  }
56 
57  // set the table of nominal photon detection efficiencies (incl. wavelength filter)
58 
59  const auto* topgp = TOPGeometryPar::Instance();
60  const auto* geo = topgp->getGeometry();
61  auto qe = geo->getNominalQE(); // get a copy
62  qe.applyFilterTransmission(geo->getWavelengthFilter());
63 
64  double minE = TOPGeometryPar::c_hc / qe.getMaxLambda();
65  double maxE = TOPGeometryPar::c_hc / qe.getMinLambda();
66  if (minE >= maxE) {
67  B2FATAL("TOP::YScanner: quantum efficiency found zero for all wavelengths");
68  return;
69  }
70  m_efficiency.set(minE, (maxE - minE) / (N - 1));
71 
72  const auto& tdc = geo->getNominalTDC();
73  for (unsigned i = 0; i < N; i++) {
74  double e = m_efficiency.getX(i);
75  double lambda = TOPGeometryPar::c_hc / e;
76  double effi = qe.getEfficiency(lambda) * tdc.getEfficiency();
77  m_efficiency.entries.push_back(TableEntry(effi, e, e * e));
78  }
79 
80  // set cosine of total reflection angle using photon mean energy for beta = 1
81 
82  double s = 0;
83  double se = 0;
84  double see = 0;
85  for (const auto& entry : m_efficiency.entries) {
86  double e = entry.x;
87  double p = entry.y * (1 - 1 / pow(topgp->getPhaseIndex(e), 2));
88  s += p;
89  se += p * e;
90  see += p * e * e;
91  }
92  if (s == 0) return;
93  m_meanE0 = se / s;
94  m_rmsE0 = sqrt(see / s - m_meanE0 * m_meanE0);
95  m_cosTotal = sqrt(1 - 1 / pow(topgp->getPhaseIndex(m_meanE0), 2));
96  }
97 
98 
99  void YScanner::clear() const
100  {
101  m_momentum = 0;
102  m_beta = 0;
103  m_length = 0;
104  m_numPhotons = 0;
105  m_meanE = 0;
106  m_rmsE = 0;
107  m_sigmaScat = 0;
108  m_sigmaAlpha = 0;
111  m_quasyEnergyDistribution = nullptr;
112  m_aboveThreshold = false;
113  m_results.clear();
114  m_scanDone = false;
115  }
116 
117 
118  void YScanner::prepare(double momentum, double beta, double length) const
119  {
120  clear();
121 
122  m_momentum = momentum;
123  m_beta = beta;
124  m_length = length;
125 
126  // set photon energy distribution, and the mean and r.m.s of photon energy
127 
128  auto area = setEnergyDistribution(beta);
129  if (area == 0) return;
130 
131  // set number of Cerenkov photons per azimuthal angle
132 
133  m_numPhotons = 370 * area / (2 * M_PI);
134 
135  // set multiple scattering and surface roughness sigmas in photon energy units
136 
137  const double radLength = 12.3; // quartz radiation length [cm]
138  double thetaScat = 13.6e-3 / beta / momentum * sqrt(length / 2 / radLength); // r.m.s of multiple scattering angle
139 
140  const auto* topgp = TOPGeometryPar::Instance();
141  double n = topgp->getPhaseIndex(m_meanE);
142  if (beta * n < 1) {
143  B2ERROR("TOP::YScanner::prepare: beta * n < 1 ==> must be a bug!");
144  return;
145  }
146  double dndE = topgp->getPhaseIndexDerivative(m_meanE);
147  double dEdTheta = n * sqrt(pow(beta * n, 2) - 1) / dndE;
148  m_sigmaScat = std::abs(thetaScat * dEdTheta); // r.m.s of multiple scattering angle converted to photon energy
149  m_sigmaAlpha = std::abs(m_bars.back().sigmaAlpha * dEdTheta); // surface roughness converted to photon energy
150 
151  // set photon energy distribution convoluted with multiple scattering
152 
154 
155  m_aboveThreshold = true;
156  }
157 
158 
159  double YScanner::setEnergyDistribution(double beta) const
160  {
161  const auto* topgp = TOPGeometryPar::Instance();
162 
164 
165  double s = 0;
166  double se = 0;
167  double see = 0;
168  for (const auto& entry : m_efficiency.entries) {
169  double e = entry.x;
170  double p = std::max(entry.y * (1 - 1 / pow(beta * topgp->getPhaseIndex(e), 2)), 0.0);
171  double ee = entry.xsq;
172  m_energyDistribution.entries.push_back(TableEntry(p, e, ee));
173  s += p;
174  se += p * e;
175  see += p * ee;
176  }
177  if (s == 0) return 0;
178 
179  for (auto& entry : m_energyDistribution.entries) entry.y /= s;
180 
181  m_meanE = se / s;
182  m_rmsE = sqrt(see / s - m_meanE * m_meanE);
183 
184  return s * m_energyDistribution.step;
185  }
186 
187 
188  void YScanner::setQuasyEnergyDistribution(double sigma) const
189  {
190  if (m_quasyEnergyDistributions.size() > 1000) {
192  B2ERROR("TOP::YScanner:setQuasyEnergyDistribution: unexpectedly large size of the std::map found, map cleared");
193  }
194 
195  double step = m_energyDistribution.step;
196  int ng = lround(3 * sigma / step);
197  auto& quasyEnergyDistribution = m_quasyEnergyDistributions[ng];
198 
199  if (quasyEnergyDistribution.entries.empty()) {
200  std::vector<double> gaus;
201  for (int i = 0; i <= ng; i++) {
202  double x = step * i / sigma;
203  gaus.push_back(exp(-0.5 * x * x));
204  }
205 
206  quasyEnergyDistribution.set(m_energyDistribution.getX(-ng), step);
207  int N = m_energyDistribution.entries.size();
208  double sum = 0;
209  for (int k = -ng; k < N + ng; k++) {
210  double s = 0;
211  double se = 0;
212  double see = 0;
213  for (int i = -ng; i <= ng; i++) {
214  double p = gaus[std::abs(i)] * m_energyDistribution.getY(k - i);
215  double e = m_energyDistribution.getX(k - i);
216  s += p;
217  se += p * e;
218  see += p * e * e;
219  }
220  if (s > 0) {
221  se /= s;
222  see /= s;
223  }
224  quasyEnergyDistribution.entries.push_back(TableEntry(s, se, see));
225  sum += s;
226  }
227  for (auto& entry : quasyEnergyDistribution.entries) entry.y /= sum;
228  }
229 
230  m_quasyEnergyDistribution = &quasyEnergyDistribution;
231  }
232 
233 
234  void YScanner::expand(unsigned col, double yB, double dydz, const Derivatives& D, int Ny, bool doScan) const
235  {
236  m_results.clear();
237 
238  if (D.dyB_de == 0) return;
239 
240  double sigma = sqrt(pow(m_sigmaScat, 2) + pow(m_sigmaAlpha, 2) * std::abs(Ny));
242 
243  double minE = m_quasyEnergyDistribution->getXmin();
244  double maxE = m_quasyEnergyDistribution->getXmax();
245  double pixDx = m_pixelPositions.get(col + 1).Dx;
246  double dely = (std::abs(D.dyB_dL) * m_length + std::abs(D.dyB_dx) * pixDx) / 2;
247  double y1 = yB - dely;
248  double y2 = yB + dely;
249  if (D.dyB_de > 0) {
250  y1 += D.dyB_de * (minE - m_meanE);
251  y2 += D.dyB_de * (maxE - m_meanE);
252  } else {
253  y1 += D.dyB_de * (maxE - m_meanE);
254  y2 += D.dyB_de * (minE - m_meanE);
255  }
256  double B = m_bars.front().B;
257  int j1 = lround(y1 / B);
258  int j2 = lround(y2 / B) + 1;
259 
260  if (doScan and j2 - j1 <= s_maxReflections) {
261  scan(col, yB, dydz, D, j1, j2);
262  m_scanDone = true;
263  } else {
264  merge(col, dydz, j1, j2);
265  m_scanDone = false;
266  }
267  }
268 
269 
270  void YScanner::scan(unsigned col, double yB, double dydz, const Derivatives& D, int j1, int j2) const
271  {
272 
273  std::map<int, EnergyMask*> masks;
274  for (unsigned row = 0; row < m_pixelPositions.getNumPixelRows(); row++) {
275 
276  int pixelID = m_pixelPositions.pixelID(row, col);
277  if (not m_pixelMasks.isActive(pixelID)) continue;
278 
279  const auto& pixel = m_pixelPositions.get(pixelID);
280  std::vector<PixelProjection> projections[2];
281  PixelProjection proj;
282  for (size_t k = 0; k < m_prism.unfoldedWindows.size(); k++) {
283  projectPixel(pixel.yc, pixel.Dy, k, dydz, proj); // for even j
284  if (proj.Dy > 0) projections[0].push_back(proj);
285  projectPixel(pixel.yc, pixel.Dy, k, -dydz, proj); // for odd j
286  proj.yc = -proj.yc;
287  if (proj.Dy > 0) projections[1].push_back(proj);
288  }
289  if (projections[0].empty() and projections[1].empty()) continue;
290 
291  for (unsigned k = 0; k < 2; k++) {
292  std::sort(projections[k].begin(), projections[k].end());
293  for (auto& projection : projections[k]) {
294  int iDy = lround(projection.Dy * 1000);
295  auto& mask = masks[iDy];
296  if (not mask) {
297  double Dy = projection.Dy;
298  double step = m_quasyEnergyDistribution->step;
299  mask = new EnergyMask(D.dyB_de, D.dyB_dL, D.dyB_dx, Dy, m_length, pixel.Dx, step);
300  }
301  projection.mask = mask;
302  }
303  }
304 
305  double Ecp_old = 0;
306  double wid_old = 1000;
307  m_results.push_back(Result(pixelID));
308  for (int j = j1; j < j2; j++) {
309  double ybar = j * m_bars.front().B - yB;
310  for (const auto& projection : projections[std::abs(j) % 2]) {
311  double Ecp = (ybar + projection.yc) / D.dyB_de + m_meanE;
312  double wid = projection.mask->getFullWidth();
313  if (std::abs(Ecp - Ecp_old) > (wid + wid_old) / 2 and m_results.back().sum > 0) {
314  m_results.push_back(Result(pixelID));
315  }
316  integrate(projection.mask, Ecp, m_results.back());
317  Ecp_old = Ecp;
318  wid_old = wid;
319  }
320  }
321 
322  if (m_results.back().sum == 0) m_results.pop_back();
323  }
324 
325  for (auto& result : m_results) result.set();
326 
327  for (auto& mask : masks) {
328  if (mask.second) delete mask.second;
329  }
330 
331  }
332 
333 
334  void YScanner::merge(unsigned col, double dydz, int j1, int j2) const
335  {
336  int Neven = func::getNumOfEven(j1, j2);
337  int Nodd = j2 - j1 - Neven;
338 
339  for (unsigned row = 0; row < m_pixelPositions.getNumPixelRows(); row++) {
340 
341  int pixelID = m_pixelPositions.pixelID(row, col);
342  if (not m_pixelMasks.isActive(pixelID)) continue;
343 
344  const auto& pixel = m_pixelPositions.get(pixelID);
345  double Dy0 = 0;
346  double Dy1 = 0;
347  PixelProjection proj;
348  for (size_t k = 0; k < m_prism.unfoldedWindows.size(); k++) {
349  projectPixel(pixel.yc, pixel.Dy, k, dydz, proj); // for even reflections
350  if (proj.Dy > 0) Dy0 += proj.Dy;
351  projectPixel(pixel.yc, pixel.Dy, k, -dydz, proj); // for odd reflections
352  if (proj.Dy > 0) Dy1 += proj.Dy;
353  }
354  if (Dy0 == 0 and Dy1 == 0) continue;
355 
356  double Dy = (Dy0 * Neven + Dy1 * Nodd) / (Neven + Nodd);
357  Result result(pixelID);
358  result.sum = Dy / m_bars.front().B;
359  result.e0 = m_meanE;
360  result.sigsq = m_rmsE * m_rmsE;
361  m_results.push_back(result);
362  }
363  }
364 
365 
366  void YScanner::integrate(const EnergyMask* energyMask, double Ecp, Result& result) const
367  {
368  const auto& mask = energyMask->getMask();
369 
370  if (mask.empty()) {
371  // direct mask calculation
372  for (size_t i = 0; i < m_quasyEnergyDistribution->entries.size(); i++) {
373  double E = m_quasyEnergyDistribution->getX(i);
374  double m = energyMask->getMask(E - Ecp);
375  if (m > 0) {
376  const auto& entry = m_quasyEnergyDistribution->entries[i];
377  double s = entry.y * m;
378  result.sum += s;
379  result.e0 += entry.x * s;
380  result.sigsq += entry.xsq * s;
381  }
382  }
383  } else {
384  // pre-calculated discrete mask w/ linear interpolation
385  int i0 = m_quasyEnergyDistribution->getIndex(Ecp);
386  double fract = -(Ecp - m_quasyEnergyDistribution->getX(i0)) / m_quasyEnergyDistribution->step;
387  if (fract < 0) {
388  i0++;
389  fract += 1;
390  }
391  int N = m_quasyEnergyDistribution->entries.size() - 1;
392  int M = mask.size() - 1;
393  int i1 = std::max(i0 - M, 0);
394  int i2 = std::min(i0 + M - 1, N);
395  for (int i = i1; i <= i2; i++) {
396  const auto& entry = m_quasyEnergyDistribution->entries[i];
397  double m = mask[std::abs(i - i0)] * (1 - fract) + mask[std::abs(i - i0 + 1)] * fract;
398  double s = entry.y * m;
399  result.sum += s;
400  result.e0 += entry.x * s;
401  result.sigsq += entry.xsq * s;
402  }
403  }
404  }
405 
406 
407  double YScanner::prismEntranceY(double y, int k, double dydz) const
408  {
409  const auto& win = m_prism.unfoldedWindows[k];
410  double dz = std::abs(m_prism.zD - m_prism.zFlat);
411  double z1 = y * win.sz + win.z0 + win.nz * dz;
412  double y1 = y * win.sy + win.y0 + win.ny * dz;
413  return y1 + dydz * (m_prism.zR - z1);
414  }
415 
416 
417  void YScanner::projectPixel(double yc, double size, int k, double dydz, PixelProjection& proj) const
418  {
419  double halfSize = (k - m_prism.k0) % 2 == 0 ? size / 2 : -size / 2;
420  double y1 = prismEntranceY(yc - halfSize, k, dydz);
421  double y2 = prismEntranceY(yc + halfSize, k, dydz);
422 
423  double Bh = m_bars.front().B / 2;
424  y1 = std::max(y1, -Bh);
425  y2 = std::min(y2, Bh);
426 
427  proj.yc = (y1 + y2) / 2;
428  proj.Dy = y2 - y1;
429  }
430 
431  } //TOP
433 } //Belle2
R E
internal precision of FFTW codelets
A mask for energy masking.
Definition: EnergyMask.h:24
const std::vector< double > & getMask() const
Returns discrete mask (note: only half of the mask is stored)
Definition: EnergyMask.h:68
Pixel relative efficiencies of a single module.
Pixel masks of a single module.
Definition: PixelMasks.h:22
bool isActive(int pixelID) const
Checks if pixel is active.
Definition: PixelMasks.h:85
Pixel positions and dimensions in module local frame.
int pixelID(unsigned row, unsigned col) const
Transforms pixel row and column to pixel ID Note: for convenience pixel row and column numbering star...
unsigned getNumPixelRows() const
Returns the number of pixel rows.
const PixelData & get(int pixelID) const
Returns pixel data for given pixelID.
Base class with geometry data.
Definition: RaytracerBase.h:27
@ c_Unified
single bar with average width and thickness
Definition: RaytracerBase.h:34
@ c_SemiLinear
semi-linear approximation
Definition: RaytracerBase.h:42
Prism m_prism
prism geometry data
std::vector< BarSegment > m_bars
geometry data of bar segments
static TOPGeometryPar * Instance()
Static method to obtain the pointer to its instance.
static const double c_hc
Planck constant times speed of light in [eV*nm].
double m_cosTotal
cosine of total reflection angle
Definition: YScanner.h:468
void prepare(double momentum, double beta, double length) const
Prepare for the PDF expansion in y for a given track mass hypothesis.
Definition: YScanner.cc:118
Table m_energyDistribution
photon energy distribution
Definition: YScanner.h:479
PixelEfficiencies m_pixelEfficiencies
pixel relative efficiencies
Definition: YScanner.h:464
double m_beta
particle beta
Definition: YScanner.h:472
bool m_scanDone
true if scan performed, false if reflections just merged
Definition: YScanner.h:487
void projectPixel(double yc, double size, int k, double dydz, PixelProjection &proj) const
Calculates a projection of a pixel to prism entrance window (going down-stream the photon).
Definition: YScanner.cc:417
Table m_efficiency
nominal photon detection efficiencies (PDE)
Definition: YScanner.h:465
double m_meanE
mean photon energy
Definition: YScanner.h:475
static int s_maxReflections
maximal number of reflections to perform scan
Definition: YScanner.h:489
double m_rmsE
r.m.s of photon energy
Definition: YScanner.h:476
bool m_aboveThreshold
true if beta is above the Cerenkov threshold
Definition: YScanner.h:483
double m_sigmaScat
r.m.s.
Definition: YScanner.h:477
void merge(unsigned col, double dydz, int j1, int j2) const
Performs expansion by merging all reflections.
Definition: YScanner.cc:334
void scan(unsigned col, double yB, double dydz, const Derivatives &D, int j1, int j2) const
Performs expansion w/ the scan over reflections.
Definition: YScanner.cc:270
double m_meanE0
mean photon energy for beta = 1
Definition: YScanner.h:466
double m_momentum
particle momentum magnitude
Definition: YScanner.h:471
Table * m_quasyEnergyDistribution
a pointer to the element in m_quasyEnergyDistributions
Definition: YScanner.h:482
PixelMasks m_pixelMasks
pixel masks
Definition: YScanner.h:463
void setQuasyEnergyDistribution(double sigma) const
Sets photon energy distribution convoluted with a normalized Gaussian.
Definition: YScanner.cc:188
void integrate(const EnergyMask *energyMask, double Ecp, Result &result) const
Integrates quasy energy distribution multiplied with energy mask.
Definition: YScanner.cc:366
void expand(unsigned col, double yB, double dydz, const Derivatives &D, int Ny, bool doScan) const
Performs the PDF expansion in y for a given pixel column using scan or merge methods.
Definition: YScanner.cc:234
YScanner(int moduleID, unsigned N=64)
Class constructor.
Definition: YScanner.cc:47
std::map< int, Table > m_quasyEnergyDistributions
photon energy distributions convoluted with Gaussian of different widths
Definition: YScanner.h:481
double prismEntranceY(double y, int k, double dydz) const
Calculates y at prism entrance from detection position, reflection number and photon slope.
Definition: YScanner.cc:407
double setEnergyDistribution(double beta) const
Sets photon energy distribution and mean photon energy according to nominal PDE and particle beta.
Definition: YScanner.cc:159
double m_length
length of particle trajectory inside quartz
Definition: YScanner.h:473
double m_rmsE0
r.m.s of photon energy for beta = 1
Definition: YScanner.h:467
std::vector< Result > m_results
results of PDF expansion in y
Definition: YScanner.h:486
PixelPositions m_pixelPositions
positions and sizes of pixels
Definition: YScanner.h:462
void clear() const
Clear mutable variables.
Definition: YScanner.cc:99
double m_numPhotons
number of photons per Cerenkov azimuthal angle per track length
Definition: YScanner.h:474
double m_sigmaAlpha
surface roughness parameter in photon energy units
Definition: YScanner.h:478
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
int getNumOfEven(int j1, int j2)
Returns number of even numbers in the range given by arguments.
Definition: func.h:92
Abstract base class for different kinds of events.
Solution of inverse ray-tracing.
double step
step for numerical derivative calculation
std::vector< TOPGeoPrism::UnfoldedWindow > unfoldedWindows
unfolded prism exit windows
double zFlat
z where flat continues to slanted surface
double zR
maximal z, i.e position of prism-bar joint
double zD
detector (photo-cathode) position
int k0
index of true prism in the vector 'unfoldedWindows'
static double dyB_d(const InverseRaytracer::Solution &sol0, const InverseRaytracer::Solution &sol1)
Calculates the derivative of unfolded y coordinate at prism entrance.
Definition: YScanner.h:504
static double dLen_d(const InverseRaytracer::Solution &sol0, const InverseRaytracer::Solution &sol1)
Calculates the derivative of propagation length.
Definition: YScanner.h:498
double dFic_de
Cerenkov azimuthal angle over photon energy.
Definition: YScanner.h:48
double dLen_de
propagation length over photon energy
Definition: YScanner.h:42
Derivatives()
Default constructor.
Definition: YScanner.h:54
double dFic_dx
Cerenkov azimuthal angle over photon detection coordinate x.
Definition: YScanner.h:47
double dFic_dL
Cerenkov azimuthal angle over running parameter of particle trajectory.
Definition: YScanner.h:49
double dLen_dL
propagation length over running parameter of particle trajectory
Definition: YScanner.h:43
double dyB_de
unfolded y coordinate at prism entrance over photon energy
Definition: YScanner.h:45
double dLen_dx
propagation length over photon detection coordinate x
Definition: YScanner.h:41
double dyB_dL
unfolded y coordinate at prism entrance over running parameter of particle trajectory
Definition: YScanner.h:46
double dyB_dx
unfolded y coordinate at prism entrance over photon detection coordinate x
Definition: YScanner.h:44
static double dFic_d(const InverseRaytracer::Solution &sol0, const InverseRaytracer::Solution &sol1)
Calculates the derivative of Cerenkov azimuthal angle.
Definition: YScanner.h:510
Down-stream projection of a pixel to prism entrance window w/ a clip on bar exit thickness.
Definition: YScanner.h:182
double yc
center in y of clipped pixel projection
Definition: YScanner.h:183
double Dy
size in y of clipped pixel projection
Definition: YScanner.h:184
Single PDF peak data.
Definition: YScanner.h:195
double step
step size
Definition: YScanner.h:120
double getY(int i) const
Returns y for a given index.
Definition: YScanner.h:549
int getIndex(double x) const
Returns index.
Definition: YScanner.h:547
std::vector< TableEntry > entries
table entries
Definition: YScanner.h:121
double getXmax() const
Returns x of the last entry.
Definition: YScanner.h:545
double getX(int i) const
Returns x for a given index.
Definition: YScanner.h:541
void clear()
Clear the content entirely.
Definition: YScanner.h:520
double getXmin() const
Returns x of the first entry.
Definition: YScanner.h:543
void set(double X0, double Step)
Sets the first x and the step, and clears the entries.
Definition: YScanner.h:527