DeltaRayPDF.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/BackgroundPDF.h>
#include <top/reconstruction_cpp/PixelPositions.h>
#include <top/reconstruction_cpp/TOPTrack.h>
#include <vector>
#include <cmath>
 
namespace Belle2 {
  namespace TOP {

    class DeltaRayPDF {
 
    public:

      explicit DeltaRayPDF(int moduleID);

      void prepare(const TOPTrack& track, const Const::ChargedStable& hypothesis);

      int getModuleID() const {return m_moduleID;}

      double getNumPhotons() const {return m_numPhotons;}

      double getFraction() const {return m_fraction;}

      double getPDFValue(int pixelID, double time) const;

      double getPDFValue(double time, double dt0 = 0, double acc = 1) const;

      double getIntegral(double minTime, double maxTime) const;
 
    private:

      struct GausXY {
        double x = 0; 
        double y = 0; 

        explicit GausXY(double X): x(X), y(exp(-0.5 * x * x))
        {}
      };

      double angularDistr(double kz) const;

      double timeDistr(double t, double t0) const;

      double smearedTimeDistr(double t, double t0) const;

      double peakFraction(double tmin, double tmax, double t0) const;

      double totalFraction(double tmin, double tmax) const;

      double directFraction(double z) const;

      double photonYield(double beta, int PDGCode) const;
 
      // variables set in constructor (slot dependent)
      int m_moduleID; 
      const BackgroundPDF* m_background = 0; 
      const PixelPositions* m_pixelPositions = 0; 
      double m_zD = 0; 
      double m_zM = 0; 
      double m_phaseIndex = 0; 
      double m_groupIndex = 0; 
      double m_dispersion = 0; 
      double m_angularNorm = 0; 
      std::vector<double> m_norms; 
      std::vector<GausXY> m_tableGaus; 
 
      // variables set in prepare method (track/hypothesis dependent)
      double m_xE = 0;  
      double m_yE = 0;  
      double m_zE = 0;  
      double m_dirFrac = 0; 
      double m_dirT0 = 0;   
      double m_reflT0 = 0;  
      double m_TOF = 0;     
      double m_fraction = 0; 
      double m_numPhotons = 0; 
      std::vector<double> m_pixelAcceptances; 
 
    };
 
 
    //--- inline functions ------------------------------------------------------------
 
    inline double DeltaRayPDF::angularDistr(double kz) const
    {
      if (std::abs(kz) < 1 / m_phaseIndex) {
        double x = sqrt((1 - 1 / pow(m_phaseIndex, 2)) / (1 - pow(kz, 2)));
        return 1 - acos(x) * 4 / M_PI;
      }
      return 1;
    }
 
    inline double DeltaRayPDF::timeDistr(double t, double t0) const
    {
      if (t < t0) return 0;
      return angularDistr(t0 / t) / m_angularNorm * t0 / pow(t, 2);
    }
 
    inline double DeltaRayPDF::getPDFValue(double time, double dt0, double acc) const
    {
      double t = time - m_TOF;
      double pdfDirect = smearedTimeDistr(t, m_dirT0 + dt0) * acc;
      double pdfReflec = smearedTimeDistr(t, m_reflT0);
      return m_dirFrac * pdfDirect + (1 - m_dirFrac) * pdfReflec;
    }
 
    inline double DeltaRayPDF::getIntegral(double minTime, double maxTime) const
    {
      return totalFraction(minTime, maxTime) / m_fraction;
    }
 
  } // namespace TOP
} // namespace Belle2