YScanner.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <top/reconstruction_cpp/RaytracerBase.h>
#include <top/reconstruction_cpp/InverseRaytracer.h>
#include <top/reconstruction_cpp/PixelPositions.h>
#include <top/reconstruction_cpp/PixelMasks.h>
#include <top/reconstruction_cpp/PixelEfficiencies.h>
#include <top/reconstruction_cpp/EnergyMask.h>
#include <vector>
#include <algorithm>
#include <cmath>
#include <map>
 
 
namespace Belle2 {
  namespace TOP {

    class YScanner : public RaytracerBase {
 
    public:

      struct Derivatives {
        double dLen_dx = 0; 
        double dLen_de = 0; 
        double dLen_dL = 0; 
        double dyB_dx = 0;  
        double dyB_de = 0;  
        double dyB_dL = 0;  
        double dFic_dx = 0; 
        double dFic_de = 0; 
        double dFic_dL = 0; 

        Derivatives()
        {}

        Derivatives(const InverseRaytracer::Solution& sol,
                    const InverseRaytracer::Solution& sol_dx,
                    const InverseRaytracer::Solution& sol_de,
                    const InverseRaytracer::Solution& sol_dL);

        static double dLen_d(const InverseRaytracer::Solution& sol0, const InverseRaytracer::Solution& sol1);

        static double dyB_d(const InverseRaytracer::Solution& sol0, const InverseRaytracer::Solution& sol1);

        static double dFic_d(const InverseRaytracer::Solution& sol0, const InverseRaytracer::Solution& sol1);
      };
 

      struct TableEntry {
        double y = 0;   
        double x = 0;   
        double xsq = 0; 

        TableEntry(double Y, double X, double Xsq):
          y(Y), x(X), xsq(Xsq)
        {}
      };
 

      struct Table {
        double x0 = 0;         
        double step = 0;       
        std::vector<TableEntry> entries; 

        Table()
        {}

        void clear();

        void set(double X0, double Step);

        void set(const Table& T);

        double getX(int i) const;

        double getXmin() const;

        double getXmax() const;

        int getIndex(double x) const;

        double getY(int i) const;
      };
 

      struct PixelProjection {
        double yc = 0; 
        double Dy = 0; 
        const EnergyMask* mask = 0; 

        bool operator<(const PixelProjection& other) const {return yc < other.yc;}
      };
 

      struct WindowProjection {
        double sy = 0; 
        double sz = 0; 
        double y0 = 0; 
        double z0 = 0; 
        bool evenReflection = false; 
      };
 

      struct Result {
        int pixelID = 0;  
        double sum = 0;   
        double e0 = 0;    
        double sigsq = 0; 

        explicit Result(int ID): pixelID(ID)
        {}

        void set();
      };
 

      explicit YScanner(int moduleID, unsigned N = 64);

      static void setScanLimits(int maxReflections)
      {
        s_maxReflections = maxReflections;
      }

      void clear() const;

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

      bool isAboveThreshold() const {return m_aboveThreshold;}

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

      const PixelPositions& getPixelPositions() const {return m_pixelPositions;}

      const PixelMasks& getPixelMasks() const {return m_pixelMasks;}

      const PixelEfficiencies& getPixelEfficiencies() const {return m_pixelEfficiencies;}

      const Table& getEfficiencies() const {return m_efficiency;}

      double getCosTotal() const {return m_cosTotal;}

      double getMomentum() const {return m_momentum;}

      double getBeta() const {return m_beta;}

      double getTrackLengthInQuartz() const {return m_length;}

      double getNumPhotonsPerLen() const {return m_numPhotons;}

      double getNumPhotons() const {return m_numPhotons * m_length;}

      double getMeanEnergy() const {return m_meanE;}

      double getRMSEnergy() const {return m_rmsE;}

      double getMeanEnergyBeta1() const {return m_meanE0;}

      double getRMSEnergyBeta1() const {return m_rmsE0;}

      double getSigmaScattering() const {return m_sigmaScat;}

      double getSigmaAlpha() const {return m_sigmaAlpha;}

      double getSurfaceReflectivity(unsigned n) const
      {
        if (n < m_surfaceReflectivities.size()) return m_surfaceReflectivities[n];
        return tabulateSurfaceReflectivity(n);
      }

      const Table& getEnergyDistribution() const {return m_energyDistribution;}

      const std::map<int, Table>& getQuasyEnergyDistributions() const {return m_quasyEnergyDistributions;}

      const std::vector<Result>& getResults() const {return m_results;}

      bool isScanDone() const {return m_scanDone;}
 
 
    private:

      double setEnergyDistribution(double beta) const;

      void setQuasyEnergyDistribution(double sigma) const;

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

      void projectPixel(double yc, double size, int k, double dydz, PixelProjection proj[2]) const
      {
        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];
        }
      }

      double tabulateSurfaceReflectivity(unsigned n) const;

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

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

      PixelMasks& pixelMasks() {return m_pixelMasks;}

      PixelEfficiencies& pixelEfficiencies() {return m_pixelEfficiencies;}
 
      // variables set in constructor or, for m_surfaceReflectivities, filled lazily on first use
      // (slot dependent, therefore not reset by clear())
      PixelPositions m_pixelPositions; 
      PixelMasks m_pixelMasks; 
      PixelEfficiencies m_pixelEfficiencies; 
      Table m_efficiency; 
      std::vector<WindowProjection> m_windowProjections; 
      mutable std::vector<double>
      m_surfaceReflectivities; 
      double m_prismZR = 0; 
      double m_halfBarThickness = 0; 
      double m_meanE0 = 0; 
      double m_rmsE0 = 0; 
      double m_cosTotal = 0; 
 
      // variables set in prepare method (track/hypothesis dependent)
      mutable double m_momentum = 0; 
      mutable double m_beta = 0; 
      mutable double m_length = 0; 
      mutable double m_numPhotons = 0; 
      mutable double m_meanE = 0; 
      mutable double m_rmsE = 0; 
      mutable double m_sigmaScat = 0; 
      mutable double m_sigmaAlpha = 0; 
      mutable Table m_energyDistribution; 
      mutable std::map<int, Table>
      m_quasyEnergyDistributions; 
      mutable Table* m_quasyEnergyDistribution = nullptr; 
      mutable bool m_aboveThreshold = false; 
 
      // results of expand method (pixel column dependent)
      mutable std::vector<Result> m_results;   
      mutable bool m_scanDone = false;  
 
      static int s_maxReflections; 
      static unsigned s_maxTabulatedPower; 
 
      friend class TOPRecoManager;
 
    };
 
 
    //--- inline functions ------------------------------------------------------------
 
    inline double YScanner::Derivatives::dLen_d(const InverseRaytracer::Solution& sol0,
                                                const InverseRaytracer::Solution& sol1)
    {
      return (sol1.len - sol0.len) / sol1.step;
    }
 
    inline double YScanner::Derivatives::dyB_d(const InverseRaytracer::Solution& sol0,
                                               const InverseRaytracer::Solution& sol1)
    {
      return (sol1.yB - sol0.yB) / sol1.step;
    }
 
    inline double YScanner::Derivatives::dFic_d(const InverseRaytracer::Solution& sol0,
                                                const InverseRaytracer::Solution& sol1)
    {
      if (std::abs(sol0.cosFic) > std::abs(sol0.sinFic)) {
        return (sol1.sinFic - sol0.sinFic) / sol0.cosFic / sol1.step;
      } else {
        return -(sol1.cosFic - sol0.cosFic) / sol0.sinFic / sol1.step;
      }
    }
 
    inline void YScanner::Table::clear()
    {
      x0 = 0;
      step = 0;
      entries.clear();
    }
 
    inline void YScanner::Table::set(double X0, double Step)
    {
      x0 = X0;
      step = Step;
      entries.clear();
    }
 
    inline void YScanner::Table::set(const Table& T)
    {
      x0 = T.x0;
      step = T.step;
      entries.clear();
    }
 
    inline double YScanner::Table::getX(int i) const {return x0 + step * i;}
 
    inline double YScanner::Table::getXmin() const {return x0;}
 
    inline double YScanner::Table::getXmax() const {return x0 + step * (entries.size() - 1);}
 
    inline int YScanner::Table::getIndex(double x) const {return func::lround((x - x0) / step);}
 
    inline double YScanner::Table::getY(int i) const
    {
      unsigned k = i;
      if (k >= entries.size()) return 0;
      return entries[k].y;
    }
 
    inline void YScanner::Result::set()
    {
      if (sum == 0) return;
      e0 /= sum;
      sigsq = std::max(sigsq / sum - e0 * e0, 0.0);
    }
 
  } // namespace TOP
} // namespace Belle2