////////////////////////////////////////////////////////////////////
/// \class RAT::DU::LightPathCalculator
///
/// \brief Calculates the refracted (or straight) light path distances
///        between two points within the detector
///
/// \author Rob Stainforth <r.stainforth@liv.ac.uk>
///
/// REVISION HISTORY:\n
///  -   02/2013  : Rob Stainforth - First Revision, new file.
///  -   12/2013  : Rob Stainforth - Updated to include ELLIE reflections
///  - 2014-03-01 : P Jones - refactor as part of ds review.
///  - 2014-04-28 : Rob Stainforth - rename class from LightPath -> LightPathCalculator
///  - 2015-01-23 : Rob Stainforth - Revised path calculation
///  - 2015-02-24 : Rob Stainforth - Error handling for passed 'nan' or 'inf' values
///  - 2015-03-05 : Rob Stainforth - Changed warning messages to debug statements
///  - 2015-03-05 : M Mottram - Updated warnings for partial fill geometry.
///  - 2015-04-20 : Rob Stainforth - Fix divide by zero check in DTheta* methods.
///  - 2015-04-28 : M Mottram - Fixed partial fill calculations; skip refraction
///                             calculations if using straight paths.
///  - 2016-08-29 : M Mottram - Made some members static to speed up copy constructors.
///  - 2017-11-03 : Logan Lebanowski - Modified solid angle calculations [DocDB 4714]
///  - 2019-04-11 : Logan Lebanowski - Neck fixes and tune-up [DocDB 5722]
///  - 2019-06-17 : Logan Lebanowski - Add GetPathValid() and augment error handling
///  - 2019-07-24 : E Leming - Autocall CalcByPositionPartial() when using partial geo
///  - 2022-06-13 : J Tseng - use Point3D for coordinate system consistency
///
/// \details Returns the refracted path through the scintillator
/// AV and water of the detector region. Currently requires single
/// value of the wavelength in MeV - if not specified defaults to
/// 400 nm equivalent.
/// If total internal reflection
/// is present at one of the material interfaces, then the
/// straight line calculation of the path is performed instead.
///
/// It can be used in RAT, ROOT or external programs that link to
/// the libRATEvent library.
///
////////////////////////////////////////////////////////////////////

#ifndef __RAT_DU_LightPathCalculator__
#define __RAT_DU_LightPathCalculator__

#include <RAT/DB.hh>

#include <TGraph.h>
#include <TObject.h>
#include <TVector3.h>
#include <RAT/DU/Point3D.hh>

#include <map>
#include <iostream>

namespace RAT
{
  namespace DU
  {

    // Light Path types
    // Type 'SAW' - Scint -> AV -> Water -> PMT
    // Type 'AW' - AV -> Water -> PMT
    // Type 'ASAW' - AV -> Scint -> AV -> Water -> PMT
    // Type 'WASAW' - Water -> AV -> Scint -> AV -> Water -> PMT
    // Type 'WAW' - Water -> AV -> Water -> PMT
    // Type 'W' - Water -> PMT
    // Type 'WRefl' - Water -> Reflection -> PMT (Reflection off of the AV)
    // Type 'Null' - Light Path Uninitialised
    enum eLightPathType { SAW, AW, ASAW, WASAW, WAW, W, WRefl, Null };

    class LightPathCalculator : public TObject
    {

    public:

      /////////////////////////////////
      ////////     METHODS     ////////
      /////////////////////////////////

      /// Default constructor
      LightPathCalculator() : TObject() { }

      /// Called at the start of a run, loads from the database
      /// Initialise the inner and outer AV radii,
      /// inner and outer neck radii, PMT Bucket radii,
      /// fill fractions (partial fill) and the refractive indices
      void BeginOfRun();

      /// Initalise the refractive indices from the database
      ///
      /// @param[in] dbTable Link to the RATDB table
      /// @param[in,out] property The TGraph to store the refractive indices in
      void LoadRefractiveIndex( DBLinkPtr dbTable,
                                TGraph& property );

      /// Reset/Initalise all the values of the private member variables. Variables
      /// are set to values for which there is no physical interpretation.
      void ResetValues();

      /// Use this to calculate the path of the light from 'eventPos' to 'pmtPos'.
      /// Refraction is modelled for values of 'localityVal' > 0.0. If 'localityVal' = 0.0,
      /// then the straight line approximation is used. For refractive calculations, the refraction
      /// is modelled based on the energy of the light provided (default 3.103125 * 1e-6 MeV = 400nm).
      ///
      /// @param[in] eventPos The starting point of the light path (typically an event position)
      /// @param[in] pmtPos The end point of the PMT (typically a PMT position)
      /// @param[in] energy The photon energy in MeV (defaults to 400 nm -> 3.103125 * 1e-6 MeV)
      /// @param[in] localityVal The accepted tolerance [mm] for how close the path is calculated to the 'pmtPos' (defaults to 0.0 -> Straight Line Calculation)
      void CalcByPosition( const Point3D& eventPos,
                           const Point3D& pmtPos,
                           const Double_t energyMeV = 3.103125 * 1e-6,
                           const Double_t localityVal = 0.0 );

      /// Used for partial fill geometry.
      /// Use this to calculate the path of the light from 'eventPos' to 'pmtPos'.
      /// Refraction is modelled for values of 'localityVal' > 0.0. If 'localityVal' = 0.0,
      /// then the straight line approximation is used. For refractive calculations, the refraction
      /// is modelled based on the energy of the light provided (default 3.103125 * 1e-6 MeV = 400nm).
      ///
      /// @param[in] eventPos The starting point of the light path (typically an event position)
      /// @param[in] pmtPos The end point of the PMT (typically a PMT position)
      /// @param[in] energy The photon energy in MeV (defaults to 400 nm -> 3.103125 * 1e-6 MeV)
      /// @param[in] localityVal The accepted tolerance [mm] for how close the path is calculated to the 'pmtPos' (defaults to 0.0 -> Straight Line Calculation)
      void CalcByPositionPartial( const Point3D& eventPos,
                                  const Point3D& pmtPos,
                                  const Double_t energy = 3.103125 * 1e-6,
                                  const Double_t localityVal = 0.0 );

      /// Calculate the solid angle for this light path, as subtended at the start point from the PMT
      /// at the end point of the light path.
      /// A call to one of the 'CalcByPosition' functions must be made first before calling this
      ///
      /// @param[in] pmtNorm The PMT-Bucket normal vector (pointing IN, towards the AV)
      /// @param[in] nVal The nVal-sided polygon superimposed onto the PMT bucket for the calculation,
      /// use nVal = 0 for spherical cap approximation. nVal must be > 2.
      void CalculateSolidAngle( const TVector3& pmtNorm,
                                const UInt_t nVal );

      /// Calculate the cos of the angle, theta, the light path makes with the bucket face.
      /// This must be called following a call to 'CalcByPosition'
      ///
      /// @param[in] pmtID The ID of the PMT for which the CosTheta value needs to be calculated
      /// @param[out] The value of the CosTheta for the light path incident on this PMT
      Double_t CalculateCosThetaPMT( const Int_t pmtID );

      /// Calculate the total Frensel transmission/reflection coefficients.
      /// Use the respective 'getters' to return these values.
      /// N.B. Cannot currently handle paths that pass through the inner AV
      /// in a partially-filled geometry.
      void CalculateFresnelTRCoeff();

      /// Calculate the parallel component of the Fresnel transmission coefficient
      ///
      /// @param[in] incRI The incident progatating medium refractive index
      /// @param[in] refRI The refracted propagating medium refractive index
      /// @param[in] incAngle The incident angle on the surface
      /// @param[out] The parallel component of the Fresnel transmission coefficient
      Double_t CalculateParallelTransmissionCoefficient( const Double_t incRI,
                                                         const Double_t refRI,
                                                         const Double_t incAngle );

      /// Calculate the perpendicular component of the Fresnel transmission coefficient
      ///
      /// @param[in] incRI The incident progatating medium refractive index
      /// @param[in] refRI The refracted propagating medium refractive index
      /// @param[in] incAngle The incident angle on the surface
      /// @param[out] The perpendicular component of the Fresnel transmission coefficient
      Double_t CalculatePerpendicularTransmissionCoefficient( const Double_t incRI,
                                                              const Double_t refRI,
                                                              const Double_t incAngle );

      /////////////////////////////////
      ////////     GETTERS     ////////
      /////////////////////////////////

      /// Return refractive index in scintillator/innerAV for a given wavelength (energy) in MeV
      /// @param[in] energy The wavelength (energy) in MeV
      ///
      /// @return The refractive index in the scintillator for this wavelength (energy)
      Double_t GetInnerAVRI( const Double_t energy ) const { return InterpolateTGraph(fInnerAVRI, energy); }

      /// Return refractive index in the upper target (partial fill) for a given wavelength (energy) in MeV
      /// @param[in] energy The wavelength (energy) in MeV
      ///
      /// @return The refractive index in the upper filled part of the detector for this wavelength (energy)
      Double_t GetUpperTargetRI( const Double_t energy ) const { return InterpolateTGraph(fUpperTargetRI, energy); }

      /// Return refractive index in the lower target (partial fill) for a given wavelength (energy) in MeV
      /// @param[in] energy The wavelength (energy) in MeV
      ///
      /// @return The refractive index in the lower filled part of the detector for this wavelength (energy)
      Double_t GetLowerTargetRI( const Double_t energy ) const { return InterpolateTGraph(fLowerTargetRI, energy); }

      /// Return refractive index in AV for a given wavelength (energy) in MeV
      /// @param[in] energy The wavelength (energy) in MeV
      ///
      /// @return The refractive index in the AV for this wavelength (energy)
      Double_t GetAVRI( const Double_t energy ) const { return InterpolateTGraph(fAVRI, energy); }

      /// Return refractive index in water for a given wavelength (energy) in MeV
      /// @param[in] energy The wavelength (energy) in MeV
      ///
      /// @return The refractive index in the water for this wavelength (energy)
      Double_t GetWaterRI( const Double_t energy ) const { return InterpolateTGraph(fWaterRI, energy); }

      /// Return the loop ceiling value
      /// (i.e. max number of possible iterations to be made for the refracted path calculation)
      ///
      /// @return The maximum number of iterations for the light path calculation
      Double_t GetLoopCeiling() const { return fLoopCeiling; }

      /// Return final loop value used by iterative scheme used to calculate the refracted path
      /// (i.e. the value of 'i' for the final, i-th iteration)
      ///
      /// @return The number of iterations made for the light path calculation
      Double_t GetFinalLoopValue() const { return fFinalLoopSize; }

      /// Return whether total internal reflection was detected for the path. If so, then the values of
      /// the distances in the scintillator, AV and water are from straight line path approximations.
      ///
      /// @return TRUE: The path encountered total internal reflection, FALSE: It didn't
      Bool_t GetTIR() const { return fIsTIR; }

      /// Return whether the given end point was not calculated to within 'localityVal' of the calculated end point.
      /// 'localityVal' being the tolerance in mm for the required end path to the PMT. (see CalcByPosition)
      /// i.e. final path location > localityVal away from PMT position
      ///
      /// @return TRUE: The end position of the path is outside the required tolerance. FALSE: It isn't
      Bool_t GetResvHit() const { return fResvHit; }

      /// Return whether the path parameters calculated by CalcByPos(Partial) are valid.
      ///
      /// @return TRUE: if the path is valid (no TIR, within tolerance, and a valid light path type) FALSE: if not
      Bool_t GetPathValid() const { return ( !fIsTIR && !fResvHit && fLightPathType!=Null ); }

      /// Return whether the light path enters the neck of the AV
      /// If TRUE: then there are also distances calculated for the scintillator, AV and water from
      /// passing through the neck, see 'GetDistInNeckInnerAV', GetDistInNeckAV' and 'GetDistInNeckWater'
      ///
      /// @return TRUE: if the light path enters the neck of the AV FALSE: if not
      Bool_t GetXAVNeck() const { return fXAVNeck; }

      /// Return whether the light path calculation used the straight line method
      ///
      /// @return TRUE: if the light path is a straight line FALSE: if not
      Bool_t GetStraightLine() const { return fStraightLine; }

      /// Return whether the ELLIE reflected distances are required for start positions outside the AV (i.e. in the water)
      /// If TRUE: then it is assumed that PMTs surrounding the starting position in the water (most likely a fibre position)
      /// are hit from reflections off of the AV. PMTs on the far side are assumed to have light paths associated
      /// with direct light traveling through the entire detector. See LightPath on docDB for how the different regions of PMTs
      /// surrounding the fibre position are determined to have occurred from reflections or direct light.
      ///
      /// @return Whether the ELLIE reflected distances are required for start positions outside the AV
      Bool_t GetELLIEEvent() const { return fELLIEEvent; }

      /// Return the current precision for the path
      /// (the 'localityVal' variable used by CalcByPosition prototype 2). This is essentially the tolerance
      /// of how close the final calculated end position of the light path needs to be to the one passed
      /// as an argument to any of the 'CalcByPosition' calls.
      ///
      /// @return The locality value specified for the current light path calculation
      Double_t GetPathPrecision() const { return fPathPrecision; }

      /// @return The wavelength value used for the light path calculation in nm.
      Double_t GetEnergy() const { return fEnergy; }

      /// @return The distance in Scintillator/InnerAV region
      Double_t GetDistInInnerAV() const { return fDistInInnerAV; }

      /// @return The distance in the acrylic of the AV
      Double_t GetDistInAV() const { return fDistInAV; }

      /// @return The distance in the water
      Double_t GetDistInWater() const { return fDistInWater; }

      /// @return The distance in the upper target of the detector (above fill - for partial fill geometry)
      Double_t GetDistInUpperTarget() const { return fDistInUpperTarget; }

      /// @return The distance in the lower target of the detector (below fill - for partial fill geometry)
      Double_t GetDistInLowerTarget() const { return fDistInLowerTarget; }

      /// @return The distance in the scintillator in the neck region
      Double_t GetDistInNeckInnerAV() const { return fDistInNeckInnerAV; }

      /// @return The distance in the acrylic in the neck region
      Double_t GetDistInNeckAV() const { return fDistInNeckAV; }

      /// @return The distance in the water following the path going through the neck region
      Double_t GetDistInNeckWater() const { return fDistInNeckWater; }

      /// @return The total distance on the light path
      Double_t GetTotalDist() { return ( fDistInInnerAV + fDistInAV + fDistInWater ); }

      /// @return The total distance of the light path for a partial fill geometry
      Double_t GetTotalDistPartial() { return ( fDistInUpperTarget + fDistInLowerTarget + fDistInAV + fDistInWater ); }

      /// @return The solid angle as calculated by 'CalculateSolidAngle'
      Double_t GetSolidAngle() const { return fSolidAngle; }

      /// @return The average CosTheta of the incident path on the PMT bucket face following the 'CalculateSolidAngle' call
      Double_t GetCosThetaAvg() const { return fCosThetaAvg; }

      /// @return The total Fresnel transmission coefficient for the light path, following a call to 'CalculateFresnelTRCoeff'
      Double_t GetFresnelTCoeff() const { return fFresnelTCoeff; }

      /// @return The total Fresnel reflectivity coefficient for the light path, following a call to 'CalculateFresnelTRCoeff'
      Double_t GetFresnelRCoeff() const { return fFresnelRCoeff; }

      /// @return the light path start point
      Point3D GetStartPos() const {
        return Point3D(fAVSystemId, fStartPos);
      }

      /// @return The 'required' light path end point
      Point3D GetEndPos() const {
        return Point3D(fAVSystemId, fEndPos);
      }

      /// @return The light path end position as calculated by 'CalcByPosition'
      Point3D GetLightPathEndPos() const {
        return Point3D(fAVSystemId, fLightPathEndPos);
      }

      /// @return The (unit normalised) incident vector at the PMT bucket, going INTO the PMT
      TVector3 GetIncidentVecOnPMT() const { return fIncidentVecOnPMT; }

      /// @return The (unit normalised) initial light vector from the source position
      TVector3 GetInitialLightVec() const { return fInitialLightVec; }

      /// NOTE: For the below points on the AV, depending on the Light Path
      /// type, it may not have 1st, 2nd, 3rd or 4th points on the AV
      /// where the path intersected. This is based on the LightPathType,
      /// see comments under 'fLightPathType' private member declaration
      /// for details. e.g. a path with starts inside the innerAV region
      /// and proceeds outwards to a PMT will only intersect the AV twice.
      /// Alternatively, a path which begins outside the AV could travel
      /// straight through the AV and therefore intersect the AV four times;
      /// twice going in, and twice going out. Points here are denoted by the
      /// interface points between material interfaces in the detector.

      /// @return The first point on the AV where the light path intersects
      Point3D GetPointOnAV1st() const {
        return Point3D(fAVSystemId, fPointOnAV1st);
      }

      /// @return The second point on the AV where the light path intersects
      Point3D GetPointOnAV2nd() const {
        return Point3D(fAVSystemId, fPointOnAV2nd);
      }

      /// @return The third point on the AV where the light path intersects
      Point3D GetPointOnAV3rd() const {
        return Point3D(fAVSystemId, fPointOnAV3rd);
      }

      /// @return The fourth point on the AV where the light path intersects
      Point3D GetPointOnAV4th() const {
        return Point3D(fAVSystemId, fPointOnAV4th);
      }

      /// @return The first point on the AV where the light path intersects
      Point3D GetPointOnNeck1st() const {
        return Point3D(fAVSystemId, fPointOnNeck1st);
      }

      /// @return The second point on the AV where the light path intersects
      Point3D GetPointOnNeck2nd() const {
        return Point3D(fAVSystemId, fPointOnNeck2nd);
      }

      /// @return The light path type
      std::string GetLightPathType() { return fLightPathTypeMap[ fLightPathType ]; }

      /// @return The incident vector (locally) at the first incident surface on the AV
      TVector3 GetIncidentVecOn1stSurf() const { return ( fPointOnAV1st - fStartPos ).Unit(); }

      /// @return The incident vector (locally) at the second incident surface on the AV
      TVector3 GetIncidentVecOn2ndSurf() const { return ( fPointOnAV2nd - fPointOnAV1st ).Unit(); }

      /// @return The incident vector (locally) at the third incident surface on the AV
      TVector3 GetIncidentVecOn3rdSurf() const { return ( fPointOnAV3rd - fPointOnAV2nd ).Unit(); }

      /// @return The incident vector (locally) at the fourth incident surface on the AV
      TVector3 GetIncidentVecOn4thSurf() const { return ( fPointOnAV4th - fPointOnAV3rd ).Unit(); }

      /// @return The inner AV radius
      Double_t GetAVInnerRadius() const { return fAVInnerRadius; }

      /// @return The outer AV radius
      Double_t GetAVOuterRadius() const { return fAVOuterRadius; }

      /// @return The inner AV Neck radius
      Double_t GetNeckInnerRadius() const { return fNeckInnerRadius; }

      /// @return The outer AV neck radius
      Double_t GetNeckOuterRadius() const { return fNeckOuterRadius; }

      /// @return The PMT bucket radius
      Double_t GetPMTRadius() const { return fPMTRadius; }

      /////////////////////////////////
      ////////     SETTERS     ////////
      /////////////////////////////////

      /// Set the starting position of the light path
      //void SetStartPos( const TVector3& startPos ) { fStartPos = startPos; }

      /// Set the end position of the light path
      //void SetEndPos( const TVector3& endPos ) { fEndPos = endPos; }

      /// For calculations where the event position is in the water region/PSUP
      /// Calculations of the light path distance can be calculated assuming it first
      /// reflected off of the AV first. In which case SetELLIEEvent must be passed
      /// a TRUE bool. ( Default: False - No reflection off of AV )
      ///
      /// @param[in] reflect TRUE: ELLIE reflected distances required FALSE: Not required
      void SetELLIEEvent( const Bool_t reflect )
      {
        fELLIEEvent = reflect;
      }

      // This ROOT macro adds dictionary methods to this class.
      // The number is 0 as this class is never, and should never be written to disc.
      // It assumes this class has no virtual methods, use ClassDef if change this.
      ClassDefNV( LightPathCalculator, 0 );

    private:

      //////////////////////////////////////////////////
      ////////     PRIVATE UTILITY ROUTINES     ////////
      //////////////////////////////////////////////////

      /// Utility Routines for refraction between [Scint/InnerAV]/ AV / Water
      /// 1-3: Inside and Outside light path types
      /// 4-5: Outside light path types only
      /// ThetaResidual: The difference between fPMTTargetTheta and the sum
      /// of ( Theta1st + Theta2nd + Theta3rd )
      /// or ( Theta1st + Theta2nd + Theta3rd + Theta4th + Theta5th )

      /// For a typical path there are various angles defined. A path is calculated
      /// in the 2D-plane containing BOTH the start(source) position and the end(pmt)
      /// position. This ensures that the path calculated is the minimum refracted
      /// path.
      /// To begin, a specific cordinate system is set up for each path, the (x,y,z) directions of this
      /// coordinate system are defined as follows:

      /// x : The direction defined by the radial vector pointing from the origin (centre of AV)
      ///     to the source position.
      /// z : The direction perpendicular to both the radial vector from the origin (centre of AV)
      ///     to the source position, and the radial vector from the origin (centre of AV) to the
      ///     end position. Mathematically speaking, the z-unit direction defines the 2D-plane
      ///     in which the path is calculated
      /// y : The direction defined by the cross product of 'z CROSS x'. The y direction
      ///     is therefore perpendicular to both x and z directions, but lies in the same
      ///     2D-plane as the x direction.

      /// These are utility rountines which calculate the angle of refraction between
      /// the material boundaries in the detector. For a typical path whose start position is
      /// inside the AV, there will be three sets of angles:

      /// 1. The angle between the source position and the scint/innerAV and AV interface point
      ///    as viewed from the origin (centre of AV).
      /// 2. The angle between the first interface point (see above [1.]) and the AV/Water interface
      ///    point as viewed from the origin (centre of AV).
      /// 3. The angle between the second interface point (see above [2.]) and the end position
      ///    of the path.

      /// In the following, 'theta' is the test value for the initial direction of the path
      /// which is minimised against in RTSafe(). It is the initial and defining angle which
      /// ultimate determines the paths course throughout the rest of the detector and consequently
      /// the values of Theta1, Theta2, Theta3 etc. which follow as a result.

      /// The angle between the source position and the first AV intersection point
      /// as viewed from the origin (centre of AV)
      ///
      /// @param[in] theta (passed by reference).
      /// @return The angle between the source position and the first intersection point
      /// for this value of theta.
      Double_t Theta1st( const Double_t theta );

      /// The derivative on this first angle (Theta1) with respect to 'theta'
      ///
      /// @param[in] theta (passed by reference).
      /// @return The derivative on this first angle (Theta1) with respect to 'theta'
      Double_t DTheta1st( const Double_t theta );

      /// The angle between the first AV intersection point and the second AV interseciton point
      /// as viewed from the origin (centre of AV)
      ///
      /// @param[in] theta (passed by reference).
      /// @return The angle between the first AV intersection point and the second AV interseciton
      /// point as viewed from the origin (centre of AV)
      Double_t Theta2nd( const Double_t theta );

      /// The derivative on this second angle (Theta2) with respect to 'theta'
      ///
      /// @param[in] theta (passed by reference).
      /// @return The derivative on this second angle (Theta2) with respect to 'theta'
      Double_t DTheta2nd( const Double_t theta );

      /// The angle between the second AV intersection point and either the PMT position (source
      /// inside the AV) or the third intersection point (source outside the AV).
      ///
      /// @param[in] theta (passed by reference).
      /// @return The angle between the second AV intersection point and either the
      /// PMT position (source inside the AV) or the third intersection point
      /// (source outside the AV).
      Double_t Theta3rd( const Double_t theta );

      /// The derivative on this third angle (Theta3) with respect to 'theta'
      ///
      /// @param[in] theta (passed by reference).
      /// @return The derivative on this third angle (Theta3) with respect to 'theta'
      Double_t DTheta3rd( const Double_t theta );

      /// The angle between the third AV intersection point and the fourth
      /// intersection point (source outside the AV).
      ///
      /// @param[in] theta (passed by reference).
      /// @return The angle between the third AV intersection point and the fourth
      /// intersection point (source outside the AV).
      Double_t Theta4th( const Double_t theta );

      /// The derivative on this fourth angle (Theta4) with respect to 'theta'
      ///
      /// @param[in] theta (passed by reference).
      /// @return The derivative on this fourth angle (Theta4) with respect to 'theta'
      Double_t DTheta4th( const Double_t theta );

      /// The angle between the fourth AV intersection point and the fifth
      /// intersection point (source outside the AV).
      ///
      /// @param[in] theta (passed by reference).
      /// @return The angle between the third AV intersection point and the fifth
      /// intersection point (source outside the AV).
      Double_t Theta5th( const Double_t theta);

      /// The derivative on this fifth angle (Theta5) with respect to 'theta'
      ///
      /// @param[in] theta (passed by reference).
      /// @return The derivative on this fifth angle (Theta5) with respect to 'theta'
      Double_t DTheta5th( const Double_t theta);

      /// Calculate the residual between the target angle between the source and PMT position
      /// and the calculated value.
      ///
      /// @param[in] theta (passed by reference).
      /// @return The residual between the target angle between the source and PMT position
      /// and the calculated value.
      Double_t ThetaResidual( const Double_t theta );

      /// Calculate the derivative on this residual with respect to 'theta'
      ///
      /// @param[in] theta (passed by reference).
      /// @return The derivative on this residual with respect to 'theta'
      Double_t DThetaResidual( const Double_t theta );

      /// Utility function used by 'RTSafe()' to perform the minimiation for
      /// the optimal value of 'theta'.
      ///
      /// @param[in] theta The test value of 'theta' for this path.
      /// @param[in,out] funcVal Computation of 'ThetaResidual()' based on 'theta'
      /// @param[in,out] dFuncVal The derivative on the above value with respect to 'theta'
      void FuncD( Double_t theta, Double_t &funcVal, Double_t &dFuncVal );

      /// Using a combination of Newton-Raphson and bisection methods,
      /// this function returns the root of the function 'LightPathCalculator::Func'
      /// defined on the domain [x1, x2] to within an acceptable accuracy of +/- xAcc
      /// It returns the minimised (optimal) value of 'theta'
      ///
      /// @param[in] x1 The minimum possible value of 'theta' required.
      /// @param[in] x2 The maximum possible value of 'theta' required.
      /// @param[in] xAcc The acceptable tolerance allowed on this value of 'theta'
      /// @return The minimised, optimal value of 'theta' for this path.
      Double_t RTSafe( Double_t x1, Double_t x2, Double_t xAcc );


      /// Calculates the distances for light paths which start inside
      /// the AV. Uses AV coordinates.
      ///
      /// @param[in] eventPos The starting location of the light path (inside the AV)
      /// @param[out] pmtPos The end location of the light path (outside the AV)
      Bool_t CalculateDistancesInnerAV( const TVector3& eventPos,
                                        const TVector3& pmtPos );

      /// Calculates the distances for light paths which start outside
      /// of the AV. Uses AV coordinates.
      ///
      /// @param[in] eventPos The starting location of the light path (outside of the AV)
      /// @param[out] pmtPos The end location of the light path (outside the AV)
      Bool_t CalculateDistancesOutsideAV( const TVector3& eventPos,
                                          const TVector3& pmtPos );

      /// Calculate the refracted path. Performs most of the work required to
      /// obtain a refracted path
      ///
      /// @param[in] initOffset The initial unit vector pointing from the starting position for the light path direction
      void PathCalculation( const TVector3& initOffset );

      /// Readjust the initial photon direction. Used if the previous path
      /// does not meet the locality conditions
      ///
      /// @param[in] distWater distWater The point in the water where the light path EXITS OUT OF (AV coordinates)
      /// @param[in,out] initOffset The initial unit vector pointing from the starting position for the light path direction
      void ReadjustOffset( const TVector3& distWater,
                           TVector3& initOffset );

      /// Test the locality conditions for the hypothesised path end point
      /// against the actual 'required' light path end point (usually a PMT position )
      ///
      /// @param[in] i The final i-th value of the iteration that successfully calculated the light path
      ///
      /// @return TRUE: if the light path is close to the required position FALSE: if not
      Bool_t LocalityCheck( const Int_t i );

      /// Calculate the maximum angle between the event position
      /// and the path direction for the path to intersect with the
      /// sphere of radius 'edgeRadius'.
      /// Uses AV coordinates.
      ///
      /// @param[in] eventPos The starting point of the light path (typically an event position)
      /// @param[in] edgeRadius The radius of the nearest sphere to the event position
      ///
      /// @return Calculate the closest angular displacement of a path close to a surface interface
      Double_t ClosestAngle( const TVector3& eventPos,
                             const Double_t edgeRadius );

      /// Calculate the maximum allowed angle between the event position
      /// and the PMT position for it to reflect off of the AV.
      /// Uses AV coordinates.
      ///
      /// @param[in] eventPos The starting point of the light path (typically an event position)
      /// @param[in] edgeRadius The radius of the nearest sphere to the event position
      ///
      /// @return Calculate the maximum allowed angle between the event position and the PMT position for it to reflect off of the AV
      Double_t ReflectionAngle( const TVector3& eventPos,
                                const Double_t edgeRadius );


      /// Calculate refracted photon vector (unit normalised)
      ///
      /// @param[in] incidentVec The unit vector incident on the surface surface
      /// @param[in] incientSurfVec The unit surface vector of the incident surface
      /// @param[in] incRIndex The incident refractive index
      /// @param[in] refRIndex The refractive index of the refractive media
      ///
      /// @return The refracted unit vector
      TVector3 PathRefraction( const TVector3& incidentVec,
                               const TVector3& incidentSurfVec,
                               const Double_t incRIndex,
                               const Double_t refRIndex );

      /// Calculate vector from some initial point ('startPos'), with an initial
      /// starting direction, 'startDir' to the edge of a sphere of given radius ('radiusFromCentre')
      ///
      /// @param[in] startPos starting position within (or outside of) a sphere
      /// @param[in] startDir Unit starting direction vector from startPos
      /// @param[in] radiusFromCentre The radius of the spherical edge required to be calculated
      /// @param[in] outside Whether the starting position is outside of the sphere in question
      ///
      /// @return The vector at the spheres edge
      TVector3 VectorToSphereEdge( const TVector3& startPos,
                                   const TVector3& startDir,
                                   const Double_t radiusFromCentre,
                                   const Bool_t outside );

      /// Calculate vector from some initial point ('startPos'), with an initial
      /// starting direction, 'startDir' to the edge of a cylinder of given radius ('cylinderRadius')
      ///
      /// @param[in] startPos starting position within (or outside of) a cylinder
      /// @param[in] startDir Unit starting direction vector from startPos
      /// @param[in] cylinderBaseOrigin The location of the origin of the base of the cylinder
      /// @param[in] cylinderRadius The radius of the cylinder
      ///
      /// @return The vector at the spheres edge
      TVector3 VectorToCylinderEdge( const TVector3& startPos,
                                     const TVector3& startDir,
                                     const TVector3& cylinderBaseOrigin,
                                     const Double_t cylinderRadius );

      /// Calculate the path through the upper and lower regions of a partial fill geometry
      ///
      /// @param[in] enterPos The entry point at the target
      /// @param[in] enterDir The entry direction at the target
      /// @param[in,out] exitPos The exit position of the target
      /// @param[in,out] exitDir The exit direction at the exit point of the target
      void PathThroughTarget( const TVector3& enterPos,
                              const TVector3& enterDir,
                              TVector3& exitPos,
                              TVector3& exitDir );

      /// Calculate the refracted path. Performs most of the work required to
      /// obtained a refracted path
      /// @param[in] initialDir The initial direction of the partial path
      ///
      /// @return The partial path
      void PathCalculationPartial( const TVector3& initialDir );

      /// Set the AV neck variables (i.e. those pertaining to whether the
      /// path entered the neck or not, and if it did, the distance it traveled in the scintillator, acrylic and water ).
      /// This is called for all light paths, which, at some point, leave the AV region
      ///
      /// @param[in] pointOnAV The point on the inner AV where the light path enters
      /// @param[in] dirVec The (unit normalised) direction vector of the light path at the above 'pointOnAV'
      void SetAVNeckInformation( const TVector3& pointOnAV,
                                 const TVector3& dirVec );

      /// Calculates the solid angle using a more rigorous method.
      /// This can be generally used by specifying an 'nVal' value to the 'CalculateSolidAngle' call (nVal > 2).
      ///
      /// @param[in] pmtNorm The PMT bucket normal (unit normalised) pointing INTO the detector, towards the AV
      /// @param[in] nVal The nVal-sided polygon to be super-imposed on the PMT bucket and used as part of the solidangle calculation
      void CalculateSolidAnglePolygon( const TVector3& pmtNorm,
                                       const UInt_t nVal );

      /// Calculates the Fresnel transmission/reflectivity coefficients for the light path, assumed 50:50 polarisation
      /// of the light (spherically and linearly)
      ///
      /// @param[in] dir The (unit normalised) direction vector incident on the material surface
      /// @param[in] norm The (unit normalised) surface normal vector pointing in the opposite direction from the incident vector above
      /// @param[in] n1 The incident material refractive index
      /// @param[in] n2 The refracted material refractive index
      /// @param[in,out] T The calculated transmission coefficient at this boundary
      /// @param[in,out] R The calculated reflectivity coefficient at this boundary
      void FresnelTRCoeff( const TVector3& dir,
                           const TVector3& norm,
                           const Double_t n1,
                           const Double_t n2,
                           Double_t& T,
                           Double_t& R );


      /// Set if Total Internal Reflection is detected
      ///
      /// @param[in] val TRUE: The path calculation encountered total internal reflection FALSE: It didn't
      void SetTIR( const Bool_t val ) { fIsTIR = val; }

      /// Set if the calculated path was difficult to resolve.
      /// (i.e. whether the path > fPathPrecision mm away from the end path)
      ///
      /// @param[in] val TRUE: The end calculated position of the path is far from the required positions FALSE: It wasn't
      void SetResvHit( const Bool_t val ) { fResvHit = val; }

      // Wrapper function to avoid TGraph using extrapolation beyond the bounding points.
      // It's possible that using the default extrapolation returns unphysical (negative) values
      /// @param[in] graph The graph containing xy datasets loaded from the db
      /// @param[in] energy The wavelength (energy) in MeV
      ///
      /// @return The refractive index in the scintillator for this wavelength (energy)
      Double_t InterpolateTGraph( const TGraph& graph, const Double_t energy) const{
          if (graph.GetN() == 0) {
            // Cannot use a graph that's empty!
            Log::Die("LightPathCalculator::InterpolateTGraph: Cannot interpolate an empty TGraph.");
          }

          Double_t* energy_values = graph.GetX();
          if(energy < energy_values[0])
              return energy_values[0];
          if(energy > energy_values[graph.GetN()-1])
              return energy_values[graph.GetN()-1];
          return graph.Eval( energy );
      }

      ///////////////////////////////////////////////////
      /////////     PRIVATE MEMBER VARIABLES     ////////
      ///////////////////////////////////////////////////

      /// All position TVector3's are in AV coordinates
      size_t fAVSystemId; // AV coordinate system id (set at BeginOfRun)

      Double_t fNeckInnerRadius;                                         ///< Radius of the inner neck region
      Double_t fNeckOuterRadius;                                         ///< Radius of the outer neck
      Double_t fAVInnerRadius;                                           ///< Radius of the scint region
      Double_t fAVOuterRadius;                                           ///< Radius of the AV region
      Double_t fPMTRadius;                                               ///< Radius of the PMT bucket
      Double_t fFillZ;                                                   ///< z position of the partial fill

      static TGraph fInnerAVRI;                                          ///< Scintillator refractive index
      static TGraph fUpperTargetRI;                                      ///< The 'filled' region of the detector refractive index
      static TGraph fLowerTargetRI;                                      ///< The 'un-filled' region of the detector refractive index
      static TGraph fAVRI;                                               ///< AV refractive index
      static TGraph fWaterRI;                                            ///< Water refractive index

      Double_t fLoopCeiling;                                             ///< Iteration Ceiling for algortithm loop
      Double_t fFinalLoopSize;                                           ///< Final loop value which meets locality conditions
      Double_t fPathPrecision;                                           ///< The accepted path proximity/tolerance to the PMT location [mm]

      Double_t fInnerAVRIVal;                                            ///< The value of the scintillator refractive index used for this path
      Double_t fUpperTargetRIVal;                                        ///< The value of the upper target volume index used for this path (partial fill)
      Double_t fLowerTargetRIVal;                                        ///< The value of the lower target volume index used for this path (partial fill)
      Double_t fAVRIVal;                                                 ///< The value of the AV refractive index used for this path
      Double_t fWaterRIVal;                                              ///< The value of the water refractive index used for this path

      TVector3 fIncidentVecOnPMT;                                        ///< Final light path direction (unit normalised)
      TVector3 fInitialLightVec;                                         ///< Initial light path direction (unit normalised)

      TVector3 fStartPos;                                                ///< Start position of the light path (AV coordinates)
      TVector3 fEndPos;                                                  ///< Required end position of the light path (AV coords)
      TVector3 fLightPathEndPos;                                         ///< Calculated end position of the light path (AV)

      Double_t fPMTTargetTheta;                                          ///< The target PMT theta angle for the light path

      Bool_t fIsTIR;                                                     ///< TRUE: Total Internal Reflection encountered FALSE: It wasn't
      Bool_t fResvHit;                                                   ///< TRUE: Difficult path to resolve and calculate FALSE: It wasn't
      Bool_t fXAVNeck;                                                   ///< TRUE: Path entered neck region FALSE: It didn't
      Bool_t fELLIEEvent;                                              ///< TRUE: Reflected distances in water off of AV on PMTs near starting position required
                                                                         ///< FALSE: Reflected distances not required
      Bool_t fStraightLine;                                              ///< TRUE: Light Path is a straight line approximation FALSE: It isn't


      // Note: Depending on the light path type (see 'fLightPathType'), the path may
      // intersect the AV/Neck once, twice or three
      // ...or four times

      TVector3 fPointOnAV1st;                                            ///< Point on AV where light path first hits the AV
      TVector3 fPointOnAV2nd;                                            ///< Point on AV where light path hits the AV a second time
      TVector3 fPointOnAV3rd;                                            ///< Point on AV where light path hits the AV a third time
      TVector3 fPointOnAV4th;                                            ///< Point on AV where light path hits the AV a fourth time

      TVector3 fPointOnNeck1st;                                          ///< Point on the Neck where the light path hits a first time
      TVector3 fPointOnNeck2nd;                                          ///< Point on the Neck where the light path hits a second time
      TVector3 fPointOnNeck3rd;                                          ///< Point on the Neck where the light path hits a third time
      TVector3 fPointOnNeck4th;                                          ///< Point on the Neck where the light path hits a fourth time

      eLightPathType fLightPathType;                                     ///< Light path type, based on what regions of the detector the path enters

      std::map< eLightPathType, std::string > fLightPathTypeMap;         ///< Map containing a descriptor for the light path type

      Double_t fDistInInnerAV;                                           ///< Distance in the scintillator region
      Double_t fDistInUpperTarget;                                       ///< Distance in the upper target region (partial fill geometry)
      Double_t fDistInLowerTarget;                                       ///< Distance in the lower target region (partial fill geometry)
      Double_t fDistInAV;                                                ///< Distance in the acrylic region of the AV
      Double_t fDistInWater;                                             ///< Distance in the water region
      Double_t fDistInNeckInnerAV;                                       ///< Distance through the scintillator region in the neck (if GetXAVNeck() = TRUE)
      Double_t fDistInNeckAV;                                            ///< Distance through the acrylic of the AV region in the neck (if GetXAVNeck() = TRUE)
      Double_t fDistInNeckWater;                                         ///< Distance through the water region in the neck (if GetXAVNeck() = TRUE)

      Double_t fEnergy;                                                  ///< The value of the wavelength in MeV

      Double_t fSolidAngle;                                              ///< The solid angle subtended by the PMT for this light path
      Double_t fCosThetaAvg;                                             ///< Average incident angle on the PMT for this path.
                                                                         ///< This is only calculated after a call to CalculateSolidAngle

      Double_t fFresnelTCoeff;                                           ///< The combined Fresnel TRANSMISSION coefficient for this path
      Double_t fFresnelRCoeff;                                           ///< The combined Fresnel REFLECTIVITY coefficient for this path

      Bool_t fIs_partial_fill;                                           ///< Whether detector is in partial-fill phase or not

    };

  } // namespace DU

} // namespace RAT

#endif