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 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;}
 
      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;
 
      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 (slot dependent)
      PixelPositions m_pixelPositions; 
      PixelMasks m_pixelMasks; 
      PixelEfficiencies m_pixelEfficiencies; 
      Table m_efficiency; 
      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; 
      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 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