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>
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");
53 B2FATAL(
"TOP::YScanner: N must be > 1");
60 const auto* geo = topgp->getGeometry();
61 auto qe = geo->getNominalQE();
62 qe.applyFilterTransmission(geo->getWavelengthFilter());
67 B2FATAL(
"TOP::YScanner: quantum efficiency found zero for all wavelengths");
72 const auto& tdc = geo->getNominalTDC();
73 for (
unsigned i = 0; i < N; i++) {
76 double effi = qe.getEfficiency(lambda) * tdc.getEfficiency();
87 double p = entry.y * (1 - 1 / pow(topgp->getPhaseIndex(e), 2));
129 if (area == 0)
return;
137 const double radLength = 12.3;
138 double thetaScat = 13.6e-3 / beta / momentum *
sqrt(length / 2 / radLength);
141 double n = topgp->getPhaseIndex(
m_meanE);
143 B2ERROR(
"TOP::YScanner::prepare: beta * n < 1 ==> must be a bug!");
146 double dndE = topgp->getPhaseIndexDerivative(
m_meanE);
147 double dEdTheta = n *
sqrt(pow(beta * n, 2) - 1) / dndE;
170 double p = std::max(entry.y * (1 - 1 / pow(beta * topgp->getPhaseIndex(e), 2)), 0.0);
171 double ee = entry.xsq;
177 if (s == 0)
return 0;
192 B2ERROR(
"TOP::YScanner:setQuasyEnergyDistribution: unexpectedly large size of the std::map found, map cleared");
196 int ng = lround(3 * sigma / step);
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));
209 for (
int k = -ng; k < N + ng; k++) {
213 for (
int i = -ng; i <= ng; i++) {
224 quasyEnergyDistribution.entries.push_back(
TableEntry(s, se, see));
227 for (
auto& entry : quasyEnergyDistribution.entries) entry.y /= sum;
238 if (D.dyB_de == 0)
return;
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;
250 y1 += D.dyB_de * (minE -
m_meanE);
251 y2 += D.dyB_de * (maxE -
m_meanE);
253 y1 += D.dyB_de * (maxE -
m_meanE);
254 y2 += D.dyB_de * (minE -
m_meanE);
256 double B =
m_bars.front().B;
257 int j1 = lround(y1 / B);
258 int j2 = lround(y2 / B) + 1;
261 scan(col, yB, dydz, D, j1, j2);
264 merge(col, dydz, j1, j2);
273 std::map<int, EnergyMask*> masks;
280 std::vector<PixelProjection> projections[2];
284 if (proj.
Dy > 0) projections[0].push_back(proj);
287 if (proj.
Dy > 0) projections[1].push_back(proj);
289 if (projections[0].empty() and projections[1].empty())
continue;
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];
297 double Dy = projection.Dy;
301 projection.mask = mask;
306 double wid_old = 1000;
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) {
325 for (
auto& result :
m_results) result.set();
327 for (
auto& mask : masks) {
328 if (mask.second)
delete mask.second;
337 int Nodd = j2 - j1 - Neven;
350 if (proj.
Dy > 0) Dy0 += proj.
Dy;
352 if (proj.
Dy > 0) Dy1 += proj.
Dy;
354 if (Dy0 == 0 and Dy1 == 0)
continue;
356 double Dy = (Dy0 * Neven + Dy1 * Nodd) / (Neven + Nodd);
358 result.sum = Dy /
m_bars.front().B;
368 const auto& mask = energyMask->
getMask();
374 double m = energyMask->
getMask(
E - Ecp);
377 double s = entry.y * m;
379 result.e0 += entry.x * s;
380 result.sigsq += entry.xsq * s;
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++) {
397 double m = mask[std::abs(i - i0)] * (1 - fract) + mask[std::abs(i - i0 + 1)] * fract;
398 double s = entry.y * m;
400 result.e0 += entry.x * s;
401 result.sigsq += entry.xsq * s;
411 double z1 = y * win.sz + win.z0 + win.nz * dz;
412 double y1 = y * win.sy + win.y0 + win.ny * dz;
419 double halfSize = (k -
m_prism.
k0) % 2 == 0 ? size / 2 : -size / 2;
423 double Bh =
m_bars.front().B / 2;
424 y1 = std::max(y1, -Bh);
425 y2 = std::min(y2, Bh);
427 proj.
yc = (y1 + y2) / 2;
A mask for energy masking.
const std::vector< double > & getMask() const
Returns discrete mask (note: only half of the mask is stored)
Pixel relative efficiencies of a single module.
Pixel masks of a single module.
bool isActive(int pixelID) const
Checks if pixel is active.
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.
@ c_Unified
single bar with average width and thickness
@ c_SemiLinear
semi-linear approximation
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
void prepare(double momentum, double beta, double length) const
Prepare for the PDF expansion in y for a given track mass hypothesis.
Table m_energyDistribution
photon energy distribution
PixelEfficiencies m_pixelEfficiencies
pixel relative efficiencies
double m_beta
particle beta
bool m_scanDone
true if scan performed, false if reflections just merged
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).
Table m_efficiency
nominal photon detection efficiencies (PDE)
double m_meanE
mean photon energy
static int s_maxReflections
maximal number of reflections to perform scan
double m_rmsE
r.m.s of photon energy
bool m_aboveThreshold
true if beta is above the Cerenkov threshold
void merge(unsigned col, double dydz, int j1, int j2) const
Performs expansion by merging all reflections.
void scan(unsigned col, double yB, double dydz, const Derivatives &D, int j1, int j2) const
Performs expansion w/ the scan over reflections.
double m_meanE0
mean photon energy for beta = 1
double m_momentum
particle momentum magnitude
Table * m_quasyEnergyDistribution
a pointer to the element in m_quasyEnergyDistributions
PixelMasks m_pixelMasks
pixel masks
void setQuasyEnergyDistribution(double sigma) const
Sets photon energy distribution convoluted with a normalized Gaussian.
void integrate(const EnergyMask *energyMask, double Ecp, Result &result) const
Integrates quasy energy distribution multiplied with energy mask.
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.
YScanner(int moduleID, unsigned N=64)
Class constructor.
std::map< int, Table > m_quasyEnergyDistributions
photon energy distributions convoluted with Gaussian of different widths
double prismEntranceY(double y, int k, double dydz) const
Calculates y at prism entrance from detection position, reflection number and photon slope.
double setEnergyDistribution(double beta) const
Sets photon energy distribution and mean photon energy according to nominal PDE and particle beta.
double m_length
length of particle trajectory inside quartz
double m_rmsE0
r.m.s of photon energy for beta = 1
std::vector< Result > m_results
results of PDF expansion in y
PixelPositions m_pixelPositions
positions and sizes of pixels
void clear() const
Clear mutable variables.
double m_numPhotons
number of photons per Cerenkov azimuthal angle per track length
double m_sigmaAlpha
surface roughness parameter in photon energy units
double sqrt(double a)
sqrt for double
int getNumOfEven(int j1, int j2)
Returns number of even numbers in the range given by arguments.
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.
static double dLen_d(const InverseRaytracer::Solution &sol0, const InverseRaytracer::Solution &sol1)
Calculates the derivative of propagation length.
double dFic_de
Cerenkov azimuthal angle over photon energy.
double dLen_de
propagation length over photon energy
Derivatives()
Default constructor.
double dFic_dx
Cerenkov azimuthal angle over photon detection coordinate x.
double dFic_dL
Cerenkov azimuthal angle over running parameter of particle trajectory.
double dLen_dL
propagation length over running parameter of particle trajectory
double dyB_de
unfolded y coordinate at prism entrance over photon energy
double dLen_dx
propagation length over photon detection coordinate x
double dyB_dL
unfolded y coordinate at prism entrance over running parameter of particle trajectory
double dyB_dx
unfolded y coordinate at prism entrance over photon detection coordinate x
static double dFic_d(const InverseRaytracer::Solution &sol0, const InverseRaytracer::Solution &sol1)
Calculates the derivative of Cerenkov azimuthal angle.
Down-stream projection of a pixel to prism entrance window w/ a clip on bar exit thickness.
double yc
center in y of clipped pixel projection
double Dy
size in y of clipped pixel projection
double getY(int i) const
Returns y for a given index.
int getIndex(double x) const
Returns index.
std::vector< TableEntry > entries
table entries
double getXmax() const
Returns x of the last entry.
double getX(int i) const
Returns x for a given index.
void clear()
Clear the content entirely.
double getXmin() const
Returns x of the first entry.
void set(double X0, double Step)
Sets the first x and the step, and clears the entries.