9#include <top/reconstruction_cpp/PDFConstructor.h>
10#include <top/reconstruction_cpp/TOPRecoManager.h>
11#include <top/reconstruction_cpp/func.h>
12#include <top/geometry/TOPGeometryPar.h>
13#include <framework/logging/Logger.h>
37 if (not track.isValid()) {
38 B2ERROR(
"TOP::PDFConstructor: TOPTrack is not valid, cannot continue");
44 B2ERROR(
"TOP::PDFConstructor: missing reconstruction objects, cannot continue");
48 m_beta = track.getBeta(hypothesis, overrideMass);
49 m_yScanner->prepare(track.getMomentumMag(),
m_beta, track.getLengthInQuartz());
50 m_tof = track.getTOF(hypothesis, 0, overrideMass);
62 const auto& pixelPositions =
m_yScanner->getPixelPositions();
63 int numPixels = pixelPositions.getNumPixels();
65 for (
int pixelID = 1; pixelID <= numPixels; pixelID++) {
66 auto pmtType = pixelPositions.get(pixelID).pmtType;
67 const auto& tts = geo->getTTS(pmtType);
94 if (
m_track.getEmissionPoint().position.Z() > prism.zR) {
119 const auto& pixelPositions =
m_yScanner->getPixelPositions();
125 double xmi = 0, xma = 0;
126 bool ok =
rangeOfX(prism.zD, xmi, xma);
128 int kmi = lround(xmi / bar.A);
129 int kma = lround(xma / bar.A);
133 for (
int k = kmi; k <= kma; k++) {
134 for (
unsigned col = 0; col < pixelPositions.getNumPixelColumns(); col++) {
135 const auto& pixel = pixelPositions.get(col + 1);
136 if (pixel.Dx == 0)
continue;
138 if (xD < xmi or xD > xma)
continue;
150 const auto& emiPoint =
m_track.getEmissionPoint().position;
154 double xmi = 0, xma = 0;
155 bool ok =
rangeOfX(mirror.zb, xmi, xma);
157 int kmi = lround(xmi / bar.A);
158 int kma = lround(xma / bar.A);
162 double xE = emiPoint.X();
163 double zE = emiPoint.Z();
164 double Ah = bar.A / 2;
165 for (
int k = kmi; k <= kma; k++) {
168 double xL =
func::clip(-Ah, k, bar.A, xmi, xma);
169 double xR =
func::clip(Ah, k, bar.A, xmi, xma);
178 const auto& pixelPositions =
m_yScanner->getPixelPositions();
184 std::vector<double> xDs;
185 double minLen = 1e10;
188 if (xDs.size() < 2)
return;
193 std::sort(xDs.begin(), xDs.end());
194 double xmi = xDs.front();
195 double xma = xDs.back();
197 int kmi = lround(xmi / bar.A);
198 int kma = lround(xma / bar.A);
202 for (
int k = kmi; k <= kma; k++) {
203 for (
unsigned col = 0; col < pixelPositions.getNumPixelColumns(); col++) {
204 const auto& pixel = pixelPositions.get(col + 1);
205 if (pixel.Dx == 0)
continue;
207 if (xD < xmi or xD > xma)
continue;
209 setSignalPDF(reflected, col, xD, prism.zD, Nxm, xmMin, xmMax);
220 if (i0 < 0)
return false;
223 for (
unsigned i = 0; i < 2; i++) {
226 const auto& sol = solutions[i0];
227 xDs.push_back(sol.xD);
228 minLen = std::min(minLen, sol.len);
238 const double precision = 0.01;
241 double y1 =
deltaXD(x1, sol, xD);
242 if (isnan(y1))
return false;
246 double step = -dFic_dx * y1;
247 for (
int i = 1; i < 20; i++) {
248 double x2 = step * i;
249 double y2 =
deltaXD(x2, sol, xD);
250 if (isnan(y2))
return false;
256 for (
int k = 0; k < 20; k++) {
257 double x = (x2 + x3) / 2;
258 double y =
deltaXD(x, sol, xD);
259 if (isnan(y))
return false;
268 if (y2 * y1 > 0)
return false;
272 for (
int k = 0; k < 20; k++) {
273 double x = (x1 + x2) / 2;
274 double y =
deltaXD(x, sol, xD);
275 if (isnan(y))
return false;
300 for (
int i = 0; i < 10; i++) {
305 if (not ok)
return false;
308 for (
int i = 0; i < 10; i++) {
313 if (not ok)
return false;
316 for (
int i = 0; i < 10; i++) {
317 ok =
raytrace(rayTracer_dFic, 0, 0, dFic);
321 if (not ok)
return false;
325 double dLen_dL = (rayTracer_dL.getPropagationLen() -
m_fastRaytracer->getPropagationLen()) / dL;
326 double dLen_de = (rayTracer_de.getPropagationLen() -
m_fastRaytracer->getPropagationLen()) / de;
327 double dLen_dFic = (rayTracer_dFic.getPropagationLen() -
m_fastRaytracer->getPropagationLen()) / dFic;
329 double dyB_dL = (rayTracer_dL.getYB() -
m_fastRaytracer->getYB()) / dL;
330 double dyB_de = (rayTracer_de.getYB() -
m_fastRaytracer->getYB()) / de;
331 double dyB_dFic = (rayTracer_dFic.getYB() -
m_fastRaytracer->getYB()) / dFic;
335 double dx_dFic = (rayTracer_dFic.getXD() -
m_fastRaytracer->getXD()) / dFic;
337 if (dx_dFic == 0)
return false;
341 D.dLen_dx = dLen_dFic / dx_dFic;
342 D.dLen_de = dLen_de - dLen_dFic * dx_de / dx_dFic;
343 D.dLen_dL = dLen_dL - dLen_dFic * dx_dL / dx_dFic;
345 D.dyB_dx = dyB_dFic / dx_dFic;
346 D.dyB_de = dyB_de - dyB_dFic * dx_de / dx_dFic;
347 D.dyB_dL = dyB_dL - dyB_dFic * dx_dL / dx_dFic;
349 D.dFic_dx = 1 / dx_dFic;
350 D.dFic_de = - dx_de / dx_dFic;
351 D.dFic_dL = - dx_dL / dx_dFic;
359 const auto& emi =
m_track.getEmissionPoint(dL);
360 const auto& trk = emi.trackAngles;
363 double fic =
m_Fic + dFic;
364 double cosFic = cos(fic);
365 double sinFic = sin(fic);
366 double a = trk.cosTh * cer.sinThc * cosFic + trk.sinTh * cer.cosThc;
367 double b = cer.sinThc * sinFic;
368 double kx = a * trk.cosFi - b * trk.sinFi;
369 double ky = a * trk.sinFi + b * trk.cosFi;
370 double kz = trk.cosTh * cer.cosThc - trk.sinTh * cer.sinThc * cosFic;
399 double dTime =
m_fastRaytracer->getPropagationLenDelta() / speedOfLightQuartz;
403 double dt_dx = D.dLen_dx / speedOfLightQuartz;
407 double sigmaScat = D.dLen_de *
m_yScanner->getSigmaScattering() / speedOfLightQuartz;
410 double sigmaAlpha = D.dLen_de *
m_yScanner->getSigmaAlpha() / speedOfLightQuartz;
413 const auto& pixel =
m_yScanner->getPixelPositions().get(col + 1);
414 double L =
m_track.getLengthInQuartz();
417 double wid0 = (pow(dt_dx * pixel.Dx, 2) + pow(dt_dL * L, 2) + pow(dTime, 2)) / 12 + pow(sigmaScat, 2) +
418 pow(sigmaAlpha, 2) * Ny_eff;
421 double wid = wid0 + pow(dt_de *
m_yScanner->getRMSEnergy(), 2);
425 const auto& atPrismEntrance = photonStates[photonStates.size() - 2];
426 double dydz = atPrismEntrance.getKy() / atPrismEntrance.getKz();
431 double time =
m_tof + Len / speedOfLightQuartz;
432 doScan =
m_track.isScanRequired(col, time, wid);
435 m_yScanner->expand(col, yB, dydz, D, Ny_eff, doScan);
437 double numPhotons =
m_yScanner->getNumPhotons() * std::abs(D.dFic_dx * pixel.Dx);
440 for (
const auto& result :
m_yScanner->getResults()) {
441 double RQE =
m_yScanner->getPixelEfficiencies().get(result.pixelID);
442 if (RQE == 0)
continue;
444 double dE = result.e0 -
m_yScanner->getMeanEnergy();
445 double propLen = Len + D.dLen_de * dE;
449 peak.
t0 =
m_tof + propLen / speedOfLight;
450 peak.
wid = wid0 + dt_de * dt_de * result.sigsq;
453 signalPDF.append(peak);
460 extra.
sige = result.sigsq;
471 const auto& firstState = photonStates.front();
472 extra.
kxE = firstState.getKx();
473 extra.
kyE = firstState.getKy();
474 extra.
kzE = firstState.getKz();
475 const auto& lastState = photonStates.back();
476 extra.
kxD = lastState.getKx();
477 extra.
kyD = lastState.getKy();
478 extra.
kzD = lastState.getKz();
480 signalPDF.append(extra);
489 double surf =
m_yScanner->getBars().front().reflectivity;
490 double p = exp(-propLen / bulk) * pow(surf, std::abs(nx) + std::abs(ny));
498 if (maxLen < 0)
return false;
500 const auto& emission =
m_track.getEmissionPoint();
501 const auto& trk = emission.trackAngles;
506 double dz = z - emission.position.Z();
507 double cosFicLimit = (trk.cosTh * cer.cosThc - dz / maxLen) / (trk.sinTh * cer.sinThc);
508 double cosLimit = (dz > 0) ? cosFicLimit : -cosFicLimit;
509 if (cosLimit < -1)
return false;
511 std::vector<double> xmima;
512 double x0 = emission.position.X();
514 xmima.push_back(x0 - maxLen);
515 xmima.push_back(x0 + maxLen);
517 double a = trk.cosTh * cer.sinThc * cosFicLimit + trk.sinTh * cer.cosThc;
518 double b = cer.sinThc *
sqrt(1 - cosFicLimit * cosFicLimit);
519 xmima.push_back(x0 + maxLen * (a * trk.cosFi - b * trk.sinFi));
520 xmima.push_back(x0 + maxLen * (a * trk.cosFi + b * trk.sinFi));
521 std::sort(xmima.begin(), xmima.end());
528 double theta = acos(trk.cosTh);
529 if (dz < 0) theta = M_PI - theta;
531 double thetaCer = acos(cer.cosThc);
532 if (theta - thetaCer >= M_PI / 2)
return false;
534 std::vector<double> dxdz;
535 double a = -cos(theta + thetaCer) * cos(theta - thetaCer);
536 double b = sin(2 * theta) * trk.cosFi;
537 double c = pow(trk.sinFi * cer.sinThc, 2) - pow(trk.cosFi, 2) * sin(theta + thetaCer) * sin(theta - thetaCer);
538 double D = b * b - 4 * a * c;
539 if (D < 0)
return true;
542 dxdz.push_back((-b - D) / 2 / a);
543 dxdz.push_back((-b + D) / 2 / a);
545 if (b == 0)
return true;
546 dxdz.push_back(-c / b);
547 dxdz.push_back(copysign(INFINITY, b));
549 std::vector<double> cosFic(2, cosLimit);
550 for (
int i = 0; i < 2; i++) {
551 if (std::abs(dxdz[i]) < INFINITY) {
552 double aa = (dxdz[i] * cos(theta) - trk.cosFi * sin(theta)) * cer.cosThc;
553 double bb = (dxdz[i] * sin(theta) + trk.cosFi * cos(theta)) * cer.sinThc;
554 double dd = trk.sinFi * cer.sinThc;
555 cosFic[i] = aa * bb / (bb * bb + dd * dd);
556 double kz = cos(theta) * cer.cosThc - sin(theta) * cer.sinThc * cosFic[i];
557 if (kz < 0) dxdz[i] = copysign(INFINITY, dxdz[1 - i] - dxdz[i]);
560 if (dxdz[0] > dxdz[1]) {
561 std::reverse(dxdz.begin(), dxdz.end());
562 std::reverse(cosFic.begin(), cosFic.end());
564 for (
int i = 0; i < 2; i++) {
565 if (cosFic[i] < cosLimit) xmima[i] = x0 + dxdz[i] * dz;
569 xmi = std::max(xmima[0], x0 - maxLen);
570 xma = std::min(xmima[1], x0 + maxLen);
578 double z =
sqrt(1 - x * x);
579 double kx = (x - xe);
580 double kz = (z - ze);
581 double s = 2 * (kx * x + kz * z);
582 double qx = kx - s * x;
583 double qz = kz - s * z;
585 double der_z = -x / z;
586 double der_s = 2 * (kx + der_z * kz);
587 double der_qx = (1 - s) - der_s * x;
588 double der_qz = (1 - s) * der_z - der_s * z;
590 return 1 - der_z * qx / qz + (zd - z) * (der_qx * qz - der_qz * qx) / (qz * qz);
604 double xe = (xE - mirror.
xc) / mirror.
R;
605 double ze = (zE - mirror.
zc) / mirror.
R;
606 double zd = (zD - mirror.
zc) / mirror.
R;
610 double x1 = (-Ah - mirror.
xc) / mirror.
R;
612 if (y1 != y1 or std::abs(y1) == INFINITY)
return -Ah;
614 double x2 = (Ah - mirror.
xc) / mirror.
R;
616 if (y2 != y2 or std::abs(y2) == INFINITY)
return -Ah;
618 if (y1 * y2 > 0)
return -Ah;
620 for (
int i = 0; i < 50; i++) {
621 double x = (x1 + x2) / 2;
623 if (y != y or std::abs(y) == INFINITY)
return -Ah;
631 double x = (x1 + x2) / 2;
633 return x * mirror.
R + mirror.
xc;
640 const auto& pixelPositions =
m_yScanner->getPixelPositions();
644 double xE =
m_track.getEmissionPoint().position.X();
645 int nxmi = (xE > 0) ? 0 : -1;
646 int nxma = (xE > 0) ? 1 : 0;
648 for (
const auto& pixel : pixelPositions.getPixels()) {
649 if (not
m_yScanner->getPixelMasks().isActive(pixel.ID))
continue;
650 double RQE =
m_yScanner->getPixelEfficiencies().get(pixel.ID);
651 if (RQE == 0)
continue;
653 for (
int Nxe = nxmi; Nxe <= nxma; Nxe++) {
654 for (
size_t k = 0; k < prism.unfoldedWindows.size(); k++) {
656 if (sol.len == 0 or std::abs(sol.L) >
m_track.getLengthInQuartz() / 2)
continue;
659 if (not ok)
continue;
660 int Nye = k - prism.k0;
666 double slope = lastState.getKy() / lastState.getKz();
667 double dz = prism.zD - prism.zFlat;
668 double y2 = std::min(pixel.yc + pixel.Dy / 2, prism.yUp + slope * dz);
669 double y1 = std::max(pixel.yc - pixel.Dy / 2, prism.yDown + slope * dz);
671 if (Dy < 0)
continue;
674 for (
int i = 0; i < 4; i++) {
680 if (not ok)
continue;
684 for (
int i = 0; i < 4; i++) {
690 if (not ok)
continue;
694 for (
int i = 0; i < 4; i++) {
700 if (not ok)
continue;
703 double dx_dL = (lastState_dL.getX() - lastState.getX()) / dL;
704 double dy_dL = (lastState_dL.getY() - lastState.getY()) / dL;
705 double dx_dFic = (lastState_dFic.getX() - lastState.getX()) / dFic;
706 double dy_dFic = (lastState_dFic.getY() - lastState.getY()) / dFic;
707 double Jacobi = dx_dL * dy_dFic - dy_dL * dx_dFic;
708 double numPhotons =
m_yScanner->getNumPhotonsPerLen() * pixel.Dx * Dy / std::abs(Jacobi) * RQE;
710 double dLen_de = (lastState_de.getPropagationLen() - lastState.getPropagationLen()) / de;
711 double dLen_dL = (lastState_dL.getPropagationLen() - lastState.getPropagationLen()) / dL;
712 double dLen_dFic = (lastState_dFic.getPropagationLen() - lastState.getPropagationLen()) / dFic;
715 double dt_dL = dLen_dL / speedOfLightQuartz;
716 double dt_dFic = dLen_dFic / speedOfLightQuartz;
718 double chromatic = pow(dt_de *
m_yScanner->getRMSEnergy(), 2);
720 double DL = (dy_dFic * pixel.Dx - dx_dFic * pixel.Dy) / Jacobi;
721 double DFic = (dx_dL * pixel.Dy - dy_dL * pixel.Dx) / Jacobi;
722 double paralax = (pow(dt_dL * DL, 2) + pow(dt_dFic * DFic, 2)) / 12;
724 double scattering = pow(dLen_de *
m_yScanner->getSigmaScattering() / speedOfLightQuartz, 2);
728 peak.
wid = chromatic + paralax + scattering;
729 peak.
nph = 1 - exp(-numPhotons);
730 peak.
fic = atan2(sol.sinFic, sol.cosFic);
731 signalPDF.append(peak);
741 extra.
xD = lastState.getXD();
742 extra.
yD = lastState.getYD();
743 extra.
zD = lastState.getZD();
744 extra.
kxE = firstState.getKx();
745 extra.
kyE = firstState.getKy();
746 extra.
kzE = firstState.getKz();
747 extra.
kxD = lastState.getKx();
748 extra.
kyD = lastState.getKy();
749 extra.
kzD = lastState.getKz();
751 signalPDF.append(extra);
762 const auto& emi =
m_track.getEmissionPoint(sol.
L + dL);
763 const auto& trk = emi.trackAngles;
772 double cosFic = sol.
cosFic * cosDFic - sol.
sinFic * sinDFic;
773 double sinFic = sol.
sinFic * cosDFic + sol.
cosFic * sinDFic;
774 double a = trk.cosTh * cer.sinThc * cosFic + trk.sinTh * cer.cosThc;
775 double b = cer.sinThc * sinFic;
776 double kx = a * trk.cosFi - b * trk.sinFi;
777 double ky = a * trk.sinFi + b * trk.cosFi;
778 double kz = trk.cosTh * cer.cosThc - trk.sinTh * cer.sinThc * cosFic;
790 const auto& win = prism.unfoldedWindows[k];
791 double dz = std::abs(prism.zD - prism.zFlat);
793 pixel.
yc * win.sy + win.y0 + win.ny * dz,
794 pixel.
yc * win.sz + win.z0 + win.nz * dz);
797 for (
int iter = 0; iter < 100; iter++) {
799 if (std::abs(sol.L) >
m_track.getLengthInQuartz() / 2)
return sol;
800 if (std::abs(sol.L - L) < 0.01)
return sol;
803 B2DEBUG(20,
"TOP::PDFConstructor::prismSolution: iterations not converging");
810 const auto& emi =
m_track.getEmissionPoint(L);
814 auto r = rD - emi.position;
815 const auto& trk = emi.trackAngles;
816 double xx = r.X() * trk.cosFi + r.Y() * trk.sinFi;
817 double y = -r.X() * trk.sinFi + r.Y() * trk.cosFi;
818 double x = xx * trk.cosTh - r.Z() * trk.sinTh;
819 double z = xx * trk.sinTh + r.Z() * trk.cosTh;
823 double rho =
sqrt(x * x + y * y);
827 sol.
len = rho / cer.sinThc;
828 sol.
L = L + z - sol.
len * cer.cosThc;
840 B2ERROR(
"TOP::PDFConstructor::getLogL(): object status is invalid - cannot provide log likelihood");
847 double f =
pdfValue(hit.pixelID, hit.time, hit.timeErr);
851 B2ERROR(
"TOP::PDFConstructor::getLogL(): PDF value is zero or negative"
853 <<
LogVar(
"pixelID", hit.pixelID) <<
LogVar(
"time", hit.time) <<
LogVar(
"PDFValue", f));
868 B2ERROR(
"TOP::PDFConstructor::getLogL(): object status is invalid - cannot provide log likelihood");
874 if (hit.time < minTime or hit.time > maxTime)
continue;
875 double f =
pdfValue(hit.pixelID, hit.time - t0, hit.timeErr, sigt);
879 B2ERROR(
"TOP::PDFConstructor::getLogL(): PDF value is zero or negative"
881 <<
LogVar(
"pixelID", hit.pixelID) <<
LogVar(
"time", hit.time) <<
LogVar(
"PDFValue", f));
896 B2ERROR(
"TOP::PDFConstructor::getBackgroundLogL(): object status is invalid - cannot provide log likelihood");
904 if (hit.time < minTime or hit.time > maxTime)
continue;
909 B2ERROR(
"TOP::PDFConstructor::getBackgroundLogL(): PDF value is zero or negative"
911 <<
LogVar(
"pixelID", hit.pixelID) <<
LogVar(
"time", hit.time) <<
LogVar(
"PDFValue", f));
922 const std::vector<PDFConstructor::LogL>&
926 B2ERROR(
"TOP::PDFConstructor::getPixelLogLs(): object status is invalid - cannot provide log likelihoods");
933 if (hit.time < minTime or hit.time > maxTime)
continue;
934 double f =
pdfValue(hit.pixelID, hit.time - t0, hit.timeErr, sigt);
938 B2ERROR(
"TOP::PDFConstructor::getPixelLogLs(): PDF value is zero or negative"
940 <<
LogVar(
"pixelID", hit.pixelID) <<
LogVar(
"time", hit.time) <<
LogVar(
"PDFValue", f));
944 unsigned k = hit.pixelID - 1;
948 LL.effectiveSignalYield +=
m_f0 / f;
960 for (
const auto* other :
m_pdfOtherTracks) bfot += other->getExpectedDeltaPhotons(minTime, maxTime);
964 unsigned k = signalPDF.getPixelID() - 1;
965 double phot = signalPDF.getIntegral(minTime, maxTime) *
m_signalPhotons + bfot * pixelPDF[k];
967 const auto& otherPDFs = other->getSignalPDF();
968 phot += otherPDFs[k].getIntegral(minTime, maxTime) * other->getExpectedSignalPhotons();
988 unsigned k = hit.pixelID - 1;
994 double wid0 = hit.timeErr * hit.timeErr;
998 for (
const auto& peak : signalPDF.getPDFPeaks()) {
999 minT0 = std::min(minT0, peak.t0);
1000 for (
const auto& gaus : signalPDF.getTTS()->getTTS()) {
1001 double sig2 = peak.wid + gaus.sigma * gaus.sigma + wid0;
1002 double x = pow(hit.time - peak.t0 - gaus.position, 2) / sig2;
1003 if (x > 100)
continue;
1004 double wt = signalFract * peak.nph * gaus.fraction /
sqrt(2 * M_PI * sig2) * exp(-x / 2);
1006 m_pulls.push_back(
Pull(hit.pixelID, hit.time, peak.t0, gaus.position,
sqrt(sig2), peak.fic - M_PI, wt));
1013 m_pulls.push_back(
Pull(hit.pixelID, hit.time, minT0, 0, 0, 0, bg));
1015 if (sum == 0)
return;
1016 for (
size_t i = i0; i <
m_pulls.size(); i++)
m_pulls[i].wt /= sum;
Provides a type-safe way to pass members of the chargedStableSet set.
static const double speedOfLight
[cm/ns]
Fast photon propagation in quartz optics.
int getNyb() const
Returns signed number of reflections in y after mirror and before prism.
const std::vector< PhotonState > & getExtraStates() const
Returns extra states.
int getNxm() const
Returns signed number of reflections in x before mirror.
int getNxb() const
Returns signed number of reflections in x after mirror and before prism.
void propagate(const PhotonState &photon, bool averaging=false) const
Propagate photon to photo-detector plane.
int getNym() const
Returns signed number of reflections in y before mirror.
int getNxe() const
Returns signed number of reflections in x inside prism.
bool getPropagationStatus() const
Returns propagation status.
int getNye() const
Returns signed number of reflections in y inside prism.
double m_cosTotal
cosine of total reflection angle
void setSignalPDF_prism()
Sets signal PDF for track crossing prism.
void setSignalPDF()
Sets signal PDF.
double m_bkgRate
estimated background hit rate
double m_beta
particle hypothesis beta
const BackgroundPDF * getBackgroundPDF() const
Returns background PDF.
const InverseRaytracer::CerenkovAngle & cerenkovAngle(double dE=0)
Returns cosine and sine of cerenkov angle.
bool detectionPositionX(double xM, int Nxm, std::vector< double > &xDs, double &minLen)
Calculates unfolded detection position from known reflection position on the mirror and emission poin...
const YScanner * m_yScanner
PDF expander in y.
std::vector< SignalPDF > m_signalPDFs
parameterized signal PDF in pixels (index = pixelID - 1)
void setSignalPDF_reflected()
Sets signal PDF for reflected photons.
LogL getBackgroundLogL() const
Returns extended log likelihood for background hypothesis using default time window.
PrismSolution prismSolution(const PixelPositions::PixelData &pixel, unsigned k, int nx)
General solution of inverse raytracing in prism: iterative procedure calling basic solution.
LogL getLogL() const
Returns extended log likelihood (using the default time window)
double m_maxTime
time window upper edge
const InverseRaytracer * m_inverseRaytracer
inverse ray-tracer
double m_minTime
time window lower edge
double getExpectedDeltaPhotons() const
Returns the expected number of delta-ray photons within the default time window.
EStoreOption m_storeOption
signal PDF storing option
double m_Fic
temporary storage for Cerenkov azimuthal angle
double derivativeOfReflectedX(double x, double xe, double ze, double zd) const
Returns the derivative of reflected position at given x.
bool setDerivatives(YScanner::Derivatives &D, double dL, double de, double dFic)
Sets the derivatives (numerically) using forward ray-tracing.
double m_groupIndex
group refractive index at mean photon energy
std::vector< const PDFConstructor * > m_pdfOtherTracks
most probable PDF's of other tracks in the module
double m_bkgPhotons
expected number of uniform background photons
const std::vector< LogL > & getPixelLogLs(double t0, double sigt=0) const
Returns extended log likelihoods in pixels for PDF shifted in time.
double propagationLosses(double E, double propLen, int nx, int ny, SignalPDF::EPeakType type) const
Returns photon propagation losses (bulk absorption, surface reflectivity, mirror reflectivity)
int getModuleID() const
Returns slot ID.
double m_deltaPhotons
expected number of delta-ray photons
const std::vector< Pull > & getPulls() const
Returns photon pulls w.r.t PDF peaks.
double getCosCerenkovAngle(double E) const
Returns cosine of Cerenkov angle at given photon energy.
double m_signalPhotons
expected number of signal photons
std::map< SignalPDF::EPeakType, int > m_ncallsExpandPDF
number of calls to expandSignalPDF
void appendPulls(const TOPTrack::SelectedHit &hit) const
Appends pulls of a photon hit.
double getExpectedPhotons() const
Returns the expected number of all photons within the default time window.
bool m_valid
cross-check flag, true if track is valid and all the pointers above are valid
DeltaRayPDF m_deltaRayPDF
delta-ray PDF
double deltaXD(double dFic, const InverseRaytracer::Solution &sol, double xD)
Returns the difference in xD between ray-traced solution rotated by dFic and input argument.
bool doRaytracingCorrections(const InverseRaytracer::Solution &sol, double dFic_dx, double xD)
Corrects the solution of inverse ray-tracing with fast ray-tracing.
bool raytrace(const FastRaytracer &rayTracer, double dL=0, double de=0, double dFic=0)
Forward ray-tracing (called by setDerivatives)
std::vector< Pull > m_pulls
photon pulls w.r.t PDF peaks
const FastRaytracer * m_fastRaytracer
fast ray-tracer
double findReflectionExtreme(double xE, double zE, double zD, int Nxm, double A, const RaytracerBase::Mirror &mirror) const
Finds the position on the mirror of the extreme reflection.
EStoreOption
Options for storing signal PDF parameters.
@ c_Reduced
only PDF peak data
const TOPTrack & m_track
temporary reference to track at TOP
std::vector< LogL > m_pixelLLs
pixel log likelihoods (index = pixelID - 1)
EPDFOption m_PDFOption
signal PDF construction option
const BackgroundPDF * m_backgroundPDF
background PDF
EPDFOption
Signal PDF construction options.
@ c_Optimal
y dependent only where necessary
@ c_Fine
y dependent everywhere
double pdfValue(int pixelID, double time, double timeErr, double sigt=0) const
Returns the value of PDF normalized to the number of expected photons.
void expandSignalPDF(unsigned col, const YScanner::Derivatives &D, SignalPDF::EPeakType type)
Expands signal PDF in y (y-scan)
std::vector< FastRaytracer > m_rayTracers
copies of fast ray-tracer used to compute derivatives
const Const::ChargedStable m_hypothesis
particle hypothesis
bool prismRaytrace(const PrismSolution &sol, double dL=0, double dFic=0, double de=0)
Do forward raytracing of inverse raytracing solution in prism.
double m_groupIndexDerivative
derivative (dn_g/dE) of group refractive index at mean photon energy
void initializePixelLogLs(double minTime, double maxTime) const
Initializes pixel log likelihoods.
bool rangeOfX(double z, double &xmi, double &xma)
Estimates range of unfolded x coordinate of the hits on given plane perpendicular to z-axis.
std::vector< TOPTrack::SelectedHit > m_selectedHits
selected photon hits
void setSignalPDF_direct()
Sets signal PDF for direct photons.
std::set< int > m_zeroPixels
collection of pixelID's with zero pdfValue
double m_f0
temporary value of signal PDF
double m_tof
time-of-flight from IP to average photon emission position
PDFConstructor(const TOPTrack &track, const Const::ChargedStable &hypothesis, EPDFOption PDFOption=c_Optimal, EStoreOption storeOption=c_Reduced, double overrideMass=0)
Class constructor.
State of the Cerenkov photon in the quartz optics.
Parametrization of signal PDF in a single pixel.
EPeakType
Enumerator for single PDF peak types.
@ c_Reflected
reflected photon
static TOPGeometryPar * Instance()
Static method to obtain the pointer to its instance.
double getAbsorptionLength(double energy) const
Returns bulk absorption lenght of quartz at given photon energy.
double getGroupIndex(double energy) const
Returns group refractive index of quartz at given photon energy.
Singleton class providing pre-constructed reconstruction objects.
static double getMaxTime()
Returns time window upper edge.
static double getMinTime()
Returns time window lower edge.
Reconstructed track at TOP.
Class to store variables with their name which were sent to the logging service.
double sqrt(double a)
sqrt for double
double unfold(double x, int nx, double A)
unfold a coordinate.
double clip(double x, int Nx, double A, double xmi, double xma)
Performs a clip on x w.r.t xmi and xma.
double within2PI(double angle)
Returns angle within 0 and 2PI.
Abstract base class for different kinds of events.
Solution of inverse ray-tracing.
Structure that enables defining a template function: direct photons.
Structure that enables defining a template function: reflected photons.
Useful data type for returning the results of log likelihood calculation.
double effectiveSignalYield
effective number of signal photons in data
unsigned numPhotons
detected number of photons
double logL
extended log likelihood
Solution of inverse raytracing in prism.
double L
emission position distance along particle trajectory
double cosFic
cosine of azimuthal Cerenkov angle
double sinFic
sine of azimuthal Cerenkov angle
double len
propagation length
Data type for storing photon pull w.r.t PDF peak.
position and size of a pixel
double yc
position of center in y
double xc
position of center in x
spherical mirror data in module local frame.
double xc
center of curvature in x
double zc
center of curvature in z
double t0
peak position [ns]
double fic
Cerenkov azimuthal angle.
double wid
peak width squared [ns^2]
double nph
normalized number of photons in a peak
selected photon hit from TOPDigits