PDFConstructor Class Reference

PDF construction and log likelihood determination for a given track and particle hypothesis. More...

#include <PDFConstructor.h>

Classes
struct	InverseRaytracerDirect
	Structure that enables defining a template function: direct photons. More...
struct	InverseRaytracerReflected
	Structure that enables defining a template function: reflected photons. More...
struct	LogL
	Useful data type for returning the results of log likelihood calculation. More...
struct	PrismSolution
	Solution of inverse raytracing in prism. More...
struct	Pull
	Data type for storing photon pull w.r.t PDF peak. More...

Public Types
enum	EPDFOption {   c_Rough = 0 ,   c_Fine = 1 ,   c_Optimal = 2 }
	Signal PDF construction options. More...
enum	EStoreOption {   c_Reduced = 0 ,   c_Full = 1 }
	Options for storing signal PDF parameters. More...

Public Member Functions
	PDFConstructor (const TOPTrack &track, const Const::ChargedStable &hypothesis, EPDFOption PDFOption=c_Optimal, EStoreOption storeOption=c_Reduced, double overrideMass=0)
	Class constructor.
bool	isValid () const
	Checks the object status.
void	switchOffDeltaRayPDF () const
	Exclude delta-ray PDF in log likelihood calculation.
void	switchOnDeltaRayPDF () const
	Include delta-ray PDF in log likelihood calculation (this is default)
void	switchDeltaRayPDF (bool deltaPDFOn) const
	Include or exclude delta-ray PDF in log likelihood calculation.
int	getModuleID () const
	Returns slot ID.
const Const::ChargedStable &	getHypothesis () const
	Returns particle hypothesis.
const std::vector< TOPTrack::SelectedHit > &	getSelectedHits () const
	Returns selected photon hits belonging to this slot.
double	getBkgRate () const
	Returns estimated background hit rate.
double	getCosTotal () const
	Returns cosine of total reflection angle.
double	getCosCerenkovAngle (double E) const
	Returns cosine of Cerenkov angle at given photon energy.
const std::vector< SignalPDF > &	getSignalPDF () const
	Returns signal PDF.
const BackgroundPDF *	getBackgroundPDF () const
	Returns background PDF.
const DeltaRayPDF &	getDeltaRayPDF () const
	Returns delta-ray PDF.
double	getExpectedSignalPhotons () const
	Returns the expected number of signal photons within the default time window.
double	getExpectedDeltaPhotons () const
	Returns the expected number of delta-ray photons within the default time window.
double	getExpectedBkgPhotons () const
	Returns the expected number of background photons within the default time window.
double	getExpectedPhotons () const
	Returns the expected number of all photons within the default time window.
double	getExpectedSignalPhotons (double minTime, double maxTime) const
	Returns the expected number of signal photons within a given time window.
double	getExpectedDeltaPhotons (double minTime, double maxTime) const
	Returns the expected number of delta-ray photons within a given time window.
double	getExpectedBkgPhotons (double minTime, double maxTime) const
	Returns the expected number of background photons within a given time window.
double	getExpectedPhotons (double minTime, double maxTime) const
	Returns the expected number of all photons within a given time window.
double	getPDFValue (int pixelID, double time, double timeErr, double sigt=0) const
	Returns PDF value.
LogL	getLogL () const
	Returns extended log likelihood (using the default time window)
LogL	getLogL (double t0, double sigt=0) const
	Returns extended log likelihood for PDF shifted in time.
LogL	getLogL (double t0, double minTime, double maxTime, double sigt=0) const
	Returns extended log likelihood for PDF shifted in time and using different time window.
LogL	getBackgroundLogL () const
	Returns extended log likelihood for background hypothesis using default time window.
LogL	getBackgroundLogL (double minTime, double maxTime) const
	Returns extended log likelihood for background hypothesis.
const std::vector< LogL > &	getPixelLogLs (double t0, double sigt=0) const
	Returns extended log likelihoods in pixels for PDF shifted in time.
const std::vector< LogL > &	getPixelLogLs (double t0, double minTime, double maxTime, double sigt=0) const
	Returns extended log likelihoods in pixels for PDF shifted in time and using diferent time window.
const std::vector< Pull > &	getPulls () const
	Returns photon pulls w.r.t PDF peaks.
int	getNCalls_setPDF (SignalPDF::EPeakType type) const
	Returns number of calls of template function setSignalPDF<T> for a given peak type.
int	getNCalls_expandPDF (SignalPDF::EPeakType type) const
	Returns number of calls of function expandSignalPDF for a given peak type.
const std::map< double, YScanner::Derivatives > &	getDerivatives () const
	Returns a collection of derivatives for debugging purposes.
void	appendPDFOther (const PDFConstructor *pdfOther)
	Append most probable PDF of other track in this module.
void	clearPDFOther ()
	Clear the container of PDF's of other tracks.

Private Member Functions
template<class T >
void	setSignalPDF (T &t, unsigned col, double xD, double zD, int Nxm=0, double xmMin=0, double xmMax=0)
	Template function: sets signal PDF in a pixel column and for specific reflection in x.
const InverseRaytracer::CerenkovAngle &	cerenkovAngle (double dE=0)
	Returns cosine and sine of cerenkov angle.
void	setSignalPDF ()
	Sets signal PDF.
void	setSignalPDF_direct ()
	Sets signal PDF for direct photons.
void	setSignalPDF_reflected ()
	Sets signal PDF for reflected photons.
void	setSignalPDF_prism ()
	Sets signal PDF for track crossing prism.
void	setSignalPDF_reflected (int Nxm, double xmMin, double xmMax)
	Sets signal PDF for reflected photons at given reflection number.
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 point.
bool	doRaytracingCorrections (const InverseRaytracer::Solution &sol, double dFic_dx, double xD)
	Corrects the solution of inverse ray-tracing with fast ray-tracing.
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	setDerivatives (YScanner::Derivatives &D, double dL, double de, double dFic)
	Sets the derivatives (numerically) using forward ray-tracing.
bool	raytrace (const FastRaytracer &rayTracer, double dL=0, double de=0, double dFic=0)
	Forward ray-tracing (called by setDerivatives)
double	propagationLosses (double E, double propLen, int nx, int ny, SignalPDF::EPeakType type) const
	Returns photon propagation losses (bulk absorption, surface reflectivity, mirror reflectivity)
void	expandSignalPDF (unsigned col, const YScanner::Derivatives &D, SignalPDF::EPeakType type)
	Expands signal PDF in y (y-scan)
bool	rangeOfX (double z, double &xmi, double &xma)
	Estimates range of unfolded x coordinate of the hits on given plane perpendicular to z-axis.
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.
double	derivativeOfReflectedX (double x, double xe, double ze, double zd) const
	Returns the derivative of reflected position at given x.
bool	prismRaytrace (const PrismSolution &sol, double dL=0, double dFic=0, double de=0)
	Do forward raytracing of inverse raytracing solution in prism.
PrismSolution	prismSolution (const PixelPositions::PixelData &pixel, unsigned k, int nx)
	General solution of inverse raytracing in prism: iterative procedure calling basic solution.
PrismSolution	prismSolution (const ROOT::Math::XYZPoint &rD, double L)
	Basic solution of inverse raytracing in prism: assuming straight line particle trajectory.
double	pdfValueSignal (int pixelID, double time, double timeErr, double sigt=0) const
	Returns the value of signal PDF normalized to the number of expected photons.
double	pdfValueSignalDelta (int pixelID, double time, double timeErr, double sigt=0) const
	Returns the value of signal + deltaRay PDF normalized to the number of expected photons.
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	initializePixelLogLs (double minTime, double maxTime) const
	Initializes pixel log likelihoods.
void	appendPulls (const TOPTrack::SelectedHit &hit) const
	Appends pulls of a photon hit.

Private Attributes
int	m_moduleID = 0
	slot ID
const TOPTrack &	m_track
	temporary reference to track at TOP
const Const::ChargedStable	m_hypothesis
	particle hypothesis
const InverseRaytracer *	m_inverseRaytracer = 0
	inverse ray-tracer
const FastRaytracer *	m_fastRaytracer = 0
	fast ray-tracer
std::vector< FastRaytracer >	m_rayTracers
	copies of fast ray-tracer used to compute derivatives
const YScanner *	m_yScanner = 0
	PDF expander in y.
const BackgroundPDF *	m_backgroundPDF = 0
	background PDF
DeltaRayPDF	m_deltaRayPDF
	delta-ray PDF
EPDFOption	m_PDFOption = c_Optimal
	signal PDF construction option
EStoreOption	m_storeOption = c_Reduced
	signal PDF storing option
bool	m_valid = false
	cross-check flag, true if track is valid and all the pointers above are valid
double	m_beta = 0
	particle hypothesis beta
double	m_tof = 0
	time-of-flight from IP to average photon emission position
double	m_groupIndex = 0
	group refractive index at mean photon energy
double	m_groupIndexDerivative = 0
	derivative (dn_g/dE) of group refractive index at mean photon energy
double	m_cosTotal = 0
	cosine of total reflection angle
double	m_minTime = 0
	time window lower edge
double	m_maxTime = 0
	time window upper edge
std::vector< TOPTrack::SelectedHit >	m_selectedHits
	selected photon hits
double	m_bkgRate = 0
	estimated background hit rate
std::vector< SignalPDF >	m_signalPDFs
	parameterized signal PDF in pixels (index = pixelID - 1)
double	m_signalPhotons = 0
	expected number of signal photons
double	m_bkgPhotons = 0
	expected number of uniform background photons
double	m_deltaPhotons = 0
	expected number of delta-ray photons
std::map< double, InverseRaytracer::CerenkovAngle >	m_cerenkovAngles
	sine and cosine of Cerenkov angles
double	m_dFic = 0
	temporary storage for dFic used in last call to deltaXD
double	m_Fic = 0
	temporary storage for Cerenkov azimuthal angle
std::map< SignalPDF::EPeakType, int >	m_ncallsSetPDF
	number of calls to setSignalPDF<T>
std::map< SignalPDF::EPeakType, int >	m_ncallsExpandPDF
	number of calls to expandSignalPDF
std::vector< LogL >	m_pixelLLs
	pixel log likelihoods (index = pixelID - 1)
std::vector< Pull >	m_pulls
	photon pulls w.r.t PDF peaks
bool	m_deltaPDFOn = true
	include/exclude delta-ray PDF in likelihood calculation
std::map< double, YScanner::Derivatives >	m_derivatives
	a map of xD and derivatives
std::set< int >	m_zeroPixels
	collection of pixelID's with zero pdfValue
std::vector< const PDFConstructor * >	m_pdfOtherTracks
	most probable PDF's of other tracks in the module
double	m_f0 = 0
	temporary value of signal PDF

Detailed Description

PDF construction and log likelihood determination for a given track and particle hypothesis.

Definition at line 34 of file PDFConstructor.h.

Member Enumeration Documentation

EPDFOption

enum EPDFOption

Signal PDF construction options.

Enumerator
c_Rough	no dependence on y
c_Fine	y dependent everywhere
c_Optimal	y dependent only where necessary

Definition at line 41 of file PDFConstructor.h.

                      {
        c_Rough = 0,   
        c_Fine = 1,    
        c_Optimal = 2  
      };

EStoreOption

enum EStoreOption

Options for storing signal PDF parameters.

Enumerator
c_Reduced	only PDF peak data
c_Full	also extra information and derivatives

Definition at line 50 of file PDFConstructor.h.

                        {
        c_Reduced = 0, 
        c_Full = 1     
      };

Constructor & Destructor Documentation

PDFConstructor()

PDFConstructor	(	const TOPTrack &	track,
		const Const::ChargedStable &	hypothesis,
		EPDFOption	PDFOption = c_Optimal,
		EStoreOption	storeOption = c_Reduced,
		double	overrideMass = 0
	)

Class constructor.

Parameters

track	track at TOP
hypothesis	particle hypothesis
PDFOption	signal PDF construction option
storeOption	signal PDF store option
overrideMass	alternative mass value to be used intestead of the one from hypothesis. Ignored if <= 0.

Definition at line 27 of file PDFConstructor.cc.

                                                                                                       :
      m_moduleID(track.getModuleID()), m_track(track), m_hypothesis(hypothesis),
      m_inverseRaytracer(TOPRecoManager::getInverseRaytracer(m_moduleID)),
      m_fastRaytracer(TOPRecoManager::getFastRaytracer(m_moduleID)),
      m_yScanner(TOPRecoManager::getYScanner(m_moduleID)),
      m_backgroundPDF(TOPRecoManager::getBackgroundPDF(m_moduleID)),
      m_deltaRayPDF(m_moduleID),
      m_PDFOption(PDFOption), m_storeOption(storeOption)
    {
      if (not track.isValid()) {
        B2ERROR("TOP::PDFConstructor: TOPTrack is not valid, cannot continue");
        return;
      }
 
      m_valid = m_inverseRaytracer != 0 and m_fastRaytracer != 0 and m_yScanner != 0 and m_backgroundPDF != 0;
      if (not m_valid) {
        B2ERROR("TOP::PDFConstructor: missing reconstruction objects, cannot continue");
        return;
      }
 
      m_beta = track.getBeta(hypothesis, overrideMass);
      m_yScanner->prepare(track.getMomentumMag(), m_beta, track.getLengthInQuartz());
      m_tof = track.getTOF(hypothesis, 0, overrideMass);
 
      m_groupIndex = TOPGeometryPar::Instance()->getGroupIndex(m_yScanner->getMeanEnergy());
      m_groupIndexDerivative = TOPGeometryPar::Instance()->getGroupIndexDerivative(m_yScanner->getMeanEnergy());
      m_cosTotal = m_yScanner->getCosTotal();
      m_minTime = TOPRecoManager::getMinTime();
      m_maxTime = TOPRecoManager::getMaxTime();
      m_selectedHits = track.getSelectedHits();
      m_bkgRate = track.getBkgRate();
 
      // prepare the memory for storing signal PDF
 
      const auto& pixelPositions = m_yScanner->getPixelPositions();
      int numPixels = pixelPositions.getNumPixels();
      const auto* geo = TOPGeometryPar::Instance()->getGeometry();
      for (int pixelID = 1; pixelID <= numPixels; pixelID++) {
        auto pmtType = pixelPositions.get(pixelID).pmtType;
        const auto& tts = geo->getTTS(pmtType);
        m_signalPDFs.push_back(SignalPDF(pixelID, tts));
      }
 
      // construct PDF
 
      if (m_yScanner->isAboveThreshold()) {
        setSignalPDF();
      }
 
      m_deltaRayPDF.prepare(track, hypothesis);
      m_deltaPhotons = m_deltaRayPDF.getNumPhotons();
 
      m_bkgPhotons = std::max(m_bkgRate * (m_maxTime - m_minTime), 0.1);
 
      // release the memory not needed anymore
 
      m_rayTracers.clear();
    }

Member Function Documentation

appendPDFOther()

void appendPDFOther ( const PDFConstructor *  pdfOther )

inline

Append most probable PDF of other track in this module.

Parameters

pdfOther

PDF of other track

Definition at line 377 of file PDFConstructor.h.

      {
        if (not pdfOther) return;
        m_pdfOtherTracks.push_back(pdfOther);
      }

appendPulls()

void appendPulls ( const TOPTrack::SelectedHit &  hit ) const

private

Appends pulls of a photon hit.

Parameters

hit	photon hit

Definition at line 986 of file PDFConstructor.cc.

    {
      unsigned k = hit.pixelID - 1;
      if (k >= m_signalPDFs.size()) return;
      const auto& signalPDF = m_signalPDFs[k];
 
      double sfot = m_signalPhotons + m_deltaPhotons + m_bkgPhotons;
      double signalFract = m_signalPhotons / sfot;
      double wid0 = hit.timeErr * hit.timeErr;
      double minT0 = m_maxTime;
      double sum = 0;
      auto i0 = m_pulls.size();
      for (const auto& peak : signalPDF.getPDFPeaks()) {
        minT0 = std::min(minT0, peak.t0);
        for (const auto& gaus : signalPDF.getTTS()->getTTS()) {
          double sig2 = peak.wid + gaus.sigma * gaus.sigma + wid0; // sigma squared!
          double x = pow(hit.time - peak.t0 - gaus.position, 2) / sig2;
          if (x > 100) continue;
          double wt = signalFract * peak.nph * gaus.fraction / sqrt(2 * M_PI * sig2) * exp(-x / 2);
          sum += wt;
          m_pulls.push_back(Pull(hit.pixelID, hit.time, peak.t0, gaus.position, sqrt(sig2), peak.fic - M_PI, wt));
        }
      }
 
      double bg = (m_deltaPhotons * m_deltaRayPDF.getPDFValue(hit.pixelID, hit.time) +
                   m_bkgPhotons * m_backgroundPDF->getPDFValue(hit.pixelID)) / sfot;
      sum += bg;
      m_pulls.push_back(Pull(hit.pixelID, hit.time, minT0, 0, 0, 0, bg));
 
      if (sum == 0) return;
      for (size_t i = i0; i < m_pulls.size(); i++) m_pulls[i].wt /= sum;
    }

cerenkovAngle()

const InverseRaytracer::CerenkovAngle & cerenkovAngle ( double  dE = 0 )

inlineprivate

Returns cosine and sine of cerenkov angle.

Parameters

dE	energy difference to mean photon energy

Returns

cosine and sine of cerenkov angle

Definition at line 762 of file PDFConstructor.h.

    {
      auto& cer = m_cerenkovAngles[dE];
      if (cer.cosThc == 0 and cer.sinThc == 0) {
        double meanE = m_yScanner->getMeanEnergy();
        double cosThc = getCosCerenkovAngle(meanE + dE);
        cer = InverseRaytracer::CerenkovAngle(cosThc);
      }
      return cer;
    }

clearPDFOther()

void clearPDFOther ( )

inline

Clear the container of PDF's of other tracks.

Definition at line 386 of file PDFConstructor.h.

{m_pdfOtherTracks.clear();}

deltaXD()

double deltaXD ( double  dFic, const InverseRaytracer::Solution &  sol, double  xD  )

inlineprivate

Returns the difference in xD between ray-traced solution rotated by dFic and input argument.

Parameters

dFic	rotation angle around track direction (delta Cerenkov azimuthal angle)
sol	solution of inverse raytracing
xD	unfolded detection position in x

Returns

difference in xD

Definition at line 748 of file PDFConstructor.h.

    {
      m_dFic = dFic;
      m_fastRaytracer->propagate(m_inverseRaytracer->getReconstructedPhoton(sol, dFic), true);
      if (not m_fastRaytracer->getPropagationStatus()) return std::numeric_limits<double>::quiet_NaN();
      return m_fastRaytracer->getXD() - xD;
    }

derivativeOfReflectedX()

double derivativeOfReflectedX ( double  x, double  xe, double  ze, double  zd  ) const

private

Returns the derivative of reflected position at given x.

This function is used to find the position of reflection extreme. All arguments must be given in the mirror frame and in units of mirror radius.

Parameters

x	position in x on the mirror
xe	unfolded emission position in x
ze	emission position in z
zd	unfolded detection position in z

Returns

the derivative

Definition at line 576 of file PDFConstructor.cc.

    {
      double z = sqrt(1 - x * x);
      double kx = (x - xe);
      double kz = (z - ze);
      double s = 2 * (kx * x + kz * z);
      double qx = kx - s * x;
      double qz = kz - s * z;
 
      double der_z = -x / z;
      double der_s = 2 * (kx + der_z * kz);
      double der_qx = (1 - s) - der_s * x;
      double der_qz = (1 - s) * der_z - der_s * z;
 
      return 1 - der_z * qx / qz + (zd - z) * (der_qx * qz - der_qz * qx) / (qz * qz);
    }

detectionPositionX()

bool detectionPositionX ( double  xM, int  Nxm, std::vector< double > &  xDs, double &  minLen  )

private

Calculates unfolded detection position from known reflection position on the mirror and emission point.

Parameters

xM	reflection position x on the mirror
Nxm	reflection number before the mirror
xDs	unfolded detection positions (appended at each call) [in/out]
minLen	minimum of propagation lengths (updated at each call) [in/out]

Returns

true on success

Definition at line 216 of file PDFConstructor.cc.

    {
      m_inverseRaytracer->clear();
      int i0 = m_inverseRaytracer->solveForReflectionPoint(xM, Nxm, m_track.getEmissionPoint(), cerenkovAngle());
      if (i0 < 0) return false;
 
      bool ok = false;
      for (unsigned i = 0; i < 2; i++) {
        if (not m_inverseRaytracer->getStatus(i)) continue;
        const auto& solutions = m_inverseRaytracer->getSolutions(i);
        const auto& sol = solutions[i0];
        xDs.push_back(sol.xD);
        minLen = std::min(minLen, sol.len);
        ok = true;
      }
 
      return ok;
    }

doRaytracingCorrections()

bool doRaytracingCorrections ( const InverseRaytracer::Solution &  sol, double  dFic_dx, double  xD  )

private

Corrects the solution of inverse ray-tracing with fast ray-tracing.

Corrected solution is available in m_fastRaytracer.

Parameters

sol	solution of inverse raytracing
dFic_dx	derivative of Cerenkov azimuthal angle over photon detection coordinate x
xD	unfolded detection position in x

Returns

true on success

Definition at line 236 of file PDFConstructor.cc.

    {
      const double precision = 0.01; // [cm]
 
      double x1 = 0;                    // x is dFic
      double y1 = deltaXD(x1, sol, xD); // y is the difference in xD
      if (isnan(y1)) return false;
      if (std::abs(y1) < precision) return m_fastRaytracer->getTotalReflStatus(m_cosTotal);
      int n1 = m_fastRaytracer->getNxm();
 
      double step = -dFic_dx * y1;
      for (int i = 1; i < 20; i++) { // search for zero-crossing interval
        double x2 = step * i;
        double y2 = deltaXD(x2, sol, xD);
        if (isnan(y2)) return false;
        int n2 = m_fastRaytracer->getNxm();
        if (n2 != n1) { // x2 is passing the discontinuity caused by different reflection number
          double x3 = x2;
          x2 = x1;
          y2 = y1;
          for (int k = 0; k < 20; k++) { // move x2 to discontinuity using bisection
            double x = (x2 + x3) / 2;
            double y = deltaXD(x, sol, xD);
            if (isnan(y)) return false;
            int n = m_fastRaytracer->getNxm();
            if (n == n1) {
              x2 = x;
              y2 = y;
            } else {
              x3 = x;
            }
          }
          if (y2 * y1 > 0) return false; // solution does not exist
        }
        if (std::abs(y2) < precision) return m_fastRaytracer->getTotalReflStatus(m_cosTotal);
        if (y2 * y1 < 0) { // zero-crossing interval is identified
          for (int k = 0; k < 20; k++) { // find zero-crossing using bisection
            double x = (x1 + x2) / 2;
            double y = deltaXD(x, sol, xD);
            if (isnan(y)) return false;
            if (std::abs(y) < precision) return m_fastRaytracer->getTotalReflStatus(m_cosTotal);
            if (y * y1 < 0) {
              x2 = x;
            } else {
              x1 = x;
              y1 = y;
            }
          }
          return m_fastRaytracer->getTotalReflStatus(m_cosTotal);
        }
        x1 = x2;
        y1 = y2;
      }
 
      return false;
    }

expandSignalPDF()

void expandSignalPDF ( unsigned  col, const YScanner::Derivatives &  D, SignalPDF::EPeakType  type  )

private

Expands signal PDF in y (y-scan)

Parameters

col	pixel column (0-based)
D	derivatives
type	peak type, e.g. direct or reflected

Definition at line 391 of file PDFConstructor.cc.

    {
      m_ncallsExpandPDF[type]++;
 
      double Len = m_fastRaytracer->getPropagationLen();
      double speedOfLightQuartz = Const::speedOfLight / m_groupIndex; // average speed of light in quartz
 
      // difference of propagation times of true and fliped prism
      double dTime = m_fastRaytracer->getPropagationLenDelta() / speedOfLightQuartz;
 
      // derivatives: dt/de, dt/dx, dt/dL
      double dt_de = (D.dLen_de + Len * m_groupIndexDerivative / m_groupIndex) / speedOfLightQuartz;
      double dt_dx = D.dLen_dx / speedOfLightQuartz;
      double dt_dL = D.dLen_dL / speedOfLightQuartz + 1 / m_beta / Const::speedOfLight;
 
      // contribution of multiple scattering in quartz
      double sigmaScat = D.dLen_de * m_yScanner->getSigmaScattering() / speedOfLightQuartz;
 
      // contribution of quartz surface roughness at single reflection and effective number of reflections
      double sigmaAlpha = D.dLen_de * m_yScanner->getSigmaAlpha() / speedOfLightQuartz;
      int Ny_eff = 2 * std::abs(m_fastRaytracer->getNym()) + std::abs(m_fastRaytracer->getNyb());
 
      const auto& pixel = m_yScanner->getPixelPositions().get(col + 1);
      double L = m_track.getLengthInQuartz();
 
      // sigma squared: pixel size, parallax, propagation time difference, multiple scattering, surface roughness
      double wid0 = (pow(dt_dx * pixel.Dx, 2) + pow(dt_dL * L, 2) + pow(dTime, 2)) / 12 + pow(sigmaScat, 2) +
                    pow(sigmaAlpha, 2) * Ny_eff;
 
      // sigma squared: adding chromatic contribution
      double wid = wid0 + pow(dt_de * m_yScanner->getRMSEnergy(), 2);
 
      double yB = m_fastRaytracer->getYB();
      const auto& photonStates = m_fastRaytracer->getPhotonStates();
      const auto& atPrismEntrance = photonStates[photonStates.size() - 2];
      double dydz = atPrismEntrance.getKy() / atPrismEntrance.getKz();
      if (m_fastRaytracer->getNyb() % 2 != 0) dydz = -dydz;
 
      bool doScan = (m_PDFOption == c_Fine);
      if (m_PDFOption == c_Optimal) {
        double time = m_tof + Len / speedOfLightQuartz;
        doScan = m_track.isScanRequired(col, time, wid);
      }
 
      m_yScanner->expand(col, yB, dydz, D, Ny_eff, doScan);
 
      double numPhotons = m_yScanner->getNumPhotons() * std::abs(D.dFic_dx * pixel.Dx);
      int nx = m_fastRaytracer->getNx();
      int ny = m_fastRaytracer->getNy();
      for (const auto& result : m_yScanner->getResults()) {
        double RQE = m_yScanner->getPixelEfficiencies().get(result.pixelID);
        if (RQE == 0) continue;
        auto& signalPDF = m_signalPDFs[result.pixelID - 1];
        double dE = result.e0 - m_yScanner->getMeanEnergy();
        double propLen = Len + D.dLen_de * dE;
        double speedOfLight = Const::speedOfLight / TOPGeometryPar::Instance()->getGroupIndex(result.e0);
 
        SignalPDF::PDFPeak peak;
        peak.t0 = m_tof + propLen / speedOfLight;
        peak.wid = wid0 + dt_de * dt_de * result.sigsq;
        peak.nph = numPhotons * result.sum * RQE * propagationLosses(result.e0, propLen, nx, ny, type);
        peak.fic = func::within2PI(m_Fic + D.dFic_de * dE);
        signalPDF.append(peak);
 
        if (m_storeOption == c_Reduced) continue;
 
        SignalPDF::PDFExtra extra;
        extra.thc = acos(getCosCerenkovAngle(result.e0));
        extra.e = result.e0;
        extra.sige = result.sigsq;
        extra.Nxm = m_fastRaytracer->getNxm();
        extra.Nxb = m_fastRaytracer->getNxb();
        extra.Nxe = m_fastRaytracer->getNxe();
        extra.Nym = m_fastRaytracer->getNym();
        extra.Nyb = m_fastRaytracer->getNyb();
        extra.Nye = m_fastRaytracer->getNye();
        extra.xD = m_fastRaytracer->getXD();
        extra.yD = m_fastRaytracer->getYD();
        extra.zD = m_fastRaytracer->getZD();
        extra.yB = m_fastRaytracer->getYB();
        const auto& firstState =  photonStates.front();
        extra.kxE = firstState.getKx();
        extra.kyE = firstState.getKy();
        extra.kzE = firstState.getKz();
        const auto& lastState =  photonStates.back();
        extra.kxD = lastState.getKx();
        extra.kyD = lastState.getKy();
        extra.kzD = lastState.getKz();
        extra.type = type;
        signalPDF.append(extra);
      }
 
    }

findReflectionExtreme()

double findReflectionExtreme ( double  xE, double  zE, double  zD, int  Nxm, double  A, const RaytracerBase::Mirror &  mirror  ) const

private

Finds the position on the mirror of the extreme reflection.

Parameters

xE	true emission position in x
zE	emission position in z
zD	detection position in z
Nxm	signed number of reflections before mirror
A	width of the mirror segment
mirror	spherical mirror data

Returns

position of the extreme if exists or -A/2

Definition at line 594 of file PDFConstructor.cc.

    {
 
      if (Nxm % 2 == 0) {
        xE = func::unfold(xE, -Nxm, A);
      } else {
        xE = func::unfold(xE, Nxm, A);
      }
 
      double xe = (xE - mirror.xc) / mirror.R;
      double ze = (zE - mirror.zc) / mirror.R;
      double zd = (zD - mirror.zc) / mirror.R;
 
      double Ah = A / 2;
 
      double x1 = (-Ah - mirror.xc) / mirror.R;
      double y1 = derivativeOfReflectedX(x1, xe, ze, zd);
      if (y1 != y1 or std::abs(y1) == INFINITY) return -Ah;
 
      double x2 = (Ah - mirror.xc) / mirror.R;
      double y2 = derivativeOfReflectedX(x2, xe, ze, zd);
      if (y2 != y2 or std::abs(y2) == INFINITY) return -Ah;
 
      if (y1 * y2 > 0) return -Ah; // no minimum or maximum
 
      for (int i = 0; i < 50; i++) {
        double x = (x1 + x2) / 2;
        double y = derivativeOfReflectedX(x, xe, ze, zd);
        if (y != y or std::abs(y) == INFINITY) return -Ah;
        if (y * y1 < 0) {
          x2 = x;
        } else {
          x1 = x;
          y1 = y;
        }
      }
      double x = (x1 + x2) / 2;
 
      return x * mirror.R + mirror.xc;
    }

getBackgroundLogL() [1/2]

LogL getBackgroundLogL ( ) const

inline

Returns extended log likelihood for background hypothesis using default time window.

Returns

log likelihood

Definition at line 317 of file PDFConstructor.h.

{return getBackgroundLogL(m_minTime, m_maxTime);}

getBackgroundLogL() [2/2]

PDFConstructor::LogL getBackgroundLogL	(	double	minTime,
		double	maxTime
	)		const

Returns extended log likelihood for background hypothesis.

Parameters

minTime	time window lower edge
maxTime	time window upper edge

Returns

log likelihood

Definition at line 893 of file PDFConstructor.cc.

    {
      if (not m_valid) {
        B2ERROR("TOP::PDFConstructor::getBackgroundLogL(): object status is invalid - cannot provide log likelihood");
        return LogL(0);
      }
 
      double bkgPhotons = m_bkgPhotons * (maxTime - minTime) / (m_maxTime - m_minTime);
 
      LogL LL(bkgPhotons);
      for (const auto& hit : m_selectedHits) {
        if (hit.time < minTime or hit.time > maxTime) continue;
        double f = bkgPhotons * m_backgroundPDF->getPDFValue(hit.pixelID);
        if (f <= 0) {
          auto ret = m_zeroPixels.insert(hit.pixelID);
          if (ret.second) {
            B2ERROR("TOP::PDFConstructor::getBackgroundLogL(): PDF value is zero or negative"
                    << LogVar("slotID", m_moduleID)
                    << LogVar("pixelID", hit.pixelID) << LogVar("time", hit.time) << LogVar("PDFValue", f));
          }
          continue;
        }
        LL.logL += log(f);
        LL.numPhotons++;
      }
      return LL;
    }

getBackgroundPDF()

const BackgroundPDF * getBackgroundPDF ( ) const

inline

Returns background PDF.

Returns

background PDF

Definition at line 179 of file PDFConstructor.h.

{return m_backgroundPDF;}

getBkgRate()

double getBkgRate ( ) const

inline

Returns estimated background hit rate.

Returns

background hit rate per module

Definition at line 154 of file PDFConstructor.h.

{return m_bkgRate;}

getCosCerenkovAngle()

double getCosCerenkovAngle ( double  E ) const

inline

Returns cosine of Cerenkov angle at given photon energy.

Parameters

E	photon energy [eV]

Returns

cosine of Cerenkov angle

Definition at line 756 of file PDFConstructor.h.

    {
      double refind = TOPGeometryPar::Instance()->getPhaseIndex(E);
      return std::min(1 / m_beta / refind, 1.0);
    }

getCosTotal()

double getCosTotal ( ) const

inline

Returns cosine of total reflection angle.

Returns

cosine of total reflection angle at mean photon energy for beta = 1

Definition at line 160 of file PDFConstructor.h.

{return m_cosTotal;}

getDeltaRayPDF()

const DeltaRayPDF & getDeltaRayPDF ( ) const

inline

Returns delta-ray PDF.

Returns

delta-ray PDF

Definition at line 185 of file PDFConstructor.h.

{return m_deltaRayPDF;}

getDerivatives()

const std::map< double, YScanner::Derivatives > & getDerivatives ( ) const

inline

Returns a collection of derivatives for debugging purposes.

The derivatives are available only for EStoreOption::c_Full

Returns

map of xD and derivatives

Definition at line 371 of file PDFConstructor.h.

{return m_derivatives;}

getExpectedBkgPhotons() [1/2]

double getExpectedBkgPhotons ( ) const

inline

Returns the expected number of background photons within the default time window.

Photons from other tracks in the module are counted here as background.

Returns

expected number of background photons

Definition at line 204 of file PDFConstructor.h.

      {
        double photons = m_bkgPhotons;
        for (const auto* other : m_pdfOtherTracks) {
          photons += other->getExpectedSignalPhotons() + other->getExpectedDeltaPhotons();
        }
        return photons;
      }

getExpectedBkgPhotons() [2/2]

double getExpectedBkgPhotons ( double  minTime, double  maxTime  ) const

inline

Returns the expected number of background photons within a given time window.

Photons from other tracks in the module are counted here as background.

Parameters

minTime	time window lower edge
maxTime	time window upper edge

Returns

expected number of background photons

Definition at line 255 of file PDFConstructor.h.

      {
        double photons = (maxTime - minTime) / (m_maxTime - m_minTime) * m_bkgPhotons;
        for (const auto* other : m_pdfOtherTracks) {
          photons += other->getExpectedSignalPhotons(minTime, maxTime) + other->getExpectedDeltaPhotons(minTime, maxTime);
        }
        return photons;
      }

getExpectedDeltaPhotons() [1/2]

double getExpectedDeltaPhotons ( ) const

inline

Returns the expected number of delta-ray photons within the default time window.

Returns

expected number of delta-ray photons (0 if modeling turned off)

Definition at line 197 of file PDFConstructor.h.

{return m_deltaPDFOn ? m_deltaPhotons : 0;}

getExpectedDeltaPhotons() [2/2]

double getExpectedDeltaPhotons ( double  minTime, double  maxTime  ) const

inline

Returns the expected number of delta-ray photons within a given time window.

Parameters

minTime	time window lower edge
maxTime	time window upper edge

Returns

expected number of delta-ray photons (0 if modeling turned off)

Definition at line 243 of file PDFConstructor.h.

      {
        return m_deltaPDFOn ? m_deltaRayPDF.getIntegral(minTime, maxTime) * m_deltaPhotons : 0;
      }

getExpectedPhotons() [1/2]

double getExpectedPhotons ( ) const

inline

Returns the expected number of all photons within the default time window.

Returns

expected number of photons (signal + background + delta-ray)

Definition at line 217 of file PDFConstructor.h.

      {
        return getExpectedSignalPhotons() + getExpectedDeltaPhotons() + getExpectedBkgPhotons();
      }

getExpectedPhotons() [2/2]

double getExpectedPhotons ( double  minTime, double  maxTime  ) const

inline

Returns the expected number of all photons within a given time window.

Parameters

minTime	time window lower edge
maxTime	time window upper edge

Returns

expected number of photons (signal + background + delta-ray)

Definition at line 270 of file PDFConstructor.h.

      {
        return getExpectedSignalPhotons(minTime, maxTime) + getExpectedDeltaPhotons(minTime, maxTime) +
               getExpectedBkgPhotons(minTime, maxTime);
      }

getExpectedSignalPhotons() [1/2]

double getExpectedSignalPhotons ( ) const

inline

Returns the expected number of signal photons within the default time window.

Returns

expected number of signal photons

Definition at line 191 of file PDFConstructor.h.

{return m_signalPhotons;}

getExpectedSignalPhotons() [2/2]

double getExpectedSignalPhotons ( double  minTime, double  maxTime  ) const

inline

Returns the expected number of signal photons within a given time window.

Parameters

minTime	time window lower edge
maxTime	time window upper edge

Returns

expected number of signal photons

Definition at line 228 of file PDFConstructor.h.

      {
        double ps = 0;
        for (const auto& signalPDF : m_signalPDFs) {
          ps += signalPDF.getIntegral(minTime, maxTime);
        }
        return ps * m_signalPhotons;
      }

getHypothesis()

const Const::ChargedStable & getHypothesis ( ) const

inline

Returns particle hypothesis.

Returns

particle hypothesis

Definition at line 142 of file PDFConstructor.h.

{return m_hypothesis;}

getLogL() [1/3]

PDFConstructor::LogL getLogL

(

)

const

Returns extended log likelihood (using the default time window)

Returns

log likelihood

Definition at line 837 of file PDFConstructor.cc.

    {
      if (not m_valid) {
        B2ERROR("TOP::PDFConstructor::getLogL(): object status is invalid - cannot provide log likelihood");
        return LogL(0);
      }
 
      LogL LL(getExpectedPhotons());
      for (const auto& hit : m_selectedHits) {
        if (hit.time < m_minTime or hit.time > m_maxTime) continue;
        double f = pdfValue(hit.pixelID, hit.time, hit.timeErr);
        if (f <= 0) {
          auto ret = m_zeroPixels.insert(hit.pixelID);
          if (ret.second) {
            B2ERROR("TOP::PDFConstructor::getLogL(): PDF value is zero or negative"
                    << LogVar("slotID", m_moduleID)
                    << LogVar("pixelID", hit.pixelID) << LogVar("time", hit.time) << LogVar("PDFValue", f));
          }
          continue;
        }
        LL.logL += log(f);
        LL.numPhotons++;
        LL.effectiveSignalYield += m_f0 / f;
      }
      return LL;
    }

getLogL() [2/3]

PDFConstructor::LogL getLogL	(	double	t0,
		double	minTime,
		double	maxTime,
		double	sigt = 0
	)		const

Returns extended log likelihood for PDF shifted in time and using different time window.

Parameters

t0	time shift
minTime	time window lower edge
maxTime	time window upper edge
sigt	additional time smearing

Returns

log likelihood

Definition at line 865 of file PDFConstructor.cc.

    {
      if (not m_valid) {
        B2ERROR("TOP::PDFConstructor::getLogL(): object status is invalid - cannot provide log likelihood");
        return LogL(0);
      }
 
      LogL LL(getExpectedPhotons(minTime - t0, maxTime - t0));
      for (const auto& hit : m_selectedHits) {
        if (hit.time < minTime or hit.time > maxTime) continue;
        double f = pdfValue(hit.pixelID, hit.time - t0, hit.timeErr, sigt);
        if (f <= 0) {
          auto ret = m_zeroPixels.insert(hit.pixelID);
          if (ret.second) {
            B2ERROR("TOP::PDFConstructor::getLogL(): PDF value is zero or negative"
                    << LogVar("slotID", m_moduleID)
                    << LogVar("pixelID", hit.pixelID) << LogVar("time", hit.time) << LogVar("PDFValue", f));
          }
          continue;
        }
        LL.logL += log(f);
        LL.numPhotons++;
        LL.effectiveSignalYield += m_f0 / f;
      }
      return LL;
    }

getLogL() [3/3]

LogL getLogL ( double  t0, double  sigt = 0  ) const

inline

Returns extended log likelihood for PDF shifted in time.

Parameters

t0	time shift
sigt	additional time smearing

Returns

log likelihood

Definition at line 301 of file PDFConstructor.h.

{return getLogL(t0, m_minTime, m_maxTime, sigt);}

getModuleID()

int getModuleID ( ) const

inline

Returns slot ID.

Returns

slot ID

Definition at line 136 of file PDFConstructor.h.

{return m_moduleID;}

getNCalls_expandPDF()

int getNCalls_expandPDF ( SignalPDF::EPeakType  type ) const

inline

Returns number of calls of function expandSignalPDF for a given peak type.

Parameters

type	PDF peak type

Returns

number of calls

Definition at line 364 of file PDFConstructor.h.

{return m_ncallsExpandPDF[type];}

getNCalls_setPDF()

int getNCalls_setPDF ( SignalPDF::EPeakType  type ) const

inline

Returns number of calls of template function setSignalPDF<T> for a given peak type.

Parameters

type	PDF peak type

Returns

number of calls

Definition at line 357 of file PDFConstructor.h.

{return m_ncallsSetPDF[type];}

getPDFValue()

double getPDFValue ( int  pixelID, double  time, double  timeErr, double  sigt = 0  ) const

inline

Returns PDF value.

Parameters

pixelID	pixel ID
time	photon hit time
timeErr	uncertainty of hit time
sigt	additional time smearing

Returns

PDF value

Definition at line 284 of file PDFConstructor.h.

      {
        return pdfValue(pixelID, time, timeErr, sigt) / getExpectedPhotons();
      }

getPixelLogLs() [1/2]

const std::vector< PDFConstructor::LogL > & getPixelLogLs	(	double	t0,
		double	minTime,
		double	maxTime,
		double	sigt = 0
	)		const

Returns extended log likelihoods in pixels for PDF shifted in time and using diferent time window.

Parameters

t0	time shift
minTime	time window lower edge
maxTime	time window upper edge
sigt	additional time smearing

Returns

pixel log likelihoods (index = pixelID - 1)

Definition at line 923 of file PDFConstructor.cc.

    {
      if (not m_valid) {
        B2ERROR("TOP::PDFConstructor::getPixelLogLs(): object status is invalid - cannot provide log likelihoods");
        return m_pixelLLs;
      }
 
      initializePixelLogLs(minTime - t0, maxTime - t0);
 
      for (const auto& hit : m_selectedHits) {
        if (hit.time < minTime or hit.time > maxTime) continue;
        double f = pdfValue(hit.pixelID, hit.time - t0, hit.timeErr, sigt);
        if (f <= 0) {
          auto ret = m_zeroPixels.insert(hit.pixelID);
          if (ret.second) {
            B2ERROR("TOP::PDFConstructor::getPixelLogLs(): PDF value is zero or negative"
                    << LogVar("slotID", m_moduleID)
                    << LogVar("pixelID", hit.pixelID) << LogVar("time", hit.time) << LogVar("PDFValue", f));
          }
          continue;
        }
        unsigned k = hit.pixelID - 1;
        auto& LL = m_pixelLLs[k];
        LL.logL += log(f);
        LL.numPhotons++;
        LL.effectiveSignalYield += m_f0 / f;
      }
 
      return m_pixelLLs;
    }

getPixelLogLs() [2/2]

const std::vector< LogL > & getPixelLogLs ( double  t0, double  sigt = 0  ) const

inline

Returns extended log likelihoods in pixels for PDF shifted in time.

Parameters

t0	time shift
sigt	additional time smearing

Returns

pixel log likelihoods (index = pixelID - 1)

Definition at line 333 of file PDFConstructor.h.

      {return getPixelLogLs(t0, m_minTime, m_maxTime, sigt);}

getPulls()

const std::vector< PDFConstructor::Pull > & getPulls

(

)

const

Returns photon pulls w.r.t PDF peaks.

Returns

pulls

Definition at line 974 of file PDFConstructor.cc.

    {
      if (m_pulls.empty() and m_valid) {
        for (const auto& hit : m_selectedHits) {
          if (hit.time < m_minTime or hit.time > m_maxTime) continue;
          appendPulls(hit);
        }
      }
 
      return m_pulls;
    }

getSelectedHits()

const std::vector< TOPTrack::SelectedHit > & getSelectedHits ( ) const

inline

Returns selected photon hits belonging to this slot.

Returns

selected photon hits

Definition at line 148 of file PDFConstructor.h.

{return m_selectedHits;}

getSignalPDF()

const std::vector< SignalPDF > & getSignalPDF ( ) const

inline

Returns signal PDF.

Returns

signal PDF in pixels (index = pixelID - 1)

Definition at line 173 of file PDFConstructor.h.

{return m_signalPDFs;}

initializePixelLogLs()

void initializePixelLogLs ( double  minTime, double  maxTime  ) const

private

Initializes pixel log likelihoods.

Parameters

minTime	time window lower edge
maxTime	time window upper edge

Definition at line 954 of file PDFConstructor.cc.

    {
      m_pixelLLs.clear();
 
      double pb = (maxTime - minTime) / (m_maxTime - m_minTime);
      double bfot = pb * m_bkgPhotons + getExpectedDeltaPhotons(minTime, maxTime);
      for (const auto* other : m_pdfOtherTracks) bfot += other->getExpectedDeltaPhotons(minTime, maxTime);
 
      const auto& pixelPDF = m_backgroundPDF->getPDF();
      for (const auto& signalPDF : m_signalPDFs) {
        unsigned k = signalPDF.getPixelID() - 1;
        double phot = signalPDF.getIntegral(minTime, maxTime) * m_signalPhotons + bfot * pixelPDF[k];
        for (const auto* other : m_pdfOtherTracks) {
          const auto& otherPDFs = other->getSignalPDF();
          phot += otherPDFs[k].getIntegral(minTime, maxTime) * other->getExpectedSignalPhotons();
        }
        m_pixelLLs.push_back(LogL(phot));
      }
    }

isValid()

bool isValid ( ) const

inline

Checks the object status.

Returns

true if track is valid and all requited reconstruction objects are available

Definition at line 114 of file PDFConstructor.h.

{return m_valid;}

pdfValue()

double pdfValue ( int  pixelID, double  time, double  timeErr, double  sigt = 0  ) const

inlineprivate

Returns the value of PDF normalized to the number of expected photons.

Parameters

pixelID	pixel ID
time	photon hit time
timeErr	uncertainty of hit time
sigt	additional time smearing

Returns

PDF value

Definition at line 734 of file PDFConstructor.h.

    {
 
      if (not m_valid) return 0;
 
      double f = pdfValueSignal(pixelID, time, timeErr, sigt); // signal
      m_f0 = f;
      if (m_deltaPDFOn) f += m_deltaPhotons * m_deltaRayPDF.getPDFValue(pixelID, time); // delta electrons
      f += m_bkgPhotons * m_backgroundPDF->getPDFValue(pixelID); // uniform background
      for (const auto* other : m_pdfOtherTracks) f += other->pdfValueSignalDelta(pixelID, time, timeErr, sigt); // other tracks
 
      return f;
    }

pdfValueSignal()

double pdfValueSignal ( int  pixelID, double  time, double  timeErr, double  sigt = 0  ) const

inlineprivate

Returns the value of signal PDF normalized to the number of expected photons.

Parameters

pixelID	pixel ID
time	photon hit time
timeErr	uncertainty of hit time
sigt	additional time smearing

Returns

PDF value

Definition at line 718 of file PDFConstructor.h.

    {
      unsigned k = pixelID - 1;
      if (k < m_signalPDFs.size() and m_valid) {
        return  m_signalPhotons * m_signalPDFs[k].getPDFValue(time, timeErr, sigt);
      }
      return 0;
    }

pdfValueSignalDelta()

double pdfValueSignalDelta ( int  pixelID, double  time, double  timeErr, double  sigt = 0  ) const

inlineprivate

Returns the value of signal + deltaRay PDF normalized to the number of expected photons.

Parameters

pixelID	pixel ID
time	photon hit time
timeErr	uncertainty of hit time
sigt	additional time smearing

Returns

PDF value

Definition at line 727 of file PDFConstructor.h.

    {
      double f = pdfValueSignal(pixelID, time, timeErr, sigt);
      if (m_deltaPDFOn) f += m_deltaPhotons * m_deltaRayPDF.getPDFValue(pixelID, time);
      return f;
    }

prismRaytrace()

bool prismRaytrace ( const PrismSolution &  sol, double  dL = 0, double  dFic = 0, double  de = 0  )

private

Do forward raytracing of inverse raytracing solution in prism.

Parameters

sol	solution of inverse raytracing in prism
dL	step in length along trajectory for derivative calculation
dFic	step in Cherenkov azimuthal angle for derivative calculation
de	step in photon energy for derivative calculation

Returns

true on success

Definition at line 760 of file PDFConstructor.cc.

    {
      const auto& emi = m_track.getEmissionPoint(sol.L + dL);
      const auto& trk = emi.trackAngles;
      const auto& cer = cerenkovAngle(de);
 
      double cosDFic = 1;
      double sinDFic = 0;
      if (dFic != 0) {
        cosDFic = cos(dFic);
        sinDFic = sin(dFic);
      }
      double cosFic = sol.cosFic * cosDFic - sol.sinFic * sinDFic;
      double sinFic = sol.sinFic * cosDFic + sol.cosFic * sinDFic;
      double a = trk.cosTh * cer.sinThc * cosFic + trk.sinTh * cer.cosThc;
      double b = cer.sinThc * sinFic;
      double kx = a * trk.cosFi - b * trk.sinFi;
      double ky = a * trk.sinFi + b * trk.cosFi;
      double kz = trk.cosTh * cer.cosThc - trk.sinTh * cer.sinThc * cosFic;
      PhotonState photon(emi.position, kx, ky, kz);
      m_fastRaytracer->propagate(photon);
 
      return m_fastRaytracer->getPropagationStatus();
    }

prismSolution() [1/2]

PDFConstructor::PrismSolution prismSolution ( const PixelPositions::PixelData &  pixel, unsigned  k, int  nx  )

private

General solution of inverse raytracing in prism: iterative procedure calling basic solution.

Parameters

pixel	pixel data
k	index of unfolded prism exit window
nx	number of reflections in x

Returns

solution

Definition at line 786 of file PDFConstructor.cc.

    {
      const auto& prism = m_inverseRaytracer->getPrism();
      const auto& win = prism.unfoldedWindows[k];
      double dz = std::abs(prism.zD - prism.zFlat);
      ROOT::Math::XYZPoint rD(func::unfold(pixel.xc, nx, prism.A),
                              pixel.yc * win.sy + win.y0 + win.ny * dz,
                              pixel.yc * win.sz + win.z0 + win.nz * dz);
 
      double L = 0;
      for (int iter = 0; iter < 100; iter++) {
        auto sol = prismSolution(rD, L);
        if (std::abs(sol.L) > m_track.getLengthInQuartz() / 2) return sol;
        if (std::abs(sol.L - L) < 0.01) return sol;
        L = sol.L;
      }
      B2DEBUG(20, "TOP::PDFConstructor::prismSolution: iterations not converging");
      return PrismSolution();
    }

prismSolution() [2/2]

PDFConstructor::PrismSolution prismSolution ( const ROOT::Math::XYZPoint &  rD, double  L  )

private

Basic solution of inverse raytracing in prism: assuming straight line particle trajectory.

Parameters

rD	unfolded pixel position
L	emission position distance along particle trajectory

Returns

solution

Definition at line 808 of file PDFConstructor.cc.

    {
      const auto& emi = m_track.getEmissionPoint(L);
 
      // transformation of detection position to system of particle
 
      auto r = rD - emi.position;
      const auto& trk = emi.trackAngles;
      double xx = r.X() * trk.cosFi + r.Y() * trk.sinFi;
      double y = -r.X() * trk.sinFi + r.Y() * trk.cosFi;
      double x = xx * trk.cosTh - r.Z() * trk.sinTh;
      double z = xx * trk.sinTh + r.Z() * trk.cosTh;
 
      // solution
 
      double rho = sqrt(x * x + y * y);
      const auto& cer = cerenkovAngle();
 
      PrismSolution sol;
      sol.len = rho / cer.sinThc;
      sol.L = L + z - sol.len * cer.cosThc;
      sol.cosFic = x / rho;
      sol.sinFic = y / rho;
 
      return sol;
    }

propagationLosses()

double propagationLosses ( double  E, double  propLen, int  nx, int  ny, SignalPDF::EPeakType  type  ) const

private

Returns photon propagation losses (bulk absorption, surface reflectivity, mirror reflectivity)

Parameters

E	photon energy
propLen	propagation length
nx	total number of reflections in x
ny	total number of reflections in y
type	peak type, e.g. direct or reflected

Returns

survival probability due to propagation losses

Definition at line 485 of file PDFConstructor.cc.

    {
      double bulk = TOPGeometryPar::Instance()->getAbsorptionLength(E);
      double surf = m_yScanner->getBars().front().reflectivity;
      double p = exp(-propLen / bulk) * pow(surf, std::abs(nx) + std::abs(ny));
      if (type == SignalPDF::c_Reflected) p *= std::min(m_yScanner->getMirror().reflectivity, 1.0);
      return p;
    }

rangeOfX()

bool rangeOfX ( double  z, double &  xmi, double &  xma  )

private

Estimates range of unfolded x coordinate of the hits on given plane perpendicular to z-axis.

Parameters

z	position of the plane
xmi	lower limit [out]
xma	upper limit [out]

Returns

true if at least some photons can reach the plane within the time window

Definition at line 495 of file PDFConstructor.cc.

    {
      double maxLen = (m_maxTime - m_tof) / m_groupIndex * Const::speedOfLight; // maximal propagation length
      if (maxLen < 0) return false;
 
      const auto& emission = m_track.getEmissionPoint();
      const auto& trk = emission.trackAngles;
      const auto& cer = cerenkovAngle();
 
      // range in x from propagation lenght limit
 
      double dz = z - emission.position.Z();
      double cosFicLimit = (trk.cosTh * cer.cosThc - dz / maxLen) / (trk.sinTh * cer.sinThc); // at maxLen
      double cosLimit = (dz > 0) ? cosFicLimit : -cosFicLimit;
      if (cosLimit < -1) return false; // photons cannot reach the plane at z within propagation lenght limit
 
      std::vector<double> xmima;
      double x0 = emission.position.X();
      if (cosLimit > 1) {
        xmima.push_back(x0 - maxLen);
        xmima.push_back(x0 + maxLen);
      } else {
        double a = trk.cosTh * cer.sinThc * cosFicLimit + trk.sinTh * cer.cosThc;
        double b = cer.sinThc * sqrt(1 - cosFicLimit * cosFicLimit);
        xmima.push_back(x0 + maxLen * (a * trk.cosFi - b * trk.sinFi));
        xmima.push_back(x0 + maxLen * (a * trk.cosFi + b * trk.sinFi));
        std::sort(xmima.begin(), xmima.end());
      }
      xmi = xmima[0];
      xma = xmima[1];
 
      // range in x from minimal/maximal possible extensions in x, if they exist (d(kx/kz)/dFic = 0)
 
      double theta = acos(trk.cosTh);
      if (dz < 0) theta = M_PI - theta; // rotation around x by 180 deg. (z -> -z, phi -> -phi)
      dz = std::abs(dz);
      double thetaCer = acos(cer.cosThc);
      if (theta - thetaCer >= M_PI / 2) return false; // photons cannot reach the plane at z
 
      std::vector<double> dxdz;
      double a = -cos(theta + thetaCer) * cos(theta - thetaCer);
      double b = sin(2 * theta) * trk.cosFi;
      double c = pow(trk.sinFi * cer.sinThc, 2) - pow(trk.cosFi, 2) * sin(theta + thetaCer) * sin(theta - thetaCer);
      double D = b * b - 4 * a * c;
      if (D < 0) return true; // minimum and maximum do not exist, range is given by propagation length limit
      if (a != 0) {
        D = sqrt(D);
        dxdz.push_back((-b - D) / 2 / a);
        dxdz.push_back((-b + D) / 2 / a);
      } else {
        if (b == 0) return true; // minimum and maximum do not exist, range is given by propagation length limit
        dxdz.push_back(-c / b);
        dxdz.push_back(copysign(INFINITY, b));
      }
      std::vector<double> cosFic(2, cosLimit);
      for (int i = 0; i < 2; i++) {
        if (std::abs(dxdz[i]) < INFINITY) {
          double aa = (dxdz[i] * cos(theta) - trk.cosFi * sin(theta)) * cer.cosThc;
          double bb = (dxdz[i] * sin(theta) + trk.cosFi * cos(theta)) * cer.sinThc;
          double dd = trk.sinFi * cer.sinThc;
          cosFic[i] = aa * bb / (bb * bb + dd * dd);
          double kz = cos(theta) * cer.cosThc - sin(theta) * cer.sinThc * cosFic[i];
          if (kz < 0) dxdz[i] = copysign(INFINITY, dxdz[1 - i] - dxdz[i]);
        }
      }
      if (dxdz[0] > dxdz[1]) {
        std::reverse(dxdz.begin(), dxdz.end());
        std::reverse(cosFic.begin(), cosFic.end());
      }
      for (int i = 0; i < 2; i++) {
        if (cosFic[i] < cosLimit) xmima[i] = x0 + dxdz[i] * dz;
      }
 
      // just to make sure xmi/xma are within the limits given by maximal propagation length
      xmi = std::max(xmima[0], x0 - maxLen);
      xma = std::min(xmima[1], x0 + maxLen);
 
      return xma > xmi;
    }

raytrace()

bool raytrace ( const FastRaytracer &  rayTracer, double  dL = 0, double  de = 0, double  dFic = 0  )

private

Forward ray-tracing (called by setDerivatives)

Parameters

rayTracer	forward ray-tracer object
dL	value to be added to trajectory running parameter [cm]
de	value to be added to photon energy [eV]
dFic	value to be added to Cherenkov azimuthal angle

Returns

true on success

Definition at line 357 of file PDFConstructor.cc.

    {
      const auto& emi = m_track.getEmissionPoint(dL);
      const auto& trk = emi.trackAngles;
      const auto& cer = cerenkovAngle(de);
 
      double fic = m_Fic + dFic;
      double cosFic = cos(fic);
      double sinFic = sin(fic);
      double a = trk.cosTh * cer.sinThc * cosFic + trk.sinTh * cer.cosThc;
      double b = cer.sinThc * sinFic;
      double kx = a * trk.cosFi - b * trk.sinFi;
      double ky = a * trk.sinFi + b * trk.cosFi;
      double kz = trk.cosTh * cer.cosThc - trk.sinTh * cer.sinThc * cosFic;
      PhotonState photon(emi.position, kx, ky, kz);
      rayTracer.propagate(photon, true);
      if (not rayTracer.getPropagationStatus()) return false;
 
      if (rayTracer.getNxm() != m_fastRaytracer->getNxm()) return false;
      if (rayTracer.getNym() != m_fastRaytracer->getNym()) return false;
      if (rayTracer.getNxe() != m_fastRaytracer->getNxe()) return false;
      if (rayTracer.getNye() != m_fastRaytracer->getNye()) return false;
      if (rayTracer.getNxb() != m_fastRaytracer->getNxb()) return false;
      if (rayTracer.getNyb() != m_fastRaytracer->getNyb()) return false;
 
      if (rayTracer.getExtraStates().empty() or m_fastRaytracer->getExtraStates().empty()) return true;
 
      if (rayTracer.getExtraStates().back().getNx() != m_fastRaytracer->getExtraStates().back().getNx()) return false;
      if (rayTracer.getExtraStates().back().getNy() != m_fastRaytracer->getExtraStates().back().getNy()) return false;
 
      return true;
    }

setDerivatives()

bool setDerivatives ( YScanner::Derivatives &  D, double  dL, double  de, double  dFic  )

private

Sets the derivatives (numerically) using forward ray-tracing.

Parameters

D	derivatives to be set
dL	step in trajectory running parameter [cm]
de	step in photon energy [eV]
dFic	step in Cherenkov azimuthal angle

Returns

true on success

Definition at line 294 of file PDFConstructor.cc.

    {
      while (m_rayTracers.size() < 3) m_rayTracers.push_back(*m_fastRaytracer); // push back a copy of the object
 
      bool ok = true;
      auto& rayTracer_dL = m_rayTracers[0];
      for (int i = 0; i < 10; i++) {
        ok = raytrace(rayTracer_dL, dL);
        if (ok) break;
        dL = - dL / 2;
      }
      if (not ok) return false;
 
      auto& rayTracer_de = m_rayTracers[1];
      for (int i = 0; i < 10; i++) {
        ok = raytrace(rayTracer_de, 0, de);
        if (ok) break;
        de = - de / 2;
      }
      if (not ok) return false;
 
      auto& rayTracer_dFic = m_rayTracers[2];
      for (int i = 0; i < 10; i++) {
        ok = raytrace(rayTracer_dFic, 0, 0, dFic);
        if (ok) break;
        dFic = - dFic / 2;
      }
      if (not ok) return false;
 
      // partial derivatives on L, e and Fic
 
      double dLen_dL = (rayTracer_dL.getPropagationLen() - m_fastRaytracer->getPropagationLen()) / dL;
      double dLen_de = (rayTracer_de.getPropagationLen() - m_fastRaytracer->getPropagationLen()) / de;
      double dLen_dFic = (rayTracer_dFic.getPropagationLen() - m_fastRaytracer->getPropagationLen()) / dFic;
 
      double dyB_dL = (rayTracer_dL.getYB() - m_fastRaytracer->getYB()) / dL;
      double dyB_de = (rayTracer_de.getYB() - m_fastRaytracer->getYB()) / de;
      double dyB_dFic = (rayTracer_dFic.getYB() - m_fastRaytracer->getYB()) / dFic;
 
      double dx_dL = (rayTracer_dL.getXD() - m_fastRaytracer->getXD()) / dL;
      double dx_de = (rayTracer_de.getXD() - m_fastRaytracer->getXD()) / de;
      double dx_dFic = (rayTracer_dFic.getXD() - m_fastRaytracer->getXD()) / dFic;
 
      if (dx_dFic == 0) return false;
 
      // derivatives on L, e, and x. Derivatives on L and e have to be given at constant x.
 
      D.dLen_dx = dLen_dFic / dx_dFic;
      D.dLen_de = dLen_de - dLen_dFic * dx_de / dx_dFic;
      D.dLen_dL = dLen_dL - dLen_dFic * dx_dL / dx_dFic;
 
      D.dyB_dx = dyB_dFic / dx_dFic;
      D.dyB_de = dyB_de - dyB_dFic * dx_de / dx_dFic;
      D.dyB_dL = dyB_dL - dyB_dFic * dx_dL / dx_dFic;
 
      D.dFic_dx = 1 / dx_dFic;
      D.dFic_de = - dx_de / dx_dFic;
      D.dFic_dL = - dx_dL / dx_dFic;
 
      return true;
    }

setSignalPDF() [1/2]

void setSignalPDF ( )

private

Sets signal PDF.

Definition at line 88 of file PDFConstructor.cc.

    {
      // construct PDF analytically
 
      const auto& prism = m_inverseRaytracer->getPrism();
 
      if (m_track.getEmissionPoint().position.Z() > prism.zR) {
        setSignalPDF_direct();
        setSignalPDF_reflected();
      } else {
        setSignalPDF_prism();
      }
 
      // count expected number of signal photons
 
      for (const auto& signalPDF : m_signalPDFs) {
        m_signalPhotons += signalPDF.getSum();
      }
      if (m_signalPhotons == 0) return;
 
      // normalize PDF
 
      for (auto& signalPDF : m_signalPDFs) {
        signalPDF.normalize(m_signalPhotons);
      }
    }

setSignalPDF() [2/2]

void setSignalPDF ( T &  t, unsigned  col, double  xD, double  zD, int  Nxm = 0, double  xmMin = 0, double  xmMax = 0  )

private

Template function: sets signal PDF in a pixel column and for specific reflection in x.

Parameters

t	template parameter of type InverseRaytracerDirect or InverseRaytracerReflected
col	pixel column (0-based)
xD	unfolded detection position in x
zD	detection position in z
Nxm	signed number of reflections in x before mirror (dummy for direct PDF)
xmMin	lower limit of the reflection positions on the mirror (dummy for direct PDF)
xmMax	upper limit of the reflection positions on the mirror (dummy for direct PDF)

Definition at line 775 of file PDFConstructor.h.

    {
      m_ncallsSetPDF[t.type]++;
 
      m_inverseRaytracer->clear();
      t.inverseRaytracer = m_inverseRaytracer;
 
      // central solutions
 
      int i0 = t.solve(xD, zD, Nxm, xmMin, xmMax, m_track.getEmissionPoint(), cerenkovAngle());
      if (i0 < 0 or not m_inverseRaytracer->getStatus()) return;
      int n = 0;
      for (unsigned i = 0; i < 2; i++) {
        if (not m_inverseRaytracer->getStatus(i)) continue;
        const auto& solutions = m_inverseRaytracer->getSolutions(i);
        const auto& sol = solutions[i0];
        double time = m_tof + sol.len * m_groupIndex / Const::speedOfLight;
        if (time > m_maxTime + 1.0) continue;
        n++;
      }
      if (n == 0) return;
 
      // solutions with xD displaced by dx
 
      double dx = 0.1; // cm
      int i_dx = t.solve(xD + dx, zD, Nxm, xmMin, xmMax, m_track.getEmissionPoint(), cerenkovAngle(), dx);
      if (i_dx < 0) return;
      int k = 0;
      while (m_inverseRaytracer->isNymDifferent()) { // get rid of discontinuities
        if (k > 8) {
          B2DEBUG(20, "TOP::PDFConstructor::setSignalPDF: failed to find the same Nym (dx)");
          return;
        }
        dx = - dx / 2;
        i_dx = t.solve(xD + dx, zD, Nxm, xmMin, xmMax, m_track.getEmissionPoint(), cerenkovAngle(), dx);
        if (i_dx < 0) return;
        k++;
      }
 
      // loop over the two solutions, do ray-tracing corrections, compute the derivatives and expand PDF in y
 
      for (unsigned i = 0; i < 2; i++) {
        if (not m_inverseRaytracer->getStatus(i)) continue;
        const auto& solutions = m_inverseRaytracer->getSolutions(i);
        const auto& sol = solutions[i0];
        const auto& sol_dx = solutions[i_dx];
 
        // compute dFic/dx needed in raytracing corrections
 
        double dFic_dx = YScanner::Derivatives::dFic_d(sol, sol_dx);
 
        // do raytracing corrections
 
        m_dFic = 0;
        bool ok = doRaytracingCorrections(sol, dFic_dx, xD);
        if (not ok) continue;
        m_Fic = sol.getFic() + m_dFic;
 
        // check if time is still within the reconstruction range
 
        double time = m_tof + m_fastRaytracer->getPropagationLen() * m_groupIndex / Const::speedOfLight;
        if (time > m_maxTime) continue;
 
        // compute the derivatives
 
        YScanner::Derivatives D;
        ok = setDerivatives(D, 0.1, 0.1, 0.001);
        if (not ok) {
          B2DEBUG(20, "TOP::PDFConstructor::setSignalPDF: failed to determine derivatives: "
                  << LogVar("track momentum", m_track.getMomentumMag())
                  << LogVar("impact local z", m_track.getEmissionPoint().position.Z())
                  << LogVar("xD", xD)
                  << LogVar("Nxm", Nxm)
                  << LogVar("time", time)
                  << LogVar("peak type", t.type));
          continue;
        }
        if (m_storeOption == c_Full) m_derivatives[xD] = D;
 
        // expand PDF in y
 
        expandSignalPDF(col, D, t.type);
      }
    }

setSignalPDF_direct()

void setSignalPDF_direct ( )

private

Sets signal PDF for direct photons.

Definition at line 117 of file PDFConstructor.cc.

    {
      const auto& pixelPositions = m_yScanner->getPixelPositions();
      const auto& bar = m_inverseRaytracer->getBars().front();
      const auto& prism = m_inverseRaytracer->getPrism();
 
      // determine the range of number of reflections in x
 
      double xmi = 0, xma = 0;
      bool ok = rangeOfX(prism.zD, xmi, xma);
      if (not ok) return;
      int kmi = lround(xmi / bar.A);
      int kma = lround(xma / bar.A);
 
      // loop over reflections in x and over pixel columns
 
      for (int k = kmi; k <= kma; k++) {
        for (unsigned col = 0; col < pixelPositions.getNumPixelColumns(); col++) {
          const auto& pixel = pixelPositions.get(col + 1);
          if (pixel.Dx == 0) continue;
          double xD = func::unfold(pixel.xc, k, bar.A);
          if (xD < xmi or xD > xma) continue;
          InverseRaytracerDirect direct;
          setSignalPDF(direct, col, xD, prism.zD);
        }
      }
    }

setSignalPDF_prism()

void setSignalPDF_prism ( )

private

Sets signal PDF for track crossing prism.

Definition at line 638 of file PDFConstructor.cc.

    {
      const auto& pixelPositions = m_yScanner->getPixelPositions();
      const auto& prism = m_inverseRaytracer->getPrism();
      double speedOfLightQuartz = Const::speedOfLight / m_groupIndex; // average speed of light in quartz
 
      double xE = m_track.getEmissionPoint().position.X();
      int nxmi = (xE > 0) ? 0 : -1;
      int nxma = (xE > 0) ? 1 : 0;
 
      for (const auto& pixel : pixelPositions.getPixels()) {
        if (not m_yScanner->getPixelMasks().isActive(pixel.ID)) continue;
        double RQE = m_yScanner->getPixelEfficiencies().get(pixel.ID);
        if (RQE == 0) continue;
        auto& signalPDF = m_signalPDFs[pixel.ID - 1];
        for (int Nxe = nxmi; Nxe <= nxma; Nxe++) {
          for (size_t k = 0; k < prism.unfoldedWindows.size(); k++) {
            const auto sol = prismSolution(pixel, k, Nxe);
            if (sol.len == 0 or std::abs(sol.L) > m_track.getLengthInQuartz() / 2) continue;
 
            bool ok = prismRaytrace(sol);
            if (not ok) continue;
            int Nye = k - prism.k0;
            if (Nye != m_fastRaytracer->getNy() or Nxe != m_fastRaytracer->getNx()) continue;
            if (not m_fastRaytracer->getTotalReflStatus(m_cosTotal)) continue;
            const auto firstState = m_fastRaytracer->getPhotonStates().front();  // a copy of
            const auto lastState = m_fastRaytracer->getPhotonStates().back();  // a copy of
 
            double slope = lastState.getKy() / lastState.getKz();
            double dz = prism.zD - prism.zFlat;
            double y2 = std::min(pixel.yc + pixel.Dy / 2, prism.yUp + slope * dz);
            double y1 = std::max(pixel.yc - pixel.Dy / 2, prism.yDown + slope * dz);
            double Dy = y2 - y1;
            if (Dy < 0) continue;
 
            double dL = 0.1; // cm
            for (int i = 0; i < 4; i++) {
              ok = prismRaytrace(sol, dL);
              ok = ok and Nye == m_fastRaytracer->getNy() and Nxe == m_fastRaytracer->getNx();
              if (ok) break;
              dL = - dL / 2;
            }
            if (not ok) continue;
            const auto lastState_dL = m_fastRaytracer->getPhotonStates().back();  // a copy of
 
            double dFic = 0.01; // rad
            for (int i = 0; i < 4; i++) {
              ok = prismRaytrace(sol, 0, dFic);
              ok = ok and Nye == m_fastRaytracer->getNy() and Nxe == m_fastRaytracer->getNx();
              if (ok) break;
              dFic = - dFic / 2;
            }
            if (not ok) continue;
            const auto lastState_dFic = m_fastRaytracer->getPhotonStates().back();  // a copy of
 
            double de = 0.1; // eV
            for (int i = 0; i < 4; i++) {
              ok = prismRaytrace(sol, 0, 0, de);
              ok = ok and Nye == m_fastRaytracer->getNy() and Nxe == m_fastRaytracer->getNx();
              if (ok) break;
              de = - de / 2;
            }
            if (not ok) continue;
            const auto lastState_de = m_fastRaytracer->getPhotonStates().back();  // a copy of
 
            double dx_dL = (lastState_dL.getX() - lastState.getX()) / dL;
            double dy_dL = (lastState_dL.getY() - lastState.getY()) / dL;
            double dx_dFic = (lastState_dFic.getX() - lastState.getX()) / dFic;
            double dy_dFic = (lastState_dFic.getY() - lastState.getY()) / dFic;
            double Jacobi = dx_dL * dy_dFic - dy_dL * dx_dFic;
            double numPhotons = m_yScanner->getNumPhotonsPerLen() * pixel.Dx * Dy / std::abs(Jacobi) * RQE;
 
            double dLen_de = (lastState_de.getPropagationLen() - lastState.getPropagationLen()) / de;
            double dLen_dL = (lastState_dL.getPropagationLen() - lastState.getPropagationLen()) / dL;
            double dLen_dFic = (lastState_dFic.getPropagationLen() - lastState.getPropagationLen()) / dFic;
 
            double dt_de = (dLen_de + sol.len * m_groupIndexDerivative / m_groupIndex) / speedOfLightQuartz;
            double dt_dL = dLen_dL / speedOfLightQuartz;
            double dt_dFic = dLen_dFic / speedOfLightQuartz;
 
            double chromatic = pow(dt_de * m_yScanner->getRMSEnergy(), 2);
 
            double DL = (dy_dFic * pixel.Dx - dx_dFic * pixel.Dy) / Jacobi;
            double DFic = (dx_dL * pixel.Dy - dy_dL * pixel.Dx) / Jacobi;
            double paralax = (pow(dt_dL * DL, 2) + pow(dt_dFic * DFic, 2)) / 12;
 
            double scattering = pow(dLen_de * m_yScanner->getSigmaScattering() / speedOfLightQuartz, 2);
 
            SignalPDF::PDFPeak peak;
            peak.t0 = m_tof + sol.L / m_beta / Const::speedOfLight + sol.len / speedOfLightQuartz;
            peak.wid = chromatic + paralax + scattering;
            peak.nph = 1 - exp(-numPhotons); // because photons that pile-up are counted as one
            peak.fic = atan2(sol.sinFic, sol.cosFic);
            signalPDF.append(peak);
 
            if (m_storeOption == c_Reduced) continue;
 
            SignalPDF::PDFExtra extra;
            extra.thc = acos(getCosCerenkovAngle(m_yScanner->getMeanEnergy()));
            extra.e = m_yScanner->getMeanEnergy();
            extra.sige = m_yScanner->getRMSEnergy();
            extra.Nxe = Nxe;
            extra.Nye = Nye;
            extra.xD = lastState.getXD();
            extra.yD = lastState.getYD();
            extra.zD = lastState.getZD();
            extra.kxE = firstState.getKx();
            extra.kyE = firstState.getKy();
            extra.kzE = firstState.getKz();
            extra.kxD = lastState.getKx();
            extra.kyD = lastState.getKy();
            extra.kzD = lastState.getKz();
            extra.type = SignalPDF::c_Direct;
            signalPDF.append(extra);
 
          } // reflections in y (unfolded prism windows)
        } // reflections in x
      } // pixels
 
    }

setSignalPDF_reflected() [1/2]

void setSignalPDF_reflected ( )

private

Sets signal PDF for reflected photons.

Definition at line 145 of file PDFConstructor.cc.

    {
      const auto& bar = m_inverseRaytracer->getBars().back();
      const auto& prism = m_inverseRaytracer->getPrism();
      const auto& mirror = m_inverseRaytracer->getMirror();
      const auto& emiPoint = m_track.getEmissionPoint().position;
 
      // determine the range of number of reflections in x before mirror
 
      double xmi = 0, xma = 0;
      bool ok = rangeOfX(mirror.zb, xmi, xma);
      if (not ok) return;
      int kmi = lround(xmi / bar.A);
      int kma = lround(xma / bar.A);
 
      // loop over reflections in x before mirror
 
      double xE = emiPoint.X();
      double zE = emiPoint.Z();
      double Ah = bar.A / 2;
      for (int k = kmi; k <= kma; k++) {
        double x0 = findReflectionExtreme(xE, zE, prism.zD, k, bar.A, mirror);
        x0 = func::clip(x0, k, bar.A, xmi, xma);
        double xL = func::clip(-Ah, k, bar.A, xmi, xma);
        double xR = func::clip(Ah, k, bar.A, xmi, xma);
        if (x0 > xL) setSignalPDF_reflected(k, xL, x0);
        if (x0 < xR) setSignalPDF_reflected(k, x0, xR);
      }
    }

setSignalPDF_reflected() [2/2]

void setSignalPDF_reflected ( int  Nxm, double  xmMin, double  xmMax  )

private

Sets signal PDF for reflected photons at given reflection number.

Parameters

Nxm	reflection number in x before mirror
xmMin	lower limit of the reflection positions on the mirror
xmMax	upper limit of the reflection positions on the mirror

Definition at line 176 of file PDFConstructor.cc.

    {
      const auto& pixelPositions = m_yScanner->getPixelPositions();
      const auto& bar = m_inverseRaytracer->getBars().back();
      const auto& prism = m_inverseRaytracer->getPrism();
 
      // determine the range of number of reflections in x after mirror
 
      std::vector<double> xDs;
      double minLen = 1e10;
      if (not detectionPositionX(xmMin, Nxm, xDs, minLen)) return;
      if (not detectionPositionX(xmMax, Nxm, xDs, minLen)) return;
      if (xDs.size() < 2) return;
 
      double minTime = m_tof + minLen * m_groupIndex / Const::speedOfLight;
      if (minTime > m_maxTime) return;
 
      std::sort(xDs.begin(), xDs.end());
      double xmi = xDs.front();
      double xma = xDs.back();
 
      int kmi = lround(xmi / bar.A);
      int kma = lround(xma / bar.A);
 
      // loop over reflections in x after mirror and over pixel columns
 
      for (int k = kmi; k <= kma; k++) {
        for (unsigned col = 0; col < pixelPositions.getNumPixelColumns(); col++) {
          const auto& pixel = pixelPositions.get(col + 1);
          if (pixel.Dx == 0) continue;
          double xD = func::unfold(pixel.xc, k, bar.A);
          if (xD < xmi or xD > xma) continue;
          InverseRaytracerReflected reflected;
          setSignalPDF(reflected, col, xD, prism.zD, Nxm, xmMin, xmMax);
        }
      }
 
    }

switchDeltaRayPDF()

void switchDeltaRayPDF ( bool  deltaPDFOn ) const

inline

Include or exclude delta-ray PDF in log likelihood calculation.

Parameters

deltaPDFOn

true = include, false = exclude

Definition at line 130 of file PDFConstructor.h.

{m_deltaPDFOn = deltaPDFOn;}

switchOffDeltaRayPDF()

void switchOffDeltaRayPDF ( ) const

inline

Exclude delta-ray PDF in log likelihood calculation.

Definition at line 119 of file PDFConstructor.h.

{m_deltaPDFOn = false;}

switchOnDeltaRayPDF()

void switchOnDeltaRayPDF ( ) const

inline

Include delta-ray PDF in log likelihood calculation (this is default)

Definition at line 124 of file PDFConstructor.h.

{m_deltaPDFOn = true;}

Member Data Documentation

m_backgroundPDF

const BackgroundPDF* m_backgroundPDF = 0

private

background PDF

Definition at line 679 of file PDFConstructor.h.

m_beta

double m_beta = 0

private

particle hypothesis beta

Definition at line 685 of file PDFConstructor.h.

m_bkgPhotons

double m_bkgPhotons = 0

private

expected number of uniform background photons

Definition at line 697 of file PDFConstructor.h.

m_bkgRate

double m_bkgRate = 0

private

estimated background hit rate

Definition at line 693 of file PDFConstructor.h.

m_cerenkovAngles

std::map<double, InverseRaytracer::CerenkovAngle> m_cerenkovAngles

private

sine and cosine of Cerenkov angles

Definition at line 700 of file PDFConstructor.h.

m_cosTotal

double m_cosTotal = 0

private

cosine of total reflection angle

Definition at line 689 of file PDFConstructor.h.

m_deltaPDFOn

bool m_deltaPDFOn = true

mutableprivate

include/exclude delta-ray PDF in likelihood calculation

Definition at line 707 of file PDFConstructor.h.

m_deltaPhotons

double m_deltaPhotons = 0

private

expected number of delta-ray photons

Definition at line 698 of file PDFConstructor.h.

m_deltaRayPDF

DeltaRayPDF m_deltaRayPDF

private

delta-ray PDF

Definition at line 680 of file PDFConstructor.h.

m_derivatives

std::map<double, YScanner::Derivatives> m_derivatives

private

a map of xD and derivatives

Definition at line 708 of file PDFConstructor.h.

m_dFic

double m_dFic = 0

private

temporary storage for dFic used in last call to deltaXD

Definition at line 701 of file PDFConstructor.h.

m_f0

double m_f0 = 0

mutableprivate

temporary value of signal PDF

Definition at line 712 of file PDFConstructor.h.

m_fastRaytracer

const FastRaytracer* m_fastRaytracer = 0

private

fast ray-tracer

Definition at line 676 of file PDFConstructor.h.

m_Fic

double m_Fic = 0

private

temporary storage for Cerenkov azimuthal angle

Definition at line 702 of file PDFConstructor.h.

m_groupIndex

double m_groupIndex = 0

private

group refractive index at mean photon energy

Definition at line 687 of file PDFConstructor.h.

m_groupIndexDerivative

double m_groupIndexDerivative = 0

private

derivative (dn_g/dE) of group refractive index at mean photon energy

Definition at line 688 of file PDFConstructor.h.

m_hypothesis

const Const::ChargedStable m_hypothesis

private

particle hypothesis

Definition at line 674 of file PDFConstructor.h.

m_inverseRaytracer

const InverseRaytracer* m_inverseRaytracer = 0

private

inverse ray-tracer

Definition at line 675 of file PDFConstructor.h.

m_maxTime

double m_maxTime = 0

private

time window upper edge

Definition at line 691 of file PDFConstructor.h.

m_minTime

double m_minTime = 0

private

time window lower edge

Definition at line 690 of file PDFConstructor.h.

m_moduleID

int m_moduleID = 0

private

slot ID

Definition at line 672 of file PDFConstructor.h.

m_ncallsExpandPDF

std::map<SignalPDF::EPeakType, int> m_ncallsExpandPDF

mutableprivate

number of calls to expandSignalPDF

Definition at line 704 of file PDFConstructor.h.

m_ncallsSetPDF

std::map<SignalPDF::EPeakType, int> m_ncallsSetPDF

mutableprivate

number of calls to setSignalPDF<T>

Definition at line 703 of file PDFConstructor.h.

m_PDFOption

EPDFOption m_PDFOption = c_Optimal

private

signal PDF construction option

Definition at line 681 of file PDFConstructor.h.

m_pdfOtherTracks

std::vector<const PDFConstructor*> m_pdfOtherTracks

private

most probable PDF's of other tracks in the module

Definition at line 711 of file PDFConstructor.h.

m_pixelLLs

std::vector<LogL> m_pixelLLs

mutableprivate

pixel log likelihoods (index = pixelID - 1)

Definition at line 705 of file PDFConstructor.h.

m_pulls

std::vector<Pull> m_pulls

mutableprivate

photon pulls w.r.t PDF peaks

Definition at line 706 of file PDFConstructor.h.

m_rayTracers

std::vector<FastRaytracer> m_rayTracers

private

copies of fast ray-tracer used to compute derivatives

Definition at line 677 of file PDFConstructor.h.

m_selectedHits

std::vector<TOPTrack::SelectedHit> m_selectedHits

private

selected photon hits

Definition at line 692 of file PDFConstructor.h.

m_signalPDFs

std::vector<SignalPDF> m_signalPDFs

private

parameterized signal PDF in pixels (index = pixelID - 1)

Definition at line 695 of file PDFConstructor.h.

m_signalPhotons

double m_signalPhotons = 0

private

expected number of signal photons

Definition at line 696 of file PDFConstructor.h.

m_storeOption

EStoreOption m_storeOption = c_Reduced

private

signal PDF storing option

Definition at line 682 of file PDFConstructor.h.

m_tof

double m_tof = 0

private

time-of-flight from IP to average photon emission position

Definition at line 686 of file PDFConstructor.h.

m_track

const TOPTrack& m_track

private

temporary reference to track at TOP

Definition at line 673 of file PDFConstructor.h.

m_valid

bool m_valid = false

private

cross-check flag, true if track is valid and all the pointers above are valid

Definition at line 683 of file PDFConstructor.h.

m_yScanner

const YScanner* m_yScanner = 0

private

PDF expander in y.

Definition at line 678 of file PDFConstructor.h.

m_zeroPixels

std::set<int> m_zeroPixels

mutableprivate

collection of pixelID's with zero pdfValue

Definition at line 709 of file PDFConstructor.h.

---

The documentation for this class was generated from the following files:

• top/reconstruction_cpp/include/PDFConstructor.h
• top/reconstruction_cpp/src/PDFConstructor.cc