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() {…}

        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 Derivatives {…};
 

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

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

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

        Table()
        {}
        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 Table {…};
 

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

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

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

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

        void set();
      };
      struct Result {…};
 

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

      static void setScanLimits(int maxReflections)
      {
        s_maxReflections = maxReflections;
      }
      static void setScanLimits(int 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;
 
    };
    class YScanner : public RaytracerBase {…};
 
 
    //--- 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::dLen_d(const InverseRaytracer::Solution& sol0, {…}
 
    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::dyB_d(const InverseRaytracer::Solution& sol0, {…}
 
    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 double YScanner::Derivatives::dFic_d(const InverseRaytracer::Solution& sol0, {…}
 
    inline void YScanner::Table::clear()
    {
      x0 = 0;
      step = 0;
      entries.clear();
    }
    inline void YScanner::Table::clear() {…}
 
    inline void YScanner::Table::set(double X0, double Step)
    {
      x0 = X0;
      step = Step;
      entries.clear();
    }
    inline void YScanner::Table::set(double X0, double Step) {…}
 
    inline void YScanner::Table::set(const Table& T)
    {
      x0 = T.x0;
      step = T.step;
      entries.clear();
    }
    inline void YScanner::Table::set(const Table& T) {…}
 
    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 double YScanner::Table::getY(int i) const {…}
 
    inline void YScanner::Result::set()
    {
      if (sum == 0) return;
      e0 /= sum;
      sigsq = std::max(sigsq / sum - e0 * e0, 0.0);
    }
    inline void YScanner::Result::set() {…}
 
  } // namespace TOP
} // namespace Belle2