InverseRaytracer.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/PhotonState.h>
#include <top/reconstruction_cpp/TOPTrack.h>
#include <vector>
#include <Math/Point3D.h>
#include <cmath>
 
namespace Belle2 {
  namespace TOP {
 
    class InverseRaytracer : public RaytracerBase {
 
    public:
 
      enum ErrorCodes {
        c_NoPhysicsSolution = -1,   
        c_NoEquationSolution = -2,  
        c_NoReflectionPoint = -3    
      };
 
 
      struct CerenkovAngle {
        double cosThc = 0; 
        double sinThc = 0; 
        CerenkovAngle()
        {}
 
        explicit CerenkovAngle(double cosTheta);
      };
 
 
      struct Solution {
        double cosFic = 0; 
        double sinFic = 0; 
        double kx = 0;     
        double ky = 0;     
        double kz = 0;     
        double len = 0;    
        double xD = 0;     
        double zD = 0;     
        double yD = 0;     
        double yB = 0;     
        int Nxm = 0;       
        int Nym = 0;       
        double step = 0;   
        bool totRefl = true; 
        Solution(double cfi, double sfi);
 
        void setDirection(const CerenkovAngle& cer, const TOPTrack::TrackAngles& trk);
 
        void setTotalReflStatus(double A, double B, double cosTotal);
 
        void setTotalReflStatus(double A, double B, double cosTotal, double Kx, double Ky);
 
        bool getStatus() const;
 
        double getFic() const {return atan2(sinFic, cosFic);}
      };
 
 
      InverseRaytracer(int moduleID, double cosTotal): RaytracerBase(moduleID, c_Unified, c_SemiLinear),
        m_cosTotal(cosTotal)
      {}
 
      static void setMaxPropagationLen(double maxLen) {s_maxLen = maxLen;}
 
      void clear() const;
 
      int solveDirect(double xD, double zD, const TOPTrack::AssumedEmission& assumedEmission,
                      const CerenkovAngle& cer, double step = 0) const;
 
      int solveReflected(double xD, double zD, int Nxm, double xmMin, double xmMax,
                         const TOPTrack::AssumedEmission& assumedEmission,
                         const CerenkovAngle& cer, double step = 0) const;
 
      int solveForReflectionPoint(double xM, int Nxm, const TOPTrack::AssumedEmission& assumedEmission,
                                  const CerenkovAngle& cer) const;
 
      double getCosTotal() const {return m_cosTotal;}
 
      std::vector<Solution>& getSolutions(unsigned i) const {return m_solutions[i % 2];}
 
      bool getStatus() const {return (m_ok[0] or m_ok[1]);}
 
      bool getStatus(unsigned i) const {return m_ok[i % 2];}
 
      bool isNymDifferent() const;
 
      PhotonState getReconstructedPhoton(const Solution& sol, double DFic = 0) const;
 
    private:
 
      bool solve(double dxdz, const CerenkovAngle& cer, const TOPTrack::TrackAngles& trk) const;
 
      bool findReflectionPoint(double xE, double zE, double xD, double zD, double xmMin, double xmMax,
                               double& xM, double& zM, double& dxdz) const;
 
      double getDeltaXE(double x, double xe, double ze, double xd, double zd) const;
 
      double m_cosTotal = 0; 
      mutable std::vector<Solution> m_solutions[2]; 
      mutable bool m_ok[2] = {false, false}; 
      mutable ROOT::Math::XYZPoint m_emiPoint; 
      mutable CerenkovAngle m_cer; 
      mutable TOPTrack::TrackAngles m_trk; 
      static double s_maxLen; 
    };
 
  } // namespace TOP
} // namespace Belle2