Utility for expanding the PDF in y direction. More...

#include <YScanner.h>

Inheritance diagram for YScanner:

Classes
struct	Derivatives
	Derivatives. More...

struct	PixelProjection
	Down-stream projection of a pixel to prism entrance window w/ a clip on bar exit thickness. More...

struct	Result
	Single PDF peak data. More...

struct	Table
	A table of equidistant entries. More...

struct	TableEntry
	Table entry. More...

struct	WindowProjection
	Geometry-only quantities needed to project a pixel to a given unfolded prism exit window. More...

Public Types
enum	EGeometry { c_Unified = 0 , c_Segmented = 1 }
	Treatment of quartz geometry. More...

enum	EOptics { c_SemiLinear = 0 , c_Exact = 1 }
	Treatment of spherical mirror optics. More...

Public Member Functions
	YScanner (int moduleID, unsigned N=64)
	Class constructor.

void	clear () const
	Clear mutable variables.

void	prepare (double momentum, double beta, double length) const
	Prepare for the PDF expansion in y for a given track mass hypothesis.

bool	isAboveThreshold () const
	Returns above Cerenkov threshold flag which is set in the prepare method.

void	expand (unsigned col, double yB, double dydz, const Derivatives &D, int Ny, bool doScan) const
	Performs the PDF expansion in y for a given pixel column using scan or merge methods.

const PixelPositions &	getPixelPositions () const
	Returns pixel positions and their sizes.

const PixelMasks &	getPixelMasks () const
	Returns pixel masks.

const PixelEfficiencies &	getPixelEfficiencies () const
	Returns pixel relative efficiencies.

const Table &	getEfficiencies () const
	Returns nominal photon detection efficiencies (PDE)

double	getCosTotal () const
	Returns cosine of total reflection angle.

double	getMomentum () const
	Returns particle momentum.

double	getBeta () const
	Returns particle beta.

double	getTrackLengthInQuartz () const
	Returns particle trajectory length inside quartz.

double	getNumPhotonsPerLen () const
	Returns number of photons per Cerenkov azimuthal angle per track length.

double	getNumPhotons () const
	Returns number of photons per Cerenkov azimuthal angle.

double	getMeanEnergy () const
	Returns mean photon energy.

double	getRMSEnergy () const
	Returns r.m.s of photon energy.

double	getMeanEnergyBeta1 () const
	Returns mean photon energy for beta = 1.

double	getRMSEnergyBeta1 () const
	Returns r.m.s of photon energy for beta = 1.

double	getSigmaScattering () const
	Returns r.m.s of multiple scattering angle in quartz converted to photon energy.

double	getSigmaAlpha () const
	Returns surface roughness parameter in units of photon energy.

double	getSurfaceReflectivity (unsigned n) const
	Returns the reflectivity of the bar surface raised to the given power.

const Table &	getEnergyDistribution () const
	Returns photon energy distribution.

const std::map< int, Table > &	getQuasyEnergyDistributions () const
	Returns photon energy distributions convoluted with multiple scattering and surface roughness.

const std::vector< Result > &	getResults () const
	Returns the results of PDF expansion in y.

bool	isScanDone () const
	Checks which expansion method was used.

int	getModuleID () const
	Returns slot ID.

EGeometry	getGeometry () const
	Returns quartz geometry treatment.

EOptics	getOptics () const
	Returns treatment of spherical mirror optics.

const std::vector< BarSegment > &	getBars () const
	Returns geometry data of bar segments.

const Mirror &	getMirror () const
	Returns geometry data of spherical mirror.

const Prism &	getPrism () const
	Returns geometry data of prism.

void	setMirrorCenter (double xc, double yc)
	Sets the mirror center-of-curvature.

Static Public Member Functions
static void	setScanLimits (int maxReflections)
	Sets parameters for selection between expand methods.

Protected Attributes
int	m_moduleID = 0
	slot ID

EGeometry	m_geometry = c_Unified
	quartz geometry

EOptics	m_optics = c_SemiLinear
	spherical mirror optics

std::vector< BarSegment >	m_bars
	geometry data of bar segments

Mirror	m_mirror
	spherical mirror geometry data

Prism	m_prism
	prism geometry data

Private Member Functions
double	setEnergyDistribution (double beta) const
	Sets photon energy distribution and mean photon energy according to nominal PDE and particle beta.

void	setQuasyEnergyDistribution (double sigma) const
	Sets photon energy distribution convoluted with a normalized Gaussian.

void	integrate (const EnergyMask *energyMask, double Ecp, Result &result) const
	Integrates quasy energy distribution multiplied with energy mask.

void	projectPixel (double yc, double size, int k, double dydz, PixelProjection proj[2]) const
	Calculates projections of a pixel to prism entrance window (going down-stream the photon).

double	tabulateSurfaceReflectivity (unsigned n) const
	Extends the table of surface reflectivity powers up to n and returns the value.

void	scan (unsigned col, double yB, double dydz, const Derivatives &D, int j1, int j2) const
	Performs expansion w/ the scan over reflections.

void	merge (unsigned col, double dydz, int j1, int j2) const
	Performs expansion by merging all reflections.

PixelMasks &	pixelMasks ()
	Returns non-const pixel masks.

PixelEfficiencies &	pixelEfficiencies ()
	Returns non-const pixel relative efficiencies.

Private Attributes
PixelPositions	m_pixelPositions
	positions and sizes of pixels

PixelMasks	m_pixelMasks
	pixel masks

PixelEfficiencies	m_pixelEfficiencies
	pixel relative efficiencies

Table	m_efficiency
	nominal photon detection efficiencies (PDE)

std::vector< WindowProjection >	m_windowProjections
	pixel projection constants of unfolded prism exit windows

std::vector< double >	m_surfaceReflectivities
	bar surface reflectivity to the power of the index; filled on demand, never cleared

double	m_prismZR = 0
	z of the prism-bar joint (copy of m_prism.zR)

double	m_halfBarThickness = 0
	half thickness of the bar at prism entrance

double	m_meanE0 = 0
	mean photon energy for beta = 1

double	m_rmsE0 = 0
	r.m.s of photon energy for beta = 1

double	m_cosTotal = 0
	cosine of total reflection angle

double	m_momentum = 0
	particle momentum magnitude

double	m_beta = 0
	particle beta

double	m_length = 0
	length of particle trajectory inside quartz

double	m_numPhotons = 0
	number of photons per Cerenkov azimuthal angle per track length

double	m_meanE = 0
	mean photon energy

double	m_rmsE = 0
	r.m.s of photon energy

double	m_sigmaScat = 0
	r.m.s.

double	m_sigmaAlpha = 0
	surface roughness parameter in photon energy units

Table	m_energyDistribution
	photon energy distribution

std::map< int, Table >	m_quasyEnergyDistributions
	photon energy distributions convoluted with Gaussian of different widths

Table *	m_quasyEnergyDistribution = nullptr
	a pointer to the element in m_quasyEnergyDistributions

bool	m_aboveThreshold = false
	true if beta is above the Cerenkov threshold

std::vector< Result >	m_results
	results of PDF expansion in y

bool	m_scanDone = false
	true if scan performed, false if reflections just merged

Static Private Attributes
static int	s_maxReflections = 16
	maximal number of reflections to perform scan

static unsigned	s_maxTabulatedPower = 1024
	maximal power of the surface reflectivity that is tabulated

Friends
class	TOPRecoManager

Detailed Description

Utility for expanding the PDF in y direction.

Definition at line 33 of file YScanner.h.

Member Enumeration Documentation

◆ EGeometry

enum EGeometry

inherited

Treatment of quartz geometry.

Enumerator
c_Unified	single bar with average width and thickness
c_Segmented	segmented bars

Definition at line 33 of file RaytracerBase.h.

                     {
        c_Unified = 0,   
        c_Segmented = 1  
      };

◆ EOptics

enum EOptics

inherited

Treatment of spherical mirror optics.

Enumerator
c_SemiLinear	semi-linear approximation
c_Exact	exact optics

Definition at line 41 of file RaytracerBase.h.

                   {
        c_SemiLinear = 0, 
        c_Exact = 1       
      };

Constructor & Destructor Documentation

◆ YScanner()

YScanner	(	int	moduleID,
		unsigned	N = 64 )

explicit

Class constructor.

Sets pixel positions of a given module and nominal photon detection efficiency.

Parameters

moduleID	slot ID
N	size of nominal photon detection efficiency table

Definition at line 48 of file YScanner.cc.

                                              : RaytracerBase(moduleID, c_Unified, c_SemiLinear),
      m_pixelPositions(PixelPositions(moduleID)),
      m_pixelMasks(PixelMasks(moduleID)),
      m_pixelEfficiencies(PixelEfficiencies(moduleID))
    {
      if (N < 2) {
        B2FATAL("TOP::YScanner: N must be > 1");
        return;
      }
 
      // set the pixel projection constants of the unfolded prism exit windows
 
      m_prismZR = m_prism.zR;
      m_halfBarThickness = m_bars.front().B / 2;
      double dz = std::abs(m_prism.zD - m_prism.zFlat);
      for (size_t k = 0; k < m_prism.unfoldedWindows.size(); k++) {
        const auto& win = m_prism.unfoldedWindows[k];
        WindowProjection projection;
        projection.sy = win.sy;
        projection.sz = win.sz;
        projection.y0 = win.y0 + win.ny * dz;
        projection.z0 = win.z0 + win.nz * dz;
        projection.evenReflection = ((static_cast<int>(k) - m_prism.k0) % 2 == 0);
        m_windowProjections.push_back(projection);
      }
 
      // set the table of nominal photon detection efficiencies (incl. wavelength filter)
 
      const auto* topgp = TOPGeometryPar::Instance();
      const auto* geo = topgp->getGeometry();
      auto qe = geo->getNominalQE(); // get a copy
      qe.applyFilterTransmission(geo->getWavelengthFilter());
 
      double minE = TOPGeometryPar::c_hc / qe.getMaxLambda();
      double maxE = TOPGeometryPar::c_hc / qe.getMinLambda();
      if (minE >= maxE) {
        B2FATAL("TOP::YScanner: quantum efficiency found zero for all wavelengths");
        return;
      }
      m_efficiency.set(minE, (maxE - minE) / (N - 1));
 
      const auto& tdc = geo->getNominalTDC();
      for (unsigned i = 0; i < N; i++) {
        double e = m_efficiency.getX(i);
        double lambda = TOPGeometryPar::c_hc / e;
        double effi = qe.getEfficiency(lambda) * tdc.getEfficiency();
        m_efficiency.entries.push_back(TableEntry(effi, e, e * e));
      }
 
      // set cosine of total reflection angle using photon mean energy for beta = 1
 
      double s = 0;
      double se = 0;
      double see = 0;
      for (const auto& entry : m_efficiency.entries) {
        double e = entry.x;
        double p = entry.y * (1 - 1 / pow(topgp->getPhaseIndex(e), 2));
        s += p;
        se += p * e;
        see += p * e * e;
      }
      if (s == 0) return;
      m_meanE0 = se / s;
      m_rmsE0 = sqrt(see / s - m_meanE0 * m_meanE0);
      m_cosTotal = sqrt(1 - 1 / pow(topgp->getPhaseIndex(m_meanE0), 2));
    }

Member Function Documentation

◆ clear()

void clear ( ) const

Clear mutable variables.

Definition at line 130 of file YScanner.cc.

    {
      m_momentum = 0;
      m_beta = 0;
      m_length = 0;
      m_numPhotons = 0;
      m_meanE = 0;
      m_rmsE = 0;
      m_sigmaScat = 0;
      m_sigmaAlpha = 0;
      m_energyDistribution.clear();
      m_quasyEnergyDistributions.clear();
      m_quasyEnergyDistribution = nullptr;
      m_aboveThreshold = false;
      m_results.clear();
      m_scanDone = false;
    }

◆ expand()

void expand	(	unsigned	col,
		double	yB,
		double	dydz,
		const Derivatives &	D,
		int	Ny,
		bool	doScan ) const

Performs the PDF expansion in y for a given pixel column using scan or merge methods.

Results accessible with getResults() method.

Parameters

col	pixel column number (0-based)
yB	unfolded coordinate y of photon at prism entrance (= Bar exit) plane
dydz	photon slope in y-z projection at prism entrance (dy/dz)
D	the derivatives
Ny	effective number of reflections
doScan	if true decide between scan and merge methods, if false always use merge method

Definition at line 274 of file YScanner.cc.

    {
      m_results.clear();
 
      if (D.dyB_de == 0) return;
 
      double sigma = sqrt(pow(m_sigmaScat, 2) + pow(m_sigmaAlpha, 2) * std::abs(Ny));
      setQuasyEnergyDistribution(sigma);
 
      double minE = m_quasyEnergyDistribution->getXmin();
      double maxE = m_quasyEnergyDistribution->getXmax();
      double pixDx  = m_pixelPositions.get(col + 1).Dx;
      double dely = (std::abs(D.dyB_dL) * m_length + std::abs(D.dyB_dx) * pixDx) / 2;
      double y1 = yB - dely;
      double y2 = yB + dely;
      if (D.dyB_de > 0) {
        y1 += D.dyB_de * (minE - m_meanE);
        y2 += D.dyB_de * (maxE - m_meanE);
      } else {
        y1 += D.dyB_de * (maxE - m_meanE);
        y2 += D.dyB_de * (minE - m_meanE);
      }
      double B = m_bars.front().B;
      int j1  = func::lround(y1 / B);
      int j2  = func::lround(y2 / B) + 1;
 
      if (doScan and j2 - j1 <= s_maxReflections) {
        scan(col, yB, dydz, D, j1, j2);
        m_scanDone = true;
      } else {
        merge(col, dydz, j1, j2);
        m_scanDone = false;
      }
    }

◆ getBars()

const std::vector< BarSegment > & getBars ( ) const

inlineinherited

Returns geometry data of bar segments.

Returns: geometry data of bar segments

Definition at line 161 of file RaytracerBase.h.

161{return m_bars;}

◆ getBeta()

double getBeta ( ) const

inline

Returns particle beta.

Returns: particle beta

Definition at line 318 of file YScanner.h.

318{return m_beta;}

◆ 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 306 of file YScanner.h.

306{return m_cosTotal;}

◆ getEfficiencies()

const Table & getEfficiencies ( ) const

inline

Returns nominal photon detection efficiencies (PDE)

Returns: nominal photon detection efficiencies

Definition at line 300 of file YScanner.h.

300{return m_efficiency;}

◆ getEnergyDistribution()

const Table & getEnergyDistribution ( ) const

inline

Returns photon energy distribution.

Returns: photon energy distribution

Definition at line 391 of file YScanner.h.

391{return m_energyDistribution;}

◆ getGeometry()

EGeometry getGeometry ( ) const

inlineinherited

Returns quartz geometry treatment.

Returns: quartz geometry treatment

Definition at line 149 of file RaytracerBase.h.

149{return m_geometry;}

◆ getMeanEnergy()

double getMeanEnergy ( ) const

inline

Returns mean photon energy.

Returns: mean photon energy

Definition at line 342 of file YScanner.h.

342{return m_meanE;}

◆ getMeanEnergyBeta1()

double getMeanEnergyBeta1 ( ) const

inline

Returns mean photon energy for beta = 1.

Returns: mean photon energy for beta = 1

Definition at line 354 of file YScanner.h.

354{return m_meanE0;}

◆ getMirror()

const Mirror & getMirror ( ) const

inlineinherited

Returns geometry data of spherical mirror.

Returns: geometry data of spherical mirror

Definition at line 167 of file RaytracerBase.h.

167{return m_mirror;}

◆ getModuleID()

int getModuleID ( ) const

inlineinherited

Returns slot ID.

Returns: slot ID

Definition at line 143 of file RaytracerBase.h.

143{return m_moduleID;}

◆ getMomentum()

double getMomentum ( ) const

inline

Returns particle momentum.

Returns: particle momentum

Definition at line 312 of file YScanner.h.

312{return m_momentum;}

◆ getNumPhotons()

double getNumPhotons ( ) const

inline

Returns number of photons per Cerenkov azimuthal angle.

Returns: number of photons per radian

Definition at line 336 of file YScanner.h.

336{return m_numPhotons * m_length;}

◆ getNumPhotonsPerLen()

double getNumPhotonsPerLen ( ) const

inline

Returns number of photons per Cerenkov azimuthal angle per track length.

Returns: number of photons per radian per centimeter

Definition at line 330 of file YScanner.h.

330{return m_numPhotons;}

◆ getOptics()

EOptics getOptics ( ) const

inlineinherited

Returns treatment of spherical mirror optics.

Returns: spherical mirror optics

Definition at line 155 of file RaytracerBase.h.

155{return m_optics;}

◆ getPixelEfficiencies()

const PixelEfficiencies & getPixelEfficiencies ( ) const

inline

Returns pixel relative efficiencies.

Returns: pixel relative efficiencies

Definition at line 294 of file YScanner.h.

294{return m_pixelEfficiencies;}

◆ getPixelMasks()

const PixelMasks & getPixelMasks ( ) const

inline

Returns pixel masks.

Returns: pixel masks

Definition at line 288 of file YScanner.h.

288{return m_pixelMasks;}

◆ getPixelPositions()

const PixelPositions & getPixelPositions ( ) const

inline

Returns pixel positions and their sizes.

Returns: pixel positions and their sizes in module local frame

Definition at line 282 of file YScanner.h.

282{return m_pixelPositions;}

◆ getPrism()

const Prism & getPrism ( ) const

inlineinherited

Returns geometry data of prism.

Returns: geometry data of prism

Definition at line 173 of file RaytracerBase.h.

173{return m_prism;}

◆ getQuasyEnergyDistributions()

const std::map< int, Table > & getQuasyEnergyDistributions ( ) const

inline

Returns photon energy distributions convoluted with multiple scattering and surface roughness.

Map entries correspond to different Gaussian widths due to different number of reflections.

Returns: photon energy distributions convoluted with multiple scattering and surface roughness

Definition at line 398 of file YScanner.h.

398{return m_quasyEnergyDistributions;}

◆ getResults()

const std::vector< Result > & getResults ( ) const

inline

Returns the results of PDF expansion in y.

Returns: results

Definition at line 404 of file YScanner.h.

404{return m_results;}

◆ getRMSEnergy()

double getRMSEnergy ( ) const

inline

Returns r.m.s of photon energy.

Returns: r.m.s of photon energy

Definition at line 348 of file YScanner.h.

348{return m_rmsE;}

◆ getRMSEnergyBeta1()

double getRMSEnergyBeta1 ( ) const

inline

Returns r.m.s of photon energy for beta = 1.

Returns: r.m.s of photon energy for beta = 1

Definition at line 360 of file YScanner.h.

360{return m_rmsE0;}

◆ getSigmaAlpha()

double getSigmaAlpha ( ) const

inline

Returns surface roughness parameter in units of photon energy.

Returns: surface roughness parameter [eV]

Definition at line 372 of file YScanner.h.

372{return m_sigmaAlpha;}

◆ getSigmaScattering()

double getSigmaScattering ( ) const

inline

Returns r.m.s of multiple scattering angle in quartz converted to photon energy.

Returns: r.m.s of multiple scattering angle in quartz [eV]

Definition at line 366 of file YScanner.h.

366{return m_sigmaScat;}

◆ getSurfaceReflectivity()

double getSurfaceReflectivity ( unsigned n ) const

inline

Returns the reflectivity of the bar surface raised to the given power.

The reflectivity is a constant of the module, therefore the powers are tabulated on first use and kept for the whole job (they do not depend on the track or the hypothesis).

Parameters

n	power, i.e. the number of reflections

Returns: reflectivity to the power of n

Definition at line 381 of file YScanner.h.

      {
        if (n < m_surfaceReflectivities.size()) return m_surfaceReflectivities[n];
        return tabulateSurfaceReflectivity(n);
      }

◆ getTrackLengthInQuartz()

double getTrackLengthInQuartz ( ) const

inline

Returns particle trajectory length inside quartz.

Returns: particle trajectory length inside quartz

Definition at line 324 of file YScanner.h.

324{return m_length;}

◆ integrate()

void integrate	(	const EnergyMask *	energyMask,
		double	Ecp,
		Result &	result ) const

private

Integrates quasy energy distribution multiplied with energy mask.

Parameters

energyMask	energy mask (pointer must be valid)
Ecp	position of the mask
result	updated result [in/out]

Definition at line 412 of file YScanner.cc.

    {
      const auto& mask = energyMask->getMask();
 
      if (mask.empty()) {
        // direct mask calculation
        for (size_t i = 0; i < m_quasyEnergyDistribution->entries.size(); i++) {
          double E = m_quasyEnergyDistribution->getX(i);
          double m = energyMask->getMask(E - Ecp);
          if (m > 0) {
            const auto& entry = m_quasyEnergyDistribution->entries[i];
            double s = entry.y * m;
            result.sum += s;
            result.e0 += entry.x * s;
            result.sigsq += entry.xsq * s;
          }
        }
      } else {
        // pre-calculated discrete mask w/ linear interpolation
        int i0 = m_quasyEnergyDistribution->getIndex(Ecp);
        double fract = -(Ecp - m_quasyEnergyDistribution->getX(i0)) / m_quasyEnergyDistribution->step;
        if (fract < 0) {
          i0++;
          fract += 1;
        }
        int N = m_quasyEnergyDistribution->entries.size() - 1;
        int M = mask.size() - 1;
        int i1 = std::max(i0 - M, 0);
        int i2 = std::min(i0 + M - 1, N);
        for (int i = i1; i <= i2; i++) {
          const auto& entry = m_quasyEnergyDistribution->entries[i];
          double m = mask[std::abs(i - i0)] * (1 - fract) + mask[std::abs(i - i0 + 1)] * fract;
          double s = entry.y * m;
          result.sum += s;
          result.e0 += entry.x * s;
          result.sigsq += entry.xsq * s;
        }
      }
    }

◆ isAboveThreshold()

bool isAboveThreshold ( ) const

inline

Returns above Cerenkov threshold flag which is set in the prepare method.

Returns: true, if beta is above the Cherenkov threshold for at least one PDE data point

Definition at line 264 of file YScanner.h.

264{return m_aboveThreshold;}

◆ isScanDone()

bool isScanDone ( ) const

inline

Checks which expansion method was used.

Returns: true if scan, false if merge

Definition at line 410 of file YScanner.h.

410{return m_scanDone;}

◆ merge()

void merge	(	unsigned	col,
		double	dydz,
		int	j1,
		int	j2 ) const

private

Performs expansion by merging all reflections.

Parameters

col	pixel column number (0-based)
dydz	photon slope in y-z projection at prism entrance (dy/dz)
j1	first reflection number in y
j2	last (exclusive) reflection number in y

Definition at line 377 of file YScanner.cc.

    {
      int Neven = func::getNumOfEven(j1, j2);
      int Nodd = j2 - j1 - Neven;
 
      const size_t numWindows = m_windowProjections.size();
      const unsigned numRows = m_pixelPositions.getNumPixelRows();
      const double barB = m_bars.front().B;
 
      for (unsigned row = 0; row < numRows; row++) {
 
        int pixelID = m_pixelPositions.pixelID(row, col);
        if (not m_pixelMasks.isActive(pixelID)) continue;
 
        const auto& pixel = m_pixelPositions.get(pixelID);
        double Dy0 = 0;
        double Dy1 = 0;
        PixelProjection proj[2];
        for (size_t k = 0; k < numWindows; k++) {
          projectPixel(pixel.yc, pixel.Dy, k, dydz, proj);
          if (proj[0].Dy > 0) Dy0 += proj[0].Dy;
          if (proj[1].Dy > 0) Dy1 += proj[1].Dy;
        }
        if (Dy0 == 0 and Dy1 == 0) continue;
 
        double Dy = (Dy0 * Neven + Dy1 * Nodd) / (Neven + Nodd);
        Result result(pixelID);
        result.sum = Dy / barB;
        result.e0 = m_meanE;
        result.sigsq = m_rmsE * m_rmsE;
        m_results.push_back(result);
      }
    }

◆ pixelEfficiencies()

PixelEfficiencies & pixelEfficiencies ( )

inlineprivate

Returns non-const pixel relative efficiencies.

Returns: pixel relative efficiencies

Definition at line 508 of file YScanner.h.

508{return m_pixelEfficiencies;}

◆ pixelMasks()

PixelMasks & pixelMasks ( )

inlineprivate

Returns non-const pixel masks.

Returns: pixel masks

Definition at line 502 of file YScanner.h.

502{return m_pixelMasks;}

◆ prepare()

void prepare	(	double	momentum,
		double	beta,
		double	length ) const

Prepare for the PDF expansion in y for a given track mass hypothesis.

Sets photon energy and quasy-energy distributions, mean and r.m.s of photon energy, number of photons and aboveThreshold flag.

Parameters

momentum	particle momentum
beta	particle beta
length	length of particle trajectory within the quartz

Definition at line 149 of file YScanner.cc.

    {
      clear();
 
      m_momentum = momentum;
      m_beta = beta;
      m_length = length;
 
      // check for Cherenkov threshold, return if below
 
      const auto* topgp = TOPGeometryPar::Instance();
      if (beta * topgp->getPhaseIndex(m_meanE0) < 1) return;
 
      // set photon energy distribution, and the mean and r.m.s of photon energy
 
      auto area = setEnergyDistribution(beta);
      if (area == 0) return;
 
      // set number of Cerenkov photons per azimuthal angle per centimeter
 
      m_numPhotons = 370 * area / (2 * M_PI);
 
      // set multiple scattering and surface roughness sigmas in photon energy units
 
      const double radLength = 12.3; // quartz radiation length [cm]
      double thetaScat = 13.6e-3 / beta / momentum * sqrt(length / 2 / radLength); // r.m.s of multiple scattering angle
 
      double n = topgp->getPhaseIndex(m_meanE);
      if (beta * n < 1) {
        B2ERROR("TOP::YScanner::prepare: beta * n < 1 ==> must be a bug!");
        return;
      }
      double dndE = topgp->getPhaseIndexDerivative(m_meanE);
      double dEdTheta = n * sqrt(pow(beta * n, 2) - 1) / dndE;
      m_sigmaScat = std::abs(thetaScat * dEdTheta); // r.m.s of multiple scattering angle converted to photon energy
      m_sigmaAlpha = std::abs(m_bars.back().sigmaAlpha * dEdTheta); // surface roughness converted to photon energy
 
      // set photon energy distribution convoluted with multiple scattering
 
      setQuasyEnergyDistribution(m_sigmaScat);
 
      m_aboveThreshold = true;
    }

◆ projectPixel()

void projectPixel	(	double	yc,
		double	size,
		int	k,
		double	dydz,
		PixelProjection	proj[2] ) const

inlineprivate

Calculates projections of a pixel to prism entrance window (going down-stream the photon).

Parameters

yc	Pixel center in y.
size	Pixel size in y.
k	Valid index of vector of unfolded prism exit windows.
dydz	Photon slope at prism entrance; dydz is used for even projections and -dydz is used for odd projections.
proj	Projections of a pixel to prism entrance (results)

Definition at line 444 of file YScanner.h.

      {
        const auto& win = m_windowProjections[k];
        double halfSize = win.evenReflection ? size / 2 : -size / 2;
        const double ypix[2] =  {yc - halfSize, yc + halfSize}; // pixel edges in y
        double yproj[2][2] = {{0}}; // pixel projections to prism entrance window (second index corresponds to pixel edges)
 
        #pragma omp simd
        for (int i = 0; i < 2; ++i) {
          /* Formerly YScanner::prismEntranceY. */
          double z = ypix[i] * win.sz + win.z0;
          double y = ypix[i] * win.sy + win.y0;
          double dy = dydz * (m_prismZR - z);
          yproj[0][i] = y + dy; // even reflections
          yproj[1][i] = y - dy; // odd reflections
        }
 
        double Bh = m_halfBarThickness;
        for (int i = 0; i < 2; ++i) {
          yproj[i][0] = std::max(yproj[i][0], -Bh);
          yproj[i][1] = std::min(yproj[i][1], Bh);
          proj[i].yc = (yproj[i][0] + yproj[i][1]) / 2;
          proj[i].Dy = yproj[i][1] - yproj[i][0];
        }
      }

◆ scan()

void scan	(	unsigned	col,
		double	yB,
		double	dydz,
		const Derivatives &	D,
		int	j1,
		int	j2 ) const

private

Performs expansion w/ the scan over reflections.

Parameters

col	pixel column number (0-based)
yB	unfolded coordinate y of photon at prism entrance (= Bar exit) plane
dydz	photon slope in y-z projection at prism entrance (dy/dz)
D	the derivatives
j1	first reflection number in y
j2	last (exclusive) reflection number in y

Definition at line 310 of file YScanner.cc.

    {
 
      const size_t numWindows = m_windowProjections.size();
      const unsigned numRows = m_pixelPositions.getNumPixelRows();
      const double barB = m_bars.front().B;
 
      std::map<int, EnergyMask*> masks;
      for (unsigned row = 0; row < numRows; row++) {
 
        int pixelID = m_pixelPositions.pixelID(row, col);
        if (not m_pixelMasks.isActive(pixelID)) continue;
 
        const auto& pixel = m_pixelPositions.get(pixelID);
        std::vector<PixelProjection> projections[2];
        PixelProjection proj[2];
        for (size_t k = 0; k < numWindows; k++) {
          projectPixel(pixel.yc, pixel.Dy, k, dydz, proj);
          if (proj[0].Dy > 0) projections[0].push_back(proj[0]);
          proj[1].yc = -proj[1].yc;
          if (proj[1].Dy > 0) projections[1].push_back(proj[1]);
        }
        if (projections[0].empty() and projections[1].empty()) continue;
 
        for (unsigned k = 0; k < 2; k++) {
          std::sort(projections[k].begin(), projections[k].end());
          for (auto& projection : projections[k]) {
            int iDy = func::lround(projection.Dy * 1000);
            auto& mask = masks[iDy];
            if (not mask) {
              double Dy = projection.Dy;
              double step = m_quasyEnergyDistribution->step;
              mask = new EnergyMask(D.dyB_de, D.dyB_dL, D.dyB_dx, Dy, m_length, pixel.Dx, step);
            }
            projection.mask = mask;
          }
        }
 
        double Ecp_old = 0;
        double wid_old = 1000;
        m_results.push_back(Result(pixelID));
        for (int j = j1; j < j2; j++) {
          double ybar = j * barB - yB;
          for (const auto& projection : projections[std::abs(j) % 2]) {
            double Ecp = (ybar + projection.yc) / D.dyB_de + m_meanE;
            double wid = projection.mask->getFullWidth();
            if (std::abs(Ecp - Ecp_old) > (wid + wid_old) / 2 and m_results.back().sum > 0) {
              m_results.push_back(Result(pixelID));
            }
            integrate(projection.mask, Ecp, m_results.back());
            Ecp_old = Ecp;
            wid_old = wid;
          }
        }
 
        if (m_results.back().sum == 0) m_results.pop_back();
      }
 
      for (auto& result : m_results) result.set();
 
      for (const auto& mask : masks) {
        if (mask.second) delete mask.second;
      }
 
    }

◆ setEnergyDistribution()

double setEnergyDistribution ( double beta ) const

private

Sets photon energy distribution and mean photon energy according to nominal PDE and particle beta.

Parameters

beta	particle beta

Returns: integral of distribution before normalization

Definition at line 194 of file YScanner.cc.

    {
      const auto* topgp = TOPGeometryPar::Instance();
 
      m_energyDistribution.set(m_efficiency);
 
      double s = 0;
      double se = 0;
      double see = 0;
      for (const auto& entry : m_efficiency.entries) {
        double e = entry.x;
        double p = std::max(entry.y * (1 - 1 / pow(beta * topgp->getPhaseIndex(e), 2)), 0.0);
        double ee = entry.xsq;
        m_energyDistribution.entries.push_back(TableEntry(p, e, ee));
        s += p;
        se += p * e;
        see += p * ee;
      }
      if (s == 0) return 0;
 
      for (auto& entry : m_energyDistribution.entries) entry.y /= s;
 
      m_meanE = se / s;
      m_rmsE = sqrt(std::max(see / s - m_meanE * m_meanE, 0.0));
 
      return s * m_energyDistribution.step;
    }

◆ setMirrorCenter()

void setMirrorCenter	(	double	xc,
		double	yc )

inherited

Sets the mirror center-of-curvature.

Parameters

xc	center of curvature in x
yc	center of curvature in y

Definition at line 100 of file RaytracerBase.cc.

    {
      m_mirror.xc = xc;
      m_mirror.yc = yc;
 
      double zc = m_mirror.zc;
      double R = m_mirror.R;
      double zb = zc + R;
      double Ah = m_bars.back().A / 2;
      double Bh = m_bars.back().B / 2;
      zb = std::min(zb, zc + sqrt(pow(R, 2) - pow(xc - Ah, 2) - pow(yc - Bh, 2)));
      zb = std::min(zb, zc + sqrt(pow(R, 2) - pow(xc + Ah, 2) - pow(yc - Bh, 2)));
      zb = std::min(zb, zc + sqrt(pow(R, 2) - pow(xc - Ah, 2) - pow(yc + Bh, 2)));
      zb = std::min(zb, zc + sqrt(pow(R, 2) - pow(xc + Ah, 2) - pow(yc + Bh, 2)));
      m_mirror.zb = zb;
    }

◆ setQuasyEnergyDistribution()

void setQuasyEnergyDistribution ( double sigma ) const

private

Sets photon energy distribution convoluted with a normalized Gaussian.

Parameters

sigma width of the Gaussian [eV]

Definition at line 223 of file YScanner.cc.

    {
      if (m_quasyEnergyDistributions.size() > 1000) {
        m_quasyEnergyDistributions.clear();
        B2ERROR("TOP::YScanner:setQuasyEnergyDistribution: unexpectedly large size of the std::map found, map cleared");
      }
 
      double step = m_energyDistribution.step;
      int ng = func::lround(3 * sigma / step);
      auto& quasyEnergyDistribution = m_quasyEnergyDistributions[ng];
 
      if (quasyEnergyDistribution.entries.empty()) {
        std::vector<double> gaus;
        for (int i = 0; i <= ng; i++) {
          double x = step * i / sigma;
          gaus.push_back(exp(-0.5 * x * x));
        }
 
        quasyEnergyDistribution.set(m_energyDistribution.getX(-ng), step);
        int N = m_energyDistribution.entries.size();
        const auto* entries = m_energyDistribution.entries.data();
        const double* gausTable = gaus.data();
        double sum = 0;
        for (int k = -ng; k < N + ng; k++) {
          // outside of this range getY(k - i) is zero, so those terms contribute nothing
          const int iMin = std::max(-ng, k - N + 1);
          const int iMax = std::min(ng, k);
          double s = 0;
          double se = 0;
          double see = 0;
          for (int i = iMin; i <= iMax; i++) {
            double p = gausTable[std::abs(i)] * entries[k - i].y;
            double e = entries[k - i].x; // equals m_energyDistribution.getX(k - i)
            s += p;
            se += p * e;
            see += p * e * e;
          }
          if (s > 0) {
            se /= s;
            see /= s;
          }
          quasyEnergyDistribution.entries.push_back(TableEntry(s, se, see));
          sum += s;
        }
        for (auto& entry : quasyEnergyDistribution.entries) entry.y /= sum;
      }
 
      m_quasyEnergyDistribution = &quasyEnergyDistribution;
    }

◆ setScanLimits()

static void setScanLimits ( int maxReflections )

inlinestatic

Sets parameters for selection between expand methods.

Parameters

maxReflections maximal number of reflections in y to perform scan

Definition at line 240 of file YScanner.h.

      {
        s_maxReflections = maxReflections;
      }

◆ tabulateSurfaceReflectivity()

double tabulateSurfaceReflectivity ( unsigned n ) const

private

Extends the table of surface reflectivity powers up to n and returns the value.

Called by getSurfaceReflectivity only when the table is not long enough.

Parameters

n	power, i.e. the number of reflections

Returns: reflectivity to the power of n

Definition at line 116 of file YScanner.cc.

    {
      double reflectivity = m_bars.front().reflectivity;
 
      // an unreasonably large number of reflections must not blow up the table
      if (n > s_maxTabulatedPower) return pow(reflectivity, n);
 
      while (m_surfaceReflectivities.size() <= n) {
        m_surfaceReflectivities.push_back(pow(reflectivity, m_surfaceReflectivities.size()));
      }
      return m_surfaceReflectivities[n];
    }

Friends And Related Symbol Documentation

◆ TOPRecoManager

friend class TOPRecoManager

friend

Definition at line 547 of file YScanner.h.

Member Data Documentation

◆ m_aboveThreshold

bool m_aboveThreshold = false

mutableprivate

true if beta is above the Cerenkov threshold

Definition at line 538 of file YScanner.h.

◆ m_bars

std::vector<BarSegment> m_bars

protectedinherited

geometry data of bar segments

Definition at line 188 of file RaytracerBase.h.

◆ m_beta

double m_beta = 0

mutableprivate

particle beta

Definition at line 527 of file YScanner.h.

◆ m_cosTotal

double m_cosTotal = 0

private

cosine of total reflection angle

Definition at line 523 of file YScanner.h.

◆ m_efficiency

Table m_efficiency

private

nominal photon detection efficiencies (PDE)

Definition at line 515 of file YScanner.h.

◆ m_energyDistribution

Table m_energyDistribution

mutableprivate

photon energy distribution

Definition at line 534 of file YScanner.h.

◆ m_geometry

EGeometry m_geometry = c_Unified

protectedinherited

quartz geometry

Definition at line 185 of file RaytracerBase.h.

◆ m_halfBarThickness

double m_halfBarThickness = 0

private

half thickness of the bar at prism entrance

Definition at line 520 of file YScanner.h.

◆ m_length

double m_length = 0

mutableprivate

length of particle trajectory inside quartz

Definition at line 528 of file YScanner.h.

◆ m_meanE

double m_meanE = 0

mutableprivate

mean photon energy

Definition at line 530 of file YScanner.h.

◆ m_meanE0

double m_meanE0 = 0

private

mean photon energy for beta = 1

Definition at line 521 of file YScanner.h.

◆ m_mirror

Mirror m_mirror

protectedinherited

spherical mirror geometry data

Definition at line 189 of file RaytracerBase.h.

◆ m_moduleID

int m_moduleID = 0

protectedinherited

slot ID

Definition at line 184 of file RaytracerBase.h.

◆ m_momentum

double m_momentum = 0

mutableprivate

particle momentum magnitude

Definition at line 526 of file YScanner.h.

◆ m_numPhotons

double m_numPhotons = 0

mutableprivate

number of photons per Cerenkov azimuthal angle per track length

Definition at line 529 of file YScanner.h.

◆ m_optics

EOptics m_optics = c_SemiLinear

protectedinherited

spherical mirror optics

Definition at line 186 of file RaytracerBase.h.

◆ m_pixelEfficiencies

PixelEfficiencies m_pixelEfficiencies

private

pixel relative efficiencies

Definition at line 514 of file YScanner.h.

◆ m_pixelMasks

PixelMasks m_pixelMasks

private

pixel masks

Definition at line 513 of file YScanner.h.

◆ m_pixelPositions

PixelPositions m_pixelPositions

private

positions and sizes of pixels

Definition at line 512 of file YScanner.h.

◆ m_prism

Prism m_prism

protectedinherited

prism geometry data

Definition at line 190 of file RaytracerBase.h.

◆ m_prismZR

double m_prismZR = 0

private

z of the prism-bar joint (copy of m_prism.zR)

Definition at line 519 of file YScanner.h.

◆ m_quasyEnergyDistribution

Table* m_quasyEnergyDistribution = nullptr

mutableprivate

a pointer to the element in m_quasyEnergyDistributions

Definition at line 537 of file YScanner.h.

◆ m_quasyEnergyDistributions

std::map<int, Table> m_quasyEnergyDistributions

mutableprivate

photon energy distributions convoluted with Gaussian of different widths

Definition at line 536 of file YScanner.h.

◆ m_results

std::vector<Result> m_results

mutableprivate

results of PDF expansion in y

Definition at line 541 of file YScanner.h.

◆ m_rmsE

double m_rmsE = 0

mutableprivate

r.m.s of photon energy

Definition at line 531 of file YScanner.h.

◆ m_rmsE0

double m_rmsE0 = 0

private

r.m.s of photon energy for beta = 1

Definition at line 522 of file YScanner.h.

◆ m_scanDone

bool m_scanDone = false

mutableprivate

true if scan performed, false if reflections just merged

Definition at line 542 of file YScanner.h.

◆ m_sigmaAlpha

double m_sigmaAlpha = 0

mutableprivate

surface roughness parameter in photon energy units

Definition at line 533 of file YScanner.h.

◆ m_sigmaScat

double m_sigmaScat = 0

mutableprivate

r.m.s.

of multiple scattering angle in photon energy units

Definition at line 532 of file YScanner.h.

◆ m_surfaceReflectivities

std::vector<double> m_surfaceReflectivities

mutableprivate

bar surface reflectivity to the power of the index; filled on demand, never cleared

Definition at line 518 of file YScanner.h.

◆ m_windowProjections

std::vector<WindowProjection> m_windowProjections

private

pixel projection constants of unfolded prism exit windows

Definition at line 516 of file YScanner.h.

◆ s_maxReflections

int s_maxReflections = 16

staticprivate

maximal number of reflections to perform scan

Definition at line 544 of file YScanner.h.

◆ s_maxTabulatedPower

unsigned s_maxTabulatedPower = 1024

staticprivate

maximal power of the surface reflectivity that is tabulated

Definition at line 545 of file YScanner.h.

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

top/reconstruction_cpp/include/YScanner.h
top/reconstruction_cpp/src/YScanner.cc

Classes

Public Types

Public Member Functions

Static Public Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Static Private Attributes

Friends

Detailed Description

Member Enumeration Documentation

◆ EGeometry

◆ EOptics

Constructor & Destructor Documentation

◆ YScanner()

Member Function Documentation

◆ clear()

◆ expand()

◆ getBars()

◆ getBeta()

◆ getCosTotal()

◆ getEfficiencies()

◆ getEnergyDistribution()

◆ getGeometry()

◆ getMeanEnergy()

◆ getMeanEnergyBeta1()

◆ getMirror()

◆ getModuleID()

◆ getMomentum()

◆ getNumPhotons()

◆ getNumPhotonsPerLen()

◆ getOptics()

◆ getPixelEfficiencies()

◆ getPixelMasks()

◆ getPixelPositions()

◆ getPrism()

◆ getQuasyEnergyDistributions()

◆ getResults()

◆ getRMSEnergy()

◆ getRMSEnergyBeta1()

◆ getSigmaAlpha()

◆ getSigmaScattering()

◆ getSurfaceReflectivity()

◆ getTrackLengthInQuartz()

◆ integrate()

◆ isAboveThreshold()

◆ isScanDone()

◆ merge()

◆ pixelEfficiencies()

◆ pixelMasks()

◆ prepare()

◆ projectPixel()

◆ scan()

◆ setEnergyDistribution()

◆ setMirrorCenter()

◆ setQuasyEnergyDistribution()

◆ setScanLimits()

◆ tabulateSurfaceReflectivity()

Friends And Related Symbol Documentation

◆ TOPRecoManager

Member Data Documentation

◆ m_aboveThreshold

◆ m_bars

◆ m_beta

◆ m_cosTotal

◆ m_efficiency

◆ m_energyDistribution

◆ m_geometry

◆ m_halfBarThickness

◆ m_length

◆ m_meanE

◆ m_meanE0

◆ m_mirror

◆ m_moduleID

◆ m_momentum

◆ m_numPhotons

◆ m_optics

◆ m_pixelEfficiencies

◆ m_pixelMasks

◆ m_pixelPositions