#include <CLHEP/Units/PhysicalConstants.h>
using CLHEP::pi; using CLHEP::twopi; using CLHEP::halfpi;

#include <RAT/DU/ShadowingCalculator.hh>
#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/Utility.hh>
using namespace RAT;
using namespace RAT::DU;

#include <TVector3.h>
#include <TMath.h>

#include <algorithm>
#include <string>

#include <iostream>
using namespace std;

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::BeginOfRun()
{

  // First clear member vectors
  fAVPipe1Points.clear();
  fAVPipe2Points.clear();
  fAVPipe3Points.clear();
  fAVPipe4Points.clear();
  fAVPipeLabSPoints.clear();
  fAVPipe6Points.clear();
  fAVPipeLabRPoints.clear();
  fAVPipePoints.clear();
  fNCDAnchorPoints.clear();
  fBellyPlateCentrePoints.clear();
  fBellyPlatePoints.clear();
  fAVHDRopeSlingPoints.clear();
  fAVHDRopeSlingRolls.clear();
  fAVHDRopeEquatorPoints.clear();
  fAVHDRopeUpperAVPoints.clear();
  fAVHDRopeLowerAVPoints.clear();
  fSupportRopeEquatorPoints.clear();
  fSupportRopeUpperAVPoints.clear();

  DB *db = DB::Get();
  fAVInnerRadius = db->GetLink( "SOLID", "acrylic_vessel_inner" )->GetD( "r_sphere" );
  fAVOuterRadius = db->GetLink( "SOLID", "acrylic_vessel_outer" )->GetD( "r_sphere" );
  fNeckOuterRadius = db->GetLink( "SOLID", "acrylic_vessel_outer" )->GetD( "r_neck" );

  // Try and load the AV hold-down rope information.
  try{
    LoadAVHDRopeInfo();
    fExistAVHDRope = true;
    warn << "ShadowingCalculator: Loaded AV hold-down rope shadowing information.\n";
  }
  catch( DBNotFoundError& e ){
    warn << "ShadowingCalculator: Cannot load AV hold-down rope parameters, these will not be loaded.\n";
    fExistAVHDRope = false;
  }

  // Try and load the support rope information.
  try{
    LoadSupportRopeInfo();
    fExistSupportRope = true;
    warn << "ShadowingCalculator: Loaded AV support rope shadowing information.\n";
  }
  catch( DBNotFoundError& e ){
    warn << "ShadowingCalculator: Cannot load support rope parameters, these will not be loaded.\n";
    fExistSupportRope = false;
  }

  // Try and load the belly plate information.
  try{
    LoadBellyPlateInfo();
    fExistBellyPlate = true;
    warn << "ShadowingCalculator: Loaded belly plate shadowing information.\n";
  }
  catch( DBNotFoundError& e ){
    warn << "ShadowingCalculator: Cannot load belly plate parameters, these will not be loaded.\n";
    fExistBellyPlate = false;
  }

  // Try and load the NCD anchor information.
  try{
    LoadNCDAnchorInfo();
    fExistNCDAnchor = true;
    warn << "ShadowingCalculator: Loaded NCD anchor shadowing information.\n";
  }
  catch( DBNotFoundError& e ){
    warn << "ShadowingCalculator: Cannot load NCD anchor parameters, these will not be loaded.\n";
    fExistNCDAnchor = false;
  }

  // Try and load the AV pipe information.
  try{
    LoadAVPipeInfo();
    fExistAVPipe = true;
    warn << "ShadowingCalculator: Loaded AV pipe shadowing information.\n";
  }
  catch( DBNotFoundError& e ){
    warn << "ShadowingCalculator: Cannot load AV pipe parameters, these will not be loaded.\n";
    fExistAVPipe = false;
  }

  // Try and load the neck boss information.
  try{
    LoadNeckBossInfo();
    fExistNeckBoss = true;
    warn << "ShadowingCalculator: Loaded neck boss shadowing information.\n";
  }
  catch( DBNotFoundError& e ){
    warn << "ShadowingCalculator: Cannot neck boss parameters, these will not be loaded.\n";
    fExistNeckBoss = false;
  }

  // Now set the default values for the tolerances about which
  // paths are considered shadowed for the individual components
  // of the geometry. All defaults equivalent to a proximity of
  // 150.0 mm. These can be set by the user
  // by using the Set*Tolerance() methods.
  fAVHDRopeTolerance = 150.0;
  fSupportRopeTolerance = 150.0;
  fBellyPlateTolerance = 150.0;
  fNCDAnchorTolerance = 150.0;
  fAVPipeTolerance = 150.0;
  fNeckBossTolerance = 150.0;

  // Set all the shadowing flags to be false to begin with
  fAVHDRopeShadowFlag = false;
  fSupportRopeShadowFlag = false;
  fBellyPlateShadowFlag = false;
  fNCDAnchorShadowFlag = false;
  fAVPipeShadowFlag = false;
  fNeckBossShadowFlag = false;

}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForAVHDRopeShadowing( RAT::DU::LightPathCalculator& lPath ){

  Bool_t avhdShadow = false;

  // Initialise values for the point on the AV to check for shadowing
  // and the starting direction from which to propagate the path from.
  TVector3 avPoint( 0.0, 0.0, 0.0 );
  TVector3 startDir( 0.0, 0.0, 0.0 );

  // For paths going OUT OF the AV, we only need to check one point on the AV
  // and one direction. If the path is going accross the detector, there will be a
  // second AV intersection point and direction vector which needs to be associated
  // with the path.
  TVector3 avPoint2( 0.0, 0.0, 0.0 );
  TVector3 startDir2( 0.0, 0.0, 0.0 );

  // Declare the radius of the ropes and a tolerance
  // equivalent to about 170.0 mm = 150.0 mm [=default fAVHDRopeTolerance] + rope radius (~20.0 mm)
  const Double_t ropeRadius = fAVHDRopeRadius;
  const Double_t ropeTolerancePhi = ( ropeRadius + fAVHDRopeTolerance ) / ( fAVHDRopeEquatorPoints[ 0 ].Mag() );

  // Determine the light path type in order to identify which intersection point with the
  // AV to check with, and also the starting direction to propagate the path
  // from to check for rope shadowing. For some light path types there may be more than one
  // crossing of the AV, and so we will need to perform two checks in these cases.

  string lightPathType = lPath.GetLightPathType();

  // One check required.
  if ( lightPathType == "Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV2nd();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // One check required.
  else if ( lightPathType == "AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV3rd();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // One check required.
  else if ( lightPathType == "AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV4th();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV2nd();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // This is a special case and is catered for later on in this method.
  else if ( lightPathType == "Water" ){
    startDir = lPath.GetInitialLightVec();
    startDir2 = -1.0 * lPath.GetInitialLightVec();
  }

  // Two checks required.
  else if ( lightPathType == "Water->Reflection->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV1st();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // Light uninitialised, throw error.
  else if ( lightPathType == "Light Path Uninitialised" ){
    warn << "ShadowingCalculator::CheckAVHDRopeShadowing: LightPathCalculator object is uninitialised, call LightPathCalculator::CalcByPosition first before checking for shadowing!\n";
    return false;
  }

  if ( lightPathType != "Water"
       && lightPathType != "LightPath Uninitialised" ){

    // Check to see what part of the rope net we need to check against.
    // If below the equator, then need to check for the anchored part of the rope net
    // using cylindrical symmetry
    if ( avPoint.Z() < fAVHDRopeEquatorPoints[ 0 ].Z() ){

      TVector3 cylinderOrigin( 0.0, 0.0, 0.0 );
      TVector3 pointOnCylinder = VectorToCylinderEdge( avPoint, startDir, cylinderOrigin, fAVHDRopeEquatorPoints[ 0 ].Mag() );

      // Find an example vector to a location on the ropes which pass from inside to outside the PSUP.
      TVector3 zAxisDown( 0.0, 0.0, -1.0 );
      TVector3 zToPSUP = VectorToSphereEdge( fAVHDRopeEquatorPoints[ 0 ], zAxisDown, 8500.0, 0.0 );

      // If the point on the cylinder is outside the PSUP, then there is no shadowing.
      if ( -1.0 * TMath::Abs( pointOnCylinder.Z() ) < zToPSUP.Z() ){ avhdShadow = false; }

      // Otherwise, we need to check with each angular location of the ropes on
      // a cylindrical surface symmetrical about the centre.
      else{

        for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){

          // Create vectors for four locations on the AV equator where a single
          // rope sling lies.
          TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
          TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

          eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
          eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

          pointOnCylinder.SetZ( 0.0 );

          // If the location on this imaginary cylinder is within an angular tolerance of
          // any of the rope points, then set the shadowing boolean to TRUE.
          if ( (TMath::Abs( pointOnCylinder.Angle( eqPointPair ) ) <  ropeTolerancePhi)
               || (TMath::Abs( pointOnCylinder.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ avhdShadow = true; }
        }
      }
      // If we have encountered shadowing already, return the boolean value. Otherwise
      // we may stil have shadowing above the equator and therefore need to proceed to the 'else'
      // statement below.
      if ( avhdShadow == true ){ return true; }
    }

    // We are above the equator, and so we just need to check the AV intersection points with the set
    // of points which charaterise the AV hold-down rope net.
    else{

      for ( Size_t iPoint = 0; iPoint < fAVHDRopeUpperAVPoints.size(); iPoint++ ){

        if ( ( avPoint - fAVHDRopeUpperAVPoints[ iPoint ] ).Mag() < fAVHDRopeRadius + fAVHDRopeTolerance ){ avhdShadow = true; }

      }
      if ( avhdShadow == true ){ return true; }
    }

    // If the light path type only involves a single intersection point, or regardless of this
    // it has found shadowing in any light path type, then we are done and can return the shadowing bool.
    if ( ( lightPathType == "Scint/InnerAV->AV->Water"
           || lightPathType == "AV->Scint/InnerAV->AV->Water"
           || lightPathType == "AV->Water" ) ){ return avhdShadow; }

    // Otherwise, for other light path types, we need to check the second AV intersection points
    // and their corresponding directions to check whether they propagate to a rope position.
    // As noted above though, we only reach this part if we haven't already encountered
    // shadowing in the first instances of the below path types.
    if ( ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water"
           || lightPathType == "Water->AV->Water"
           || lightPathType == "Water->Reflection->Water" ) ){

      // Check whether the second point is above or below the equator.
      if ( avPoint2.Z() < fAVHDRopeEquatorPoints[ 0 ].Z() ){

        TVector3 cylinderOrigin2( 0.0, 0.0, 0.0 );
        TVector3 pointOnCylinder2 = VectorToCylinderEdge( avPoint2, startDir2, cylinderOrigin2, fAVHDRopeEquatorPoints[ 0 ].Mag() );

        // Find an example vector to a location on the ropes which pass from inside to outside the PSUP.
        TVector3 zAxisDown( 0.0, 0.0, -1.0 );
        TVector3 zToPSUP = VectorToSphereEdge( fAVHDRopeEquatorPoints[ 0 ], zAxisDown, 8500.0, 0.0 );

        // If the point on the cylinder is outside the PSUP, then there is no shadowing.
        if ( -1.0 * TMath::Abs( pointOnCylinder2.Z() ) < zToPSUP.Z() ){ avhdShadow = false; }

        // Otherwise, we need to check with each angular location of the ropes on
        // a cylindrical surface symmetrical about the centre.
        else{

          for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){
            TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
            TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

            eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
            eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

            pointOnCylinder2.SetZ( 0.0 );

            if ( (TMath::Abs( pointOnCylinder2.Angle( eqPointPair ) ) <  ropeTolerancePhi)
                 || (TMath::Abs( pointOnCylinder2.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ avhdShadow = true; }
          }
        }
        if ( avhdShadow == true ){ return true; }
      }

      // We are above the equator, and so we just need to check the AV intersection points with the set
      // of points which charaterise the AV hold-dwon rope net.
      else{

        for ( Size_t iPoint = 0; iPoint < fAVHDRopeUpperAVPoints.size(); iPoint++ ){

          if ( ( avPoint2 - fAVHDRopeUpperAVPoints[ iPoint ] ).Mag() < fAVHDRopeRadius + fAVHDRopeTolerance ){ avhdShadow = true; }

        }
      }
      if ( avhdShadow == true ){ return true; }
      else { return false; }
    }  // End of 'Water->AV->Scint/InnerAV->AV->Water' // 'Water->AV->Water' // 'Water->Reflection->Water'
  }

  // Finally, in the special case that the light path is from one PMT to another (with zero intersection)
  // then we use the following calculation.
  if ( lightPathType == "Water" ){

    TVector3 startPos = lPath.GetStartPos();
    TVector3 endPos = lPath.GetEndPos();

    TVector3 startPosXY = lPath.GetStartPos();
    startPosXY.SetZ( 0.0 );
    TVector3 endPosXY = lPath.GetEndPos();
    endPosXY.SetZ( 0.0 );

    TVector3 origin( 0.0, 0.0, 0.0 );

    // Get vector which points to the point on one of the support ropes where it exits the PSUP
    TVector3 zAxisDown( 0.0, 0.0, -1.0 );
    TVector3 zToPSUP = VectorToSphereEdge( fAVHDRopeEquatorPoints[ 0 ], zAxisDown, 8500.0, 0.0 );

    // If the starting and end position radii (best imagined in the (x,y) plane) are smaller
    // than the radial positions of the support ropes, then no intersection is possible.
    // therefore return false.
    if ( startPosXY.Mag() < fAVHDRopeEquatorPoints[ 0 ].Mag()
         && endPosXY.Mag() < fAVHDRopeEquatorPoints[ 0 ].Mag() ){ return false; }

    // If either one of the starting or end position radii (best imagined in the (x,y) plane) is
    // greater than the radial positions of the hold down ropes, and the other position is smaller,
    // then a single intersection is possible with the ropes.
    else if ( ( startPosXY.Mag() < fAVHDRopeEquatorPoints[ 0 ].Mag()
                && endPosXY.Mag() > fAVHDRopeEquatorPoints[ 0 ].Mag() )
              ||
              ( startPosXY.Mag() > fAVHDRopeEquatorPoints[ 0 ].Mag()
                && endPosXY.Mag() < fAVHDRopeEquatorPoints[ 0 ].Mag() ) ){


      // Check which of the two instances is the case and find the corresponsing point
      // on the imaginary cylindrical surface. Surface of radii 'fAVHDRopeEquatorPoints[ 0 ].Mag()',
      // centred at the 'origin' from a starting posiiton of startPos/endPos with an initial
      // direction startDir/startDir2
      TVector3 pointOnCylinder( 0.0, 0.0, 0.0 );

      if ( startPosXY.Mag() < fAVHDRopeEquatorPoints[ 0 ].Mag() ){
        pointOnCylinder = VectorToCylinderEdge( startPos, startDir, origin,
                                                fAVHDRopeEquatorPoints[ 0 ].Mag() );
      }
      else{
        pointOnCylinder = VectorToCylinderEdge( endPos, startDir2, origin,
                                                fAVHDRopeEquatorPoints[ 0 ].Mag() );
      }


      // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
      // of shadowing, therefore return 'false'.
      if ( -1.0 * TMath::Abs( pointOnCylinder.Z() ) < zToPSUP.Z() ){ return false; }

      // Otherwise check whether the angular position of the point on the cylinder
      // is coincident (within tolerance) with any of the angular distributions
      // of the ropes
      else{

        for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){
          TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
          TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

          eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
          eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

          pointOnCylinder.SetZ( 0.0 );

          if ( (TMath::Abs( pointOnCylinder.Angle( eqPointPair ) ) <  ropeTolerancePhi)
               || (TMath::Abs( pointOnCylinder.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ return true; }

        }
      }
      return false;
    }

    // In the case where both the start and end position radii (best imagined in the (x,y) plane)
    // are outside of the circle defined by the radii of the hold-down rope locations, then there are two
    // possibilities: either they don't intersect the ropes, and pass outside, or they can potentially
    // intersect twice.
    else if ( startPosXY.Mag() > fAVHDRopeEquatorPoints[ 0 ].Mag()
              && endPosXY.Mag() > fAVHDRopeEquatorPoints[ 0 ].Mag() ){

      // Angle threshold for whether the path between the two points coincides with the cylinder
      if ( ClosestAngle( startPosXY, fAVHDRopeEquatorPoints[ 0 ].Mag() ) > startPosXY.Angle( endPosXY ) ){

        TVector3 pointOnCylinder = VectorToCylinderEdge( startPos, startDir, origin,
                                                         fAVHDRopeEquatorPoints[ 0 ].Mag() );
        TVector3 pointOnCylinder2 = VectorToCylinderEdge( endPos, startDir2, origin,
                                                          fAVHDRopeEquatorPoints[ 0 ].Mag() );

        // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
        // of shadowing, therefore set the boolean to 'false' (remember, it could intersect at the second
        // instance, so we don't want to 'return' anything just yet).
        if ( -1.0 * TMath::Abs( pointOnCylinder.Z() ) < zToPSUP.Mag() ){ avhdShadow = false; }
        else{

          // Otherwise check whether the angular position of the point on the cylinder
          // is coincident (within tolerance) with any of the angular distributions
          // of the ropes
          for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){
            TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
            TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

            eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
            eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

            pointOnCylinder.SetZ( 0.0 );

            if ( (TMath::Abs( pointOnCylinder.Angle( eqPointPair ) ) <  ropeTolerancePhi)
                 || (TMath::Abs( pointOnCylinder.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ return true; }

          }
        }

        // For the second possible intersection...

        // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
        // of shadowing, therefore return 'false' (now we can 'return' here).
        if ( TMath::Abs( pointOnCylinder2.Z() ) > zToPSUP.Mag() ){ return false; }
        else{

          for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){
            TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
            TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

            eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
            eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

            pointOnCylinder.SetZ( 0.0 );

            if ( (TMath::Abs( pointOnCylinder2.Angle( eqPointPair ) ) <  ropeTolerancePhi)
                 || (TMath::Abs( pointOnCylinder2.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ return true; }

          }
        }
      }

      // If the path doesn't intersect the cylindrical surface, then there is no chance
      // to intersect. Therefore return false.
      else{
        return false;
      }
    }

  } // End of Light Path Type: 'Water'

    // Warn Message, Light Path Uninitialised or something has gone wrong.
  else{
    warn << "ShadowingCalculator::CheckForAVHDRopeShadowing: Error! Cannot calculate shadowing for this light light path type: ";
    warn << lightPathType << "\n";
  }

  // If somehow we reach here, return false.
  return false;

}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForSupportRopeShadowing( RAT::DU::LightPathCalculator& lPath){

  Bool_t supportRopeShadow = false;

  // Initialise values for the point on the AV to check for shadowing
  // and the starting direction from which to propagate the path from.
  TVector3 avPoint( 0.0, 0.0, 0.0 );
  TVector3 startDir( 0.0, 0.0, 0.0 );

  // For paths going OUT OF the AV, we only need to check one point on the AV
  // and one direction. If the path is going accross the detector, there will be a
  // second AV intersection point and direction vector which needs to be associated
  // with the path.
  TVector3 avPoint2( 0.0, 0.0, 0.0 );
  TVector3 startDir2( 0.0, 0.0, 0.0 );

  const Double_t ropeRadius = fSupportRopeRadius;
  const Double_t ropeTolerancePhi = ( ropeRadius +  fSupportRopeTolerance ) / ( fSupportRopeEquatorPoints[ 0 ].Mag() );

  // Determine the light path type in order to identify which intersection point with the
  // AV to check with, and also the starting direction to propagate the path
  // from to check for rope shadowing. For some light path types there may be more than one
  // crossing of the AV, and so we will need to perform two checks in these cases.

  string lightPathType = lPath.GetLightPathType();

  // One check required.
  if ( lightPathType == "Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV2nd();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // One check required.
  else if ( lightPathType == "AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV3rd();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // One check required.
  else if ( lightPathType == "AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV4th();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV2nd();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // This is a special case and is catered for later on in this method.
  else if ( lightPathType == "Water" ){
    startDir = lPath.GetInitialLightVec();
    startDir2 = -1.0 * lPath.GetInitialLightVec();
  }


  // Two checks required
  else if ( lightPathType == "Water->Reflection->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV1st();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // Uninitialised light path - throw error.
  else if ( lightPathType == "Light Path Uninitialised" ){
    warn << "ShadowingCalculator::CheckForSupportRopeShadowing: LightPath is uninitialised, call CalcByPosition first before checking for shadowing!\n";
    return false;
  }

  if ( lightPathType != "Water"
       && lightPathType != "LightPath Uninitialised" ){

    // If above the equator, then check to see if the light path hits any of the support ropes
    if ( avPoint.Z() > fSupportRopeEquatorPoints[ 0 ].Z() ){

      // Find the point on an imaginary cylindrical surface with radius 'fSupportRopeEquatorPoints[ 0 ].Mag()'
      // centred at 'cylinderOrigin', where the point 'avPoint' propagates to with a direction vector 'startDir'
      TVector3 cylinderOrigin( 0.0, 0.0, 0.0 );
      TVector3 pointOnCylinder = VectorToCylinderEdge( avPoint, startDir, cylinderOrigin, fSupportRopeEquatorPoints[ 0 ].Mag() );

      // Get vector which points to the point on one of the support ropes where it exits the PSUP
      TVector3 zAxis( 0.0, 0.0, 1.0 );
      TVector3 zToPSUP = VectorToSphereEdge( fSupportRopeEquatorPoints[ 0 ], 1.0 * zAxis, 8500.0, 0.0 );

      // If the calculated 'Z' position of the point on the cylinder is above this aforementioned
      // vector location (i.e. it is actually outside of the PSUP) then there is no chance
      // of shadowing, therefore return 'false'.
      if ( TMath::Abs( pointOnCylinder.Z() ) > zToPSUP.Z() ){ supportRopeShadow = false; }
      else{

        // Otherwise check whether the angular position of the point on the cylinder
        // is coincident (within tolerance) with any of the angular distributions
        // of the ropes
        for ( Size_t iEq = 0; iEq < fSupportRopeEquatorPoints.size(); iEq++ ){

          TVector3 eqPoint = fSupportRopeEquatorPoints[ iEq ];
          pointOnCylinder.SetZ( 0.0 );

          if ( ( TMath::Abs( pointOnCylinder.Angle( eqPoint ) ) < ropeTolerancePhi ) ){ supportRopeShadow =  true; }
        }
      }
      // If we have encountered shadowing already, return the boolean value. Otherwise
      // we may have other instances to check against.
      if ( supportRopeShadow == true ){ return true; }
    }

    // If the light path type only involves a single intersection point
    // then we are done and can return the shadowing bool.
    if ( ( lightPathType == "Scint/InnerAV->AV->Water"
           || lightPathType == "AV->Scint/InnerAV->AV->Water"
           || lightPathType == "AV->Water" ) ){ return supportRopeShadow; }

    // Otherwise, for other light path types, we need to check the second AV intersection points
    // and their corresponding directions to check whether they propagate to a rope position.
    // As noted above though, we only reach this part if we haven't already encountered
    // shadowing in the first instances of the below path types.
    if ( ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water"
           || lightPathType == "Water->AV->Water"
           || lightPathType == "Water->Reflection->Water" ) ){

      // Check whether the second point is above or below the equator.
      if ( avPoint2.Z() > fSupportRopeEquatorPoints[ 0 ].Z() ){

        // Find the point on an imaginary cylindrical surface with radius 'fSupportRopeEquatorPoints[ 0 ].Mag()'
        // centred at 'cylinderOrigin2', where the point 'avPoint2' propagates to with a direction vector 'startDir2'
        TVector3 cylinderOrigin2( 0.0, 0.0, 0.0 );
        TVector3 pointOnCylinder2 = VectorToCylinderEdge( avPoint2, startDir2, cylinderOrigin2, fSupportRopeEquatorPoints[ 0 ].Mag() );

        // Get vector which points to the point on one of the support ropes where it exits the PSUP
        TVector3 zAxis( 0.0, 0.0, 1.0 );
        TVector3 zToPSUP = VectorToSphereEdge( fSupportRopeEquatorPoints[ 0 ], 1.0 * zAxis, 8500.0, 0.0 );

        // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
        // of shadowing, therefore return 'false'.
        if ( TMath::Abs( pointOnCylinder2.Z() ) > zToPSUP.Z() ){ supportRopeShadow = false; }
        else{

          // Otherwise check whether the angular position of the point on the cylinder
          // is coincident (within tolerance) with any of the angular distributions
          // of the ropes
          for ( Size_t iEq = 0; iEq < fSupportRopeEquatorPoints.size(); iEq++ ){

            TVector3 eqPoint = fSupportRopeEquatorPoints[ iEq ];
            pointOnCylinder2.SetZ( 0.0 );

            if ( ( TMath::Abs( pointOnCylinder2.Angle( eqPoint ) ) < ropeTolerancePhi ) ){ supportRopeShadow =  true; }
          }
        }
      }
      // If we have encountered shadowing already, return the boolean value.
      if ( supportRopeShadow == true ){ return true; }
      else { return false; }
    } // End of 'Water->AV->Scint/InnerAV->AV->Water' // 'Water->AV->Water' // 'Water->Reflection->Water' instances

  }

  // Finally, in the special case that the light path is from one PMT to another (with zero intersection)
  // then we use the following calculation.
  if ( lightPathType == "Water" ){

    TVector3 startPos = lPath.GetStartPos();
    TVector3 endPos = lPath.GetEndPos();

    TVector3 startPosXY = lPath.GetStartPos();
    startPosXY.SetZ( 0.0 );
    TVector3 endPosXY = lPath.GetEndPos();
    endPosXY.SetZ( 0.0 );

    TVector3 origin( 0.0, 0.0, 0.0 );

    // Get vector which points to the point on one of the support ropes where it exits the PSUP
    TVector3 zAxis( 0.0, 0.0, 1.0 );
    TVector3 zToPSUP = VectorToSphereEdge( fSupportRopeEquatorPoints[ 0 ], zAxis, 8500.0, 0.0 );

    // If the starting and end position radii (best imagined in the (x,y) plane) are smaller
    // than the radial positions of the support ropes, then no intersection is possible.
    // therefore return false.
    if ( startPosXY.Mag() < fSupportRopeEquatorPoints[ 0 ].Mag()
         && endPosXY.Mag() < fSupportRopeEquatorPoints[ 0 ].Mag() ){ return false; }

    // If either one of the starting or end position radii (best imagined in the (x,y) plane) is
    // greater than the radial positions of the support ropes, and the other position is smaller,
    // then a single intersection is possible with the ropes.
    else if ( ( startPosXY.Mag() < fSupportRopeEquatorPoints[ 0 ].Mag()
                && endPosXY.Mag() > fSupportRopeEquatorPoints[ 0 ].Mag() )
              ||
              ( startPosXY.Mag() > fSupportRopeEquatorPoints[ 0 ].Mag()
                && endPosXY.Mag() < fSupportRopeEquatorPoints[ 0 ].Mag() ) ){



      // Check which of the two instances is the case and find the corresponsing point
      // on the imaginary cylindrical surface. Surface of radii 'fSupportRopeEquatorPoints[ 0 ].Mag()',
      // centred at the 'origin' from a starting posiiton of startPos/endPos with an initial
      // direction startDir/startDir2
      TVector3 pointOnCylinder( 0.0, 0.0, 0.0 );

      if ( startPosXY.Mag() < fSupportRopeEquatorPoints[ 0 ].Mag() ){
        pointOnCylinder = VectorToCylinderEdge( startPos, startDir, origin,
                                                fSupportRopeEquatorPoints[ 0 ].Mag() );
      }
      else{
        pointOnCylinder = VectorToCylinderEdge( endPos, startDir2, origin,
                                                fSupportRopeEquatorPoints[ 0 ].Mag() );
      }

      // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
      // of shadowing, therefore return 'false'.
      if ( TMath::Abs( pointOnCylinder.Z() ) > zToPSUP.Z() ){ return false; }
      else{

        // Otherwise check whether the angular position of the point on the cylinder
        // is coincident (within tolerance) with any of the angular distributions
        // of the ropes
        for ( Size_t iEq = 0; iEq < fSupportRopeEquatorPoints.size(); iEq++ ){

          TVector3 eqPoint = fSupportRopeEquatorPoints[ iEq ];
          pointOnCylinder.SetZ( 0.0 );

          if ( ( TMath::Abs( pointOnCylinder.Angle( eqPoint ) ) < ropeTolerancePhi ) ){ supportRopeShadow =  true; }
        }

      }
      if ( supportRopeShadow == true ){ return true; }
      else { return false; }
    }

    // In the case where both the start and end position radii (best imagined in the (x,y) plane)
    // are outside of the circle defined by the radii of the support rope locations, then there are two
    // possibilities: either they don't intersect the ropes, and pass outside, or they can potentially
    // intersect twice.
    else if ( startPosXY.Mag() > fAVHDRopeEquatorPoints[ 0 ].Mag()
            && endPosXY.Mag() > fAVHDRopeEquatorPoints[ 0 ].Mag() ){

      // Angle threshold for whether the path between the two points coincides with the cylinder
      if ( ClosestAngle( startPosXY, fAVHDRopeEquatorPoints[ 0 ].Mag() ) > startPosXY.Angle( endPosXY ) ){

        TVector3 pointOnCylinder = VectorToCylinderEdge( startPos, startDir, origin,
                                                         fAVHDRopeEquatorPoints[ 0 ].Mag() );
        TVector3 pointOnCylinder2 = VectorToCylinderEdge( endPos, startDir2, origin,
                                                          fAVHDRopeEquatorPoints[ 0 ].Mag() );

        // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
        // of shadowing, therefore set the boolean to 'false' (remember, it could intersect at the second
        // instance, so we don't want to 'return' anything just yet).
        if ( TMath::Abs( pointOnCylinder.Z() ) > zToPSUP.Mag() ){ supportRopeShadow = false; }
        else{

          // Otherwise check whether the angular position of the point on the cylinder
          // is coincident (within tolerance) with any of the angular distributions
          // of the ropes
          for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){
            TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
            TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

            eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
            eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

            pointOnCylinder.SetZ( 0.0 );

            if ( (TMath::Abs( pointOnCylinder.Angle( eqPointPair ) ) <  ropeTolerancePhi)
                 || (TMath::Abs( pointOnCylinder.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ return true; }

          }
        }

        // For the second possible intersection...

        // Check vector location (i.e. it is actually outside of the PSUP) then there is no chance
        // of shadowing, therefore return 'false' (now we can 'return' here).
        if ( TMath::Abs( pointOnCylinder2.Z() ) > zToPSUP.Mag() ){ return false; }
        else{

          // Otherwise check whether the angular position of the point on the cylinder
          // is coincident (within tolerance) with any of the angular distributions
          // of the ropes
          for ( Size_t iEq = 0; iEq < fAVHDRopeSlingRolls.size(); iEq++ ){
            TVector3 eqPoint = fAVHDRopeSlingPoints[ 1 ];
            TVector3 eqPointPair = fAVHDRopeSlingPoints[ 4 ];

            eqPointPair.RotateZ( fAVHDRopeSlingRolls[ iEq ] );
            eqPoint.RotateZ( TMath::Pi() + fAVHDRopeSlingRolls[ iEq ]  );

            pointOnCylinder.SetZ( 0.0 );

            if ( (TMath::Abs( pointOnCylinder2.Angle( eqPointPair ) ) <  ropeTolerancePhi)
                 || (TMath::Abs( pointOnCylinder2.Angle( eqPoint ) ) <  ropeTolerancePhi) ){ return true; }

          }
        }
      }

      // If the path doesn't intersect the cylindrical surface, then there is no chance
      // to intersect. Therefore return false.
      else{
        return false;
      }
    }

  } // End of Light Path Type: 'Water'


   // Warn Message, Light Path Uninitialised or something has gone wrong.
  else{
    warn << "ShadowingCalculator::CheckForSupportRopeShadowing: Error! Cannot calculate shadowing for this light light path type: ";
    warn << lightPathType << "\n";
  }

  // If somehow we reach here, return false.
  return false;

}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForBellyPlateShadowing( RAT::DU::LightPathCalculator& lPath ){

  Bool_t bellyPlateShadow = false;

  // Initialise values for the point on the AV to check for shadowing
  // and the starting direction from which to propagate the path from.
  TVector3 avPoint( 0.0, 0.0, 0.0 );
  TVector3 startDir( 0.0, 0.0, 0.0 );

  // For paths going OUT OF the AV, we only need to check one point on the AV
  // and one direction. If the path is going accross the detector, there will be a
  // second AV intersection point and direction vector which needs to be associated
  // with the path.
  TVector3 avPoint2( 0.0, 0.0, 0.0 );
  TVector3 startDir2( 0.0, 0.0, 0.0 );

  // Determine the light path type in order to identify which intersection point with the
  // AV to check with, and also the starting direction to propagate the path
  // from to check for rope shadowing. For some light path types there may be more than one
  // crossing of the AV, and so we will need to perform two checks in these cases.

  string lightPathType = lPath.GetLightPathType();

  // One check required.
  if ( lightPathType == "Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV2nd();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // One check required.
  else if ( lightPathType == "AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV3rd();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // One check required.
  else if ( lightPathType == "AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = lPath.GetIncidentVecOnPMT();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV4th();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV2nd();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // This is a special case and is catered for later on in this method.
  else if ( lightPathType == "Water" ){
    startDir = lPath.GetInitialLightVec();
    startDir2 = -1.0 * lPath.GetInitialLightVec();
  }


  // Two checks required
  else if ( lightPathType == "Water->Reflection->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    startDir = -1.0 * lPath.GetInitialLightVec();
    avPoint2 = lPath.GetPointOnAV1st();
    startDir2 = lPath.GetIncidentVecOnPMT();
  }

  // Uninitialised light path - throw error.
  else if ( lightPathType == "Light Path Uninitialised" ){
    warn << "ShadowingCalculator::CheckForBellyPlateShadowing: LightPath is uninitialised, call CalcByPosition first before checking for shadowing!\n";
    return false;
  }

  // The tolerance in the z-direction of the belly plates.
  // The belly plates are centrered on z=0 mm. So we can check half height
  // + a little extra (fBellyPlateHalfHeight + 150.0 [=default fBellyPlateTolerance]) to see the AV intersection
  // point is coincident with the belly plate location.
  const Double_t bpZTol = fBellyPlateHalfHeight + fBellyPlateTolerance;

  // Angular displacement equivalent to 150.0 mm at the equator point
  const Double_t phiTol = fBellyPlateTolerance / fAVOuterRadius;
  const TVector3 zAxis( 0.0, 0.0, 1.0 );

  if ( lightPathType != "Water"
       && lightPathType != "LightPath Uninitialised" ){

    // If the AV intersection point is above or below (in the z-direction)
    // the bellow plate, then no shadowing will occur, so set the boolean to false.
    if ( avPoint.Z() > bpZTol || avPoint.Z() < -1.0 * bpZTol ){ bellyPlateShadow = false; }

    // Otherwise, in the z-direction we are in the right z-region.
    // so now need to perform angular checks.
    else if ( avPoint.Z() < bpZTol && avPoint.Z() > -1.0 * bpZTol ){

      for ( Size_t iBP = 0; iBP < fBellyPlateCentrePoints.size(); iBP++ ){

        // Get the centre locations for each belly plate.
        TVector3 bpPos = fBellyPlateCentrePoints[ iBP ];

        // Find a vector to the position of the top and bottom of the belly plate
        TVector3 bpTop = bpPos + ( fBellyPlateHalfHeight * zAxis );
        bpTop.SetMag( bpPos.Mag() );
        TVector3 bpBottom = bpPos - ( fBellyPlateHalfHeight * zAxis );
        bpBottom.SetMag( bpPos.Mag() );

        // Then find the angle between the vectors to the top and the bottom
        // of the belly plate and divide by two to find the bisector.
        Double_t bpPhi = bpTop.Angle( bpBottom );
        Double_t startPhi = bpPos.Phi() - ( bpPhi / 2.0 );

        // Next define two vectors left and right of the centre point
        // of the belly plate.
        TVector3 bpLeft = bpPos;
        bpLeft.SetPhi( startPhi - phiTol );
        TVector3 bpRight = bpPos;
        bpRight.SetPhi( startPhi + bpPhi + phiTol );

        TVector3 avPointCpy = avPoint;
        avPointCpy.SetZ( bpPos.Z() );

        // Now check to see if the angle of the AV intersection point
        // is within the angular region of the belly plate centre and
        // the left and right side of the plate itself.
        if ( avPointCpy.Angle( bpPos ) < bpLeft.Angle( bpPos )
             && avPointCpy.Angle( bpPos ) < bpRight.Angle( bpPos ) ){ return true; }

      }
      if ( bellyPlateShadow == true ){ return true; }
    }

    // If the light path type only involves a single AV intersection point
    // then we are done and can return the shadowing bool.
    if ( ( lightPathType == "Scint/InnerAV->AV->Water"
           || lightPathType == "AV->Scint/InnerAV->AV->Water"
           || lightPathType == "AV->Water" ) ){ return bellyPlateShadow; }

    // Otherwise, for other light path types, we need to check the second AV intersection points
    // and their corresponding directions to check whether they coincide with a belly plate position.
    // As noted above though, we only reach this part if we haven't already encountered
    // shadowing in the first instances of the below path types.
    if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water"
         || lightPathType == "Water->AV->Water"
         || lightPathType == "Water->Reflection->Water" ){

      // If the AV intersection point is above or below (in the z-direction)
      // the bellow plate, then no shadowing will occur, so set the boolean to false.
      if ( avPoint2.Z() > bpZTol || avPoint2.Z() < -1.0 * bpZTol ){ bellyPlateShadow = false; }

      // Otherwise, in the z-direction we are in the right z-region.
      // so now need to perform angular checks.
      else if ( avPoint2.Z() < bpZTol && avPoint2.Z() > -1.0 * bpZTol ){

        for ( Size_t iBP = 0; iBP < fBellyPlateCentrePoints.size(); iBP++ ){

          // Get the centre locations for each belly plate.
          TVector3 bpPos = fBellyPlateCentrePoints[ iBP ];

          // Find a vector to the position of the top and bottom of the belly plate
          TVector3 bpTop = bpPos + ( fBellyPlateHalfHeight * zAxis );
          bpTop.SetMag( bpPos.Mag() );
          TVector3 bpBottom = bpPos - ( fBellyPlateHalfHeight * zAxis );
          bpBottom.SetMag( bpPos.Mag() );

          // Then find the angle between the vectors to the top and the bottom
          // of the belly plate and divide by two to find the bisector.
          Double_t bpPhi = bpTop.Angle( bpBottom );
          Double_t startPhi = bpPos.Phi() - ( bpPhi / 2.0 );

          // Next define two vectors left and right of the centre point
          // of the belly plate.
          TVector3 bpLeft = bpPos;
          bpLeft.SetPhi( startPhi - phiTol );
          TVector3 bpRight = bpPos;
          bpRight.SetPhi( startPhi + bpPhi + phiTol );

          TVector3 avPoint2Cpy = avPoint2;
          avPoint2Cpy.SetZ( bpPos.Z() );

          // Now check to see if the angle of the AV intersection point
          // is within the angular region of the belly plate centre and
          // the left and right side of the plate itself.
          if ( avPoint2Cpy.Angle( bpPos ) < bpLeft.Angle( bpPos )
               && avPoint2Cpy.Angle( bpPos ) < bpRight.Angle( bpPos ) ){ return true; }

        }
      }
      if ( bellyPlateShadow == true ){ return true; }
      else{ return false; }
    }
  }

  // If the light path type is 'Water' only then it does not touch the AV, and therefore
  // cannot be shadowed due to belly plates.
  if ( lightPathType == "Water" ){ return false; }

  // Warn Message, Light Path Uninitialised or something has gone wrong.
  else{
    warn << "ShadowingCalculator::CheckForBellyPlateShadowing: Error! Cannot calculate shadowing for this light light path type: ";
    warn << lightPathType << "\n";
  }


  return false;
}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForNCDAnchorShadowing( RAT::DU::LightPathCalculator& lPath ){

  // Initialise values for the point on the AV to check for shadowing.
  TVector3 avPoint( 0.0, 0.0, 0.0 );

  // For paths going OUT OF the AV, we only need to check one point on the AV.
  // If the path is going accross the detector, there will be a second AV intersection point.
  TVector3 avPoint2( 0.0, 0.0, 0.0 );

  // Determine the light path type in order to identify which intersection point with the
  // AV to check with. For some light path types there may be more than one
  // crossing of the AV, and so we will need to perform two checks in these cases.

  // Notice that in the below, a light path type which just goes through the water
  // will not intersect with the AV at all, and hence doesn't need checking.

  string lightPathType = lPath.GetLightPathType();

  // One check required.
  if ( lightPathType == "Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV2nd();
  }

  // One check required.
  else if ( lightPathType == "AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV3rd();
  }

  // One check required.
  else if ( lightPathType == "AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV4th();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV2nd();
  }

  // Two checks required.
  else if ( lightPathType == "Water->Reflection->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV1st();
  }

  // Water type, no checks required.
  else if ( lightPathType == "Water" ){
    return false;
  }

  // Uninitialised light path - throw error.
  else if ( lightPathType == "Light Path Uninitialised" ){
    warn << "ShadowingCalculator::CheckForSupportRopeShadowing: LightPath is uninitialised, call CalcByPosition first before checking for shadowing!\n";
    return false;
  }

  // Set a tolerance equal to the NCD Anchor half-height and a little extra
  // (i.e. fNCDAnchorHalfHeight + 150.0 [=default fNCDAnchorTolerance mm) about the centre of the NCD
  const Double_t ncdTol = fNCDAnchorHalfHeight + fNCDAnchorTolerance;

  for ( Size_t iNCD = 0; iNCD < fNCDAnchorPoints.size(); iNCD++ ){
    if ( ( avPoint - fNCDAnchorPoints[ iNCD ] ).Mag() < ncdTol ){ return true; }

  }

  // If we reach here, then no shadowing was found for the first AV point.
  // If the light path type only involves a single intersection point
  // then we are done and can return the (should be 'false') shadowing bool.
  if ( ( lightPathType == "Scint/InnerAV->AV->Water"
         || lightPathType == "AV->Scint/InnerAV->AV->Water"
         || lightPathType == "AV->Water" ) ){ return false; }

  // Otherwise, for other light path types, we need to check the second AV intersection points.
  // As noted above though, we only reach this part if we haven't already encountered
  // shadowing in the first instances of the below path types.
  if ( ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water"
         || lightPathType == "Water->AV->Water"
         || lightPathType == "Water->Reflection->Water" ) ){

    for ( Size_t iNCD = 0; iNCD < fNCDAnchorPoints.size(); iNCD++ ){

      if ( ( avPoint2 - fNCDAnchorPoints[ iNCD ] ).Mag() < ncdTol ){ return true; }

    }
  }

  // if we reach here then there is no shadowing, so we just return 'false'
  return false;

}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForAVPipeShadowing( RAT::DU::LightPathCalculator& lPath ){

  // Initialise values for the point on the AV to check for shadowing.
  TVector3 avPoint( 0.0, 0.0, 0.0 );

  // For paths going OUT OF the AV, we only need to check one point on the AV.
  // If the path is going accross the detector, there will be a
  // second AV intersection point.
  TVector3 avPoint2( 0.0, 0.0, 0.0 );

  // Determine the light path type in order to identify which intersection point with the
  // AV to check with. For some light path types there may be more than one
  // crossing of the AV, and so we will need to perform two checks in these cases.

  string lightPathType = lPath.GetLightPathType();

  // One check required.
  if ( lightPathType == "Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV2nd();
  }

  // One check required.
  else if ( lightPathType == "AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV3rd();
  }

  // One check required.
  else if ( lightPathType == "AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV4th();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV2nd();
  }

  // These light path types will not intersect with the neck boss region
  else if ( lightPathType == "Water" ){ return false; }

  // Two check required.
  else if ( lightPathType == "Water->Reflection->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV1st();
  }

  // Light uninitialised, throw error.
  else if ( lightPathType == "Light Path Uninitialised" ){
    warn << "ShadowingCalculator::CheckForAVHDRopeShadowing: LightPath is uninitialised, call CalcByPosition first before checking for shadowing!\n";
    return false;
  }

  // Define the (theta, phi) coordinates in degrees of the
  // points on the AV where the light path intersects the AV
  Double_t avPointTheta = avPoint.Theta() * ( 180.0 / pi );
  Double_t avPointPhi = avPoint.Phi() * ( 180.0 / pi );

  // Check the positioning of the pipes in the (theta, phi) coordinates
  if ( avPointPhi < 75.0 || avPointPhi > 105.0 ){ return false; }
  if ( avPointTheta < 30.0
       && ( avPointPhi < 75.0 || avPointPhi > 105.0 ) ){ return false; }
  if ( avPointTheta >= 30
       && avPointTheta < 120.0
       && ( avPointPhi < 80.0 || avPointPhi > 100.0 ) ){ return false; }
  if ( avPointTheta >= 120.0
       && ( avPointPhi < 85.0 || avPointPhi > 95.0 ) ){ return false; }

  // If the light path is one of the below types, then there is only
  // one intersection point with the AV, therefore return false
  else if ( lightPathType == "Scint/InnerAV->AV->Water"
            || lightPathType == "AV->Scint/InnerAV->AV->Water"
            || lightPathType == "AV->Water"
            || lightPathType == "Water" ){ return true; }

  // otherwise we have a second point on the AV to consider, so check for
  // shadowing and return the appropriate value
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water"
            || lightPathType == "Water->AV->Water"
            || lightPathType == "Water->Reflection->Water" ){

    Double_t avPoint2Theta = avPoint.Theta() * ( 180.0 / pi );
    Double_t avPoint2Phi = avPoint.Phi() * ( 180.0 / pi );

    // Check the positioning of the pipes in the (theta, phi) coordinates
    // in relation to the second AV intersection point
    if ( avPoint2Phi < 75.0 || avPoint2Phi > 105.0 ){ return false; }
    if ( avPoint2Theta < 30.0
         && ( avPoint2Phi < 75.0 || avPoint2Phi > 105.0 ) ){ return false; }
    if ( avPoint2Theta >= 30
         && avPoint2Theta < 120.0
         && ( avPoint2Phi < 80.0 || avPoint2Phi > 100.0 ) ){ return false; }
    if ( avPoint2Theta >= 120.0
         && ( avPoint2Phi < 85.0 || avPoint2Phi > 95.0 ) ){ return false; }
    else { return true; }

  }

  return false;

}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForNeckBossShadowing( RAT::DU::LightPathCalculator& lPath )
{

  // Initialise values for the point on the AV to check for shadowing
  // and the starting direction from which to propagate the path from.
  TVector3 avPoint( 0.0, 0.0, 0.0 );

  // For paths going OUT OF the AV, we only need to check one point on the AV
  // and one direction. If the path is going accross the detector, there will be a
  // second AV intersection point.
  TVector3 avPoint2( 0.0, 0.0, 0.0 );

  // Determine the light path type in order to identify which intersection point with the
  // AV to check with. For some light path types there may be more than one
  // crossing of the AV, and so we will need to perform two checks in these cases.

  string lightPathType = lPath.GetLightPathType();

  // One check required.
  if ( lightPathType == "Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV2nd();
  }

  // One check required.
  else if ( lightPathType == "AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV3rd();
  }

  // One check required.
  else if ( lightPathType == "AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV4th();
  }

  // Two checks required.
  else if ( lightPathType == "Water->AV->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV2nd();
  }

  // These light path types will not intersect with the neck boss region
  else if ( lightPathType == "Water" ){ return false; }

  // Two check required.
  else if ( lightPathType == "Water->Reflection->Water" ){
    avPoint = lPath.GetPointOnAV1st();
    avPoint2 = lPath.GetPointOnAV1st();
  }

  // Light uninitialised, throw error.
  else if ( lightPathType == "Light Path Uninitialised" ){
    warn << "ShadowingCalculator::CheckForAVHDRopeShadowing: LightPath is uninitialised, call CalcByPosition first before checking for shadowing!\n";
    return false;
  }

  // If the point on the AV is above the theta position of
  // the neck boss, then consider it shadowed
  if ( avPoint.Theta() < fNeckBossTheta ){ return true; }


  // If the light path is one of the below types, then there is only
  // one intersection point with the AV, therefore return false
  else if ( lightPathType == "Scint/InnerAV->AV->Water"
            || lightPathType == "AV->Scint/InnerAV->AV->Water"
            || lightPathType == "AV->Water"
            || lightPathType == "Water" ){ return false; }

  // otherwise we have a second point on the AV to consider, so check for
  // shadowing and return the appropriate value
  else if ( lightPathType == "Water->AV->Scint/InnerAV->AV->Water"
            || lightPathType == "Water->AV->Water"
            || lightPathType == "Water->Reflection->Water" ){
    if ( avPoint2.Theta() < fNeckBossTheta ){ return true; }
    else { return false; }
  }

  // Else return false
  else{ return false; }

}

//////////////////////////////////////
//////////////////////////////////////

Bool_t
ShadowingCalculator::CheckForShadowing( RAT::DU::LightPathCalculator& lPath )
{

  // Initialise the shadowing boolean flags
  fBellyPlateShadowFlag = false;        // Shadowing from the belly plates
  fAVHDRopeShadowFlag = false;          // Shadowing from the AV hold-down ropes
  fNCDAnchorShadowFlag = false;         // Shadowing from NCD anchors
  fAVPipeShadowFlag = false;            // Shadowing from the AV Pipe Region
  fSupportRopeShadowFlag = false;       // Shadowing from the support ropes
  fNeckBossShadowFlag = false;          // Shadowing from the neck boss

  // Overall shadowing status. If any of the above are true, then this
  // is true. If all of the above remain false after shadowing checks,
  // this is false.
  fShadowFlag = false;

  if ( fExistAVHDRope && CheckForAVHDRopeShadowing( lPath ) == true ){
    fShadowFlag = true;
    fAVHDRopeShadowFlag = true;
    return true;
  }

  // Check for support rope shadowing
  else if ( fExistSupportRope && CheckForSupportRopeShadowing( lPath ) == true ){
    fShadowFlag = true;
    fSupportRopeShadowFlag = true;
    return true;
  }

  // Check for belly plate shadowing
  else if ( fExistBellyPlate && CheckForBellyPlateShadowing( lPath ) == true ){
    fShadowFlag = true;
    fBellyPlateShadowFlag = true;
    return true;
  }

  // Check for NCD anchor shadowing
  else if ( fExistNCDAnchor && CheckForNCDAnchorShadowing( lPath ) == true ){
    fShadowFlag = true;
    fNCDAnchorShadowFlag = true;
    return true;
  }

  // Check for AV pipe shadowing
  else if ( fExistAVPipe && CheckForAVPipeShadowing( lPath ) == true ){
    fShadowFlag = true;
    fAVPipeShadowFlag = true;
    return true;
  }

  // Check for neck boss shadowing
  else if ( fExistNeckBoss && CheckForNeckBossShadowing( lPath ) == true ){
    fShadowFlag = true;
    fNeckBossShadowFlag = true;
    return true;
  }

  // Otherwise no shadowing has been detected if we get this far
  // so we can set the overall shadowing boolean to false
  else { fShadowFlag = false; return false; }

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::SetAllGeometryTolerances( const Double_t geoTol )
{

  // Set all the tolerances for the...
  SetAVHDRopeTolerance( geoTol );    // AV hold-down ropes
  SetSupportRopeTolerance( geoTol ); // AV support (hold-up) Ropes
  SetBellyPlateTolerance( geoTol );  // Belly plates
  SetNCDAnchorTolerance( geoTol );   // NCD anchors
  SetAVPipeTolerance( geoTol );      // AV pipes
  SetNeckBossTolerance( geoTol );    // Neck boss

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::LoadAVHDRopeInfo(){

  // The AV Hold-Down rope basket is a strange geometrical creature...
  // ...so we have to be careful about what points are where in the geometry.

  // The basket itself is formed by a series of 5 'slings', which all start at the
  // bottom of the cavity floor, move up to the equator of the AV, loop over the
  // top of the AV and back down to the cavity floor.

  // A sling itself can be described as a pair of two
  // U-shaped ropes, near parallel to one another which are connected together
  // at two locations, each by an inter-connecting piece of rope. These 'inter-connection'
  // points are located at the top of the AV.

  // The coordinates for a single U-shaped rope in each sling are split up into
  // 6 (x,y,z) coordinates:
  // Cavity Floor -> AV Equator -> Top-AV -> Top-AV -> AV Equator -> Cavity Floor.

  // It is at both these 'Top-AV' locations that the two U-shaped ropes are inter-connected
  // forming a 'sling'

  // The defining orientation of each sling is based on its respective 'roll', i.e.
  // the angular displacement about the z-axis (points directly upwards through the AV neck)
  // In spherical polar coordinates, this could also classically be seen as the azimuthal angle,
  // where (0 < azimuthal angle < 2*Pi)

  // Because of this, we can't simply cut on the (X,Y,Z) positions, but we can perform a cylindrical
  // cut as we do for the support (hold-up) ropes for the sling's first, second, fifth and sixth
  // coordinates (i.e. equator and below). This is where the rope is below the equator. Below the equator of the AV the AVHD
  // ropes display cylindrical symmetry. Above the equator however we have to build up a series
  // of points which characterise the distinctive rope basket and use brute force to check against
  // each point to check for shadowing.

  DB *db = DB::Get();
  // Obtain the radius of the AV hold down rope
  db->GetLink( "GEO", "hold_down_ropes" )->GetS( "solid_definition" );
  fAVHDRopeRadius = db->GetLink( "SOLID",
                                 db->GetLink( "GEO", "hold_down_ropes" )->GetS( "solid_definition" ) )->GetD( "r_rope" );

  // Obtain the positions of an individual sling
  FillXYZPositions( fAVHDRopeSlingPoints, "SOLID",
                    db->GetLink( "GEO", "hold_down_ropes" )->GetS( "solid_definition" ),
                    "rope_x", "rope_y", "rope_z" );

  // Obtain all the 'roll' orientations of the 10 U-shaped ropes in degrees
  fAVHDRopeSlingRolls = db->GetLink( "GEO_LOCATIONS",
                                     db->GetLink( "GEO", "hold_down_ropes" )->GetS( "locations" ) )->GetDArray( "roll" );

  // Convert the degrees to radians
  for ( Size_t iRoll = 0; iRoll < fAVHDRopeSlingRolls.size(); iRoll++ ){
    fAVHDRopeSlingRolls[ iRoll ] /= ( 180.0 / TMath::Pi() );
  }

  // Take the 2nd and 5th coordinates of the u-shaped rope (where it passes the equator).
  TVector3 avEquator1st = fAVHDRopeSlingPoints[ 1 ];
  TVector3 avEquator2nd = fAVHDRopeSlingPoints[ 4 ];

  // Take the 1st and 6th coordinates of the u-shaped rope (where it is anchored to the cavity floor )
  TVector3 ropeAnchor1st = fAVHDRopeSlingPoints[ 0 ];
  TVector3 ropeAnchor2nd = fAVHDRopeSlingPoints[ 5 ];

  // For each u-shaped rope, locate the points on the AV equator where they hang down towards
  // the cavity floor. As each sling does this twice, we expect there to then be 4 * 5 = 20 locations
  // where the slings meet the AV equator.
  for ( Size_t iSling = 0; iSling < fAVHDRopeSlingRolls.size(); iSling++ ){

    TVector3 tmp1st = avEquator1st;
    TVector3 tmp2nd = avEquator2nd;

    tmp1st.RotateZ( fAVHDRopeSlingRolls[ iSling ]  );
    tmp2nd.RotateZ( fAVHDRopeSlingRolls[ iSling ]  );

    fAVHDRopeEquatorPoints.push_back( tmp1st );
    fAVHDRopeEquatorPoints.push_back( tmp2nd );

  }

  // Next, for each U-shaped rope, and it's own implicit set of six coordinates,
  // rotate each point about the z-axis so that each U-shaped rope and it's respective six points
  // are all located where they would be in the goemetry.

  for ( Size_t iRoll = 0; iRoll < fAVHDRopeSlingRolls.size(); iRoll++ ){

    // The points of the U-shaped rope on the AV equator and upper AV.
    // rotated about the z-axis with their respective 'roll' value.
    TVector3 pointUpEq1 = fAVHDRopeSlingPoints[ 1 ];
    pointUpEq1.RotateZ( fAVHDRopeSlingRolls[ iRoll ] );
    TVector3 pointUpEq2 = fAVHDRopeSlingPoints[ 4 ];
    pointUpEq2.RotateZ( fAVHDRopeSlingRolls[ iRoll ] );
    TVector3 pointUpAV1 = fAVHDRopeSlingPoints[ 2 ];
    pointUpAV1.RotateZ( fAVHDRopeSlingRolls[ iRoll ] );
    TVector3 pointUpAV2 = fAVHDRopeSlingPoints[ 3 ];
    pointUpAV2.RotateZ( fAVHDRopeSlingRolls[ iRoll ] );

    // NOTE: Whilst each inter-connecting point on each rope
    // isn't directly opposite it's equivalent on the other side,
    // it is directly opposite the inter-connecting point on its
    // partner rope. So we can use this to identify the locations
    // on the AV of each rope's partnering rope. (below - rotation by pi.)

    // Now match up with the sling with its matching partner
    TVector3 pointUpEq2Pair = pointUpEq1;
    pointUpEq2Pair.RotateZ( TMath::Pi() );
    TVector3 pointUpEq1Pair = pointUpEq2;
    pointUpEq1Pair.RotateZ( TMath::Pi() );
    TVector3 pointUpAV2Pair = pointUpAV1;
    pointUpAV2Pair.RotateZ( TMath::Pi() );
    TVector3 pointUpAV1Pair = pointUpAV2;
    pointUpAV1Pair.RotateZ( TMath::Pi() );

    // So, keeping with the description, for each sling (each consisting of two U-shaped ropes)
    // we have four locations on the AV equator and four locations on the upper AV. Eight in total.

    // Put all these points into the (x,y) plane and set their magnitudes
    // to the outer radius of the AV.
    pointUpEq1.SetZ( 0.0 );
    pointUpEq1.SetMag( fAVOuterRadius );
    pointUpEq2.SetZ( 0.0 );
    pointUpEq2.SetMag( fAVOuterRadius );

    pointUpAV1.SetZ( 0.0 );
    pointUpAV2.SetZ( 0.0 );

    pointUpEq1Pair.SetZ( 0.0 );
    pointUpEq1Pair.SetMag( fAVOuterRadius );
    pointUpEq2Pair.SetZ( 0.0 );
    pointUpEq2Pair.SetMag( fAVOuterRadius );

    pointUpAV1Pair.SetZ( 0.0 );
    pointUpAV2Pair.SetZ( 0.0 );

    // Now, in the (x,y) plane we define direction vectors from the equator going inwards.
    // We then look to project these up to points onto the AV (i.e. extend in the z-direciton).

    // Outer Equatorial Rope parts on the upper AV
    TVector3 eqToAV1 = ( pointUpAV1 - pointUpEq1 );
    TVector3 eqToAV2 = ( pointUpAV2 - pointUpEq2 );
    TVector3 eqToAV1Pair = ( pointUpAV1Pair - pointUpEq1Pair );
    TVector3 eqToAV2Pair = ( pointUpAV2Pair - pointUpEq2Pair );

    Double_t eqToAV1Mag = ( eqToAV1.Mag() );
    Double_t eqToAV2Mag = ( eqToAV2.Mag() );

    // The above accounted for direction vectors in the (x,y) plane from the equator
    // in towards the first points on the upper AV (projected into the (x,y) plane).
    // Now however, we need to go further and account for the direction vectors in the
    // (x,y) plane to the neck region (in the same plane). Note that the ropes do not
    // cross directly over the centre of the neck, they wrap around the outside of it.
    // Because of this we need to perturb the crossing point of each rope at the top
    // of the AV by some offset away from the neck (below).

    // Inner Equatorial Rope parts on the upper AV

    // Start by a direct line across the two upper AV points, we want to 'kink' it
    // about the neck position.
    Double_t magAcr = ( pointUpEq2 - pointUpEq1 ).Mag();
    Double_t magAcrDiv2 = magAcr / 2.0;
    // The 'kink'.
    TVector3 deltaOrigin = pointUpEq1 + ( magAcrDiv2 * ( pointUpEq2 - pointUpEq1 ).Unit() );
    deltaOrigin.SetMag( 250.0 + fNeckOuterRadius );
    // And the same 'kink' on the other side for the ropes' partner.
    TVector3 deltaOriginPair = deltaOrigin;
    deltaOriginPair.RotateZ( TMath::Pi() );
    // NOTE: May need to change to this

    // Now define the direction vectors from the upper equator points
    // (of which there are four per sling) to the top of the AV (in the (x,y) plane).
    TVector3 av1ToTop = ( deltaOrigin - pointUpAV1 );
    TVector3 av2ToTop = ( deltaOrigin - pointUpAV2 );
    TVector3 av1ToTopPair = ( deltaOriginPair - pointUpAV1Pair );
    TVector3 av2ToTopPair = ( deltaOriginPair - pointUpAV2Pair );

    Double_t av1ToTopMag = av1ToTop.Mag();
    Double_t av2ToTopMag = av2ToTop.Mag();

    // Finally, setup the connecing sling pieces.
    TVector3 slingTmpPos1 = pointUpAV1;
    slingTmpPos1.SetZ( 0.0 );

    TVector3 slingTmpPos2 = pointUpAV2;
    slingTmpPos2.SetZ( 0.0 );

    // Now we will scan along the direction vectors in the (x,y) plane
    // (from equator to equator) and project upwards to extend the z-component
    Double_t maxIter = 1000.0;
    TVector3 zAxis( 0.0, 0.0, 1.0 );

    Double_t angleReg = pointUpAV1.Angle( pointUpAV1Pair );

    for ( Double_t iSeg = 0.1; iSeg < maxIter; iSeg = iSeg + 2.0 ){

      //////////////////////////////////////////////////////////////////////////
      // Scan between the interconnecting area and populate with points first.
      //////////////////////////////////////////////////////////////////////////
      slingTmpPos1.SetPhi( ( iSeg / maxIter ) * ( 2.0 * TMath::Pi() ) );
      if ( slingTmpPos1.Angle( pointUpAV1 ) < angleReg
           && slingTmpPos1.Angle( pointUpAV1Pair ) < angleReg ){

        TVector3 pointOnConnect = VectorToSphereEdge( slingTmpPos1, zAxis, fAVOuterRadius, 0.0 );
        fAVHDRopeUpperAVPoints.push_back( pointOnConnect );

      }

      slingTmpPos2.SetPhi( ( iSeg / maxIter ) * ( 2.0 * TMath::Pi() ) );
      if ( slingTmpPos2.Angle( pointUpAV2 ) < angleReg
           && slingTmpPos2.Angle( pointUpAV2Pair ) < angleReg ){

        TVector3 pointOnConnect = VectorToSphereEdge( slingTmpPos2, zAxis, fAVOuterRadius, 0.0 );
        fAVHDRopeUpperAVPoints.push_back( pointOnConnect );

      }
      //////////////////////////////////////////////////////////////////////////

      //////////////////////////////////////////////////////////////////////////
      // Now populate the rest of the u-shaped rope components for each pair, forming the basket.
      //////////////////////////////////////////////////////////////////////////
      TVector3 startPos1 = pointUpEq1 + ( ( iSeg / maxIter ) * eqToAV1Mag * eqToAV1.Unit() );
      TVector3 point1 = VectorToSphereEdge( startPos1, zAxis, fAVOuterRadius, 0.0 );
      TVector3 startPos2 = pointUpEq2 + ( ( iSeg / maxIter ) * eqToAV2Mag * eqToAV2.Unit() );
      TVector3 point2 = VectorToSphereEdge( startPos2, zAxis, fAVOuterRadius, 0.0 );

      TVector3 avStartPos1 = pointUpAV1 + ( ( iSeg / maxIter ) * av1ToTopMag * av1ToTop.Unit() );
      TVector3 avPoint1 = VectorToSphereEdge( avStartPos1, zAxis, fAVOuterRadius, 0.0 );
      TVector3 avStartPos2 = pointUpAV2 + ( ( iSeg / maxIter ) * av2ToTopMag * av2ToTop.Unit() );
      TVector3 avPoint2 = VectorToSphereEdge( avStartPos2, zAxis, fAVOuterRadius, 0.0 );

      fAVHDRopeUpperAVPoints.push_back( point1 );
      fAVHDRopeUpperAVPoints.push_back( point2 );

      fAVHDRopeUpperAVPoints.push_back( avPoint1 );
      fAVHDRopeUpperAVPoints.push_back( avPoint2 );
      //////////////////////////////////////////////////////////////////////////

    }
  }


  // For each sling, locate the points below the AV equator where the
  // slings become anchored to the cavity floor

  for ( Size_t iSling = 0; iSling < fAVHDRopeSlingRolls.size(); iSling++ ){

    // Locate the point on the equator and find the direction down towards the cavity floor.
    // using this vector, populate with points.
    TVector3 anchor1st = ropeAnchor1st;
    TVector3 equator1st = avEquator1st;

    TVector3 anchor2nd = ropeAnchor2nd;
    TVector3 equator2nd = avEquator2nd;

    // Rotate for each of the u-shaped rope pieces.
    anchor1st.RotateZ( fAVHDRopeSlingRolls[ iSling ]  );
    equator1st.RotateZ( fAVHDRopeSlingRolls[ iSling ]  );

    anchor2nd.RotateZ( fAVHDRopeSlingRolls[ iSling ]  );
    equator2nd.RotateZ( fAVHDRopeSlingRolls[ iSling ]  );

    // Define the direction from the equator to the cavity floor.
    Float_t eqToAnchorDist1st = ( anchor1st - equator1st ).Mag();
    Float_t eqToAnchorDist2nd = ( anchor2nd - equator2nd ).Mag();

    // Divide the region into 100 points to iterate and populate with points.
    eqToAnchorDist1st /= 100.0;
    eqToAnchorDist2nd /= 100.0;

    TVector3 eqToAnchor1st = ( anchor1st - equator1st ).Unit();
    TVector3 eqToAnchor2nd = ( anchor2nd - equator2nd ).Unit();

    // Fill in the points.
    for ( Int_t nSections = 1; nSections < 100.0; nSections++ ){

      TVector3 pointBelowAV1st = equator1st + ( nSections * eqToAnchorDist1st * eqToAnchor1st );
      TVector3 pointBelowAV2nd = equator2nd + ( nSections * eqToAnchorDist2nd * eqToAnchor2nd );

      fAVHDRopeLowerAVPoints.push_back( pointBelowAV1st );
      fAVHDRopeLowerAVPoints.push_back( pointBelowAV2nd );

    }

  }

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::LoadSupportRopeInfo(){

  DB *db = DB::Get();

  // Obtain the radius of the support rope object and it's height in
  // it's locale frame (i.e. it's own length)
  fSupportRopeRadius = db->GetLink( "SOLID",
                                    db->GetLink( "GEO", "hold_up_ropes" )->GetS( "solid_definition" ) )->GetD( "r_rope" );
  fSupportRopeHeight = db->GetLink( "SOLID",
                                    db->GetLink( "GEO", "hold_up_ropes" )->GetS( "solid_definition" ) )->GetD( "z_rope" );

  Double_t bellyPlateSplit =  db->GetLink( "SOLID", "belly_plate_snoman" )->GetD( "octahedron_half_height" );

  std::vector< TVector3 > equatorPoints;
  // Obtain the (X,Y,Z) coordinates of all the support rope locations along the equator in the detector geometry
  FillXYZPositions( equatorPoints, "GEO_LOCATIONS",
                    db->GetLink( "GEO", "hold_up_ropes" )->GetS( "locations" ) );

  // In the following 'for' loop, this is what's going on:
  // The central location for the support rope u-shape loop is inbetween
  // the two actual ropes locations themselves. The ropes actually are slung
  // through either side of this central point. From the centre outwards they
  // are about 2/3s of the belly plate half-width(height) across. Therefore the
  // rope locations themselves are (2/3)*(half-width) eitherside of the
  // belly plate centre (which is the centre of the support rope u-shape.

  for ( Size_t iVal = 0; iVal < equatorPoints.size(); iVal++ ){

    // Define a z-unit vector and two vectors to the equator.
    TVector3 zAxis( 0.0, 0.0, 1.0 );
    TVector3 pointOnEquator1 = equatorPoints[ iVal ];
    TVector3 pointOnEquator2 = pointOnEquator1;

    // Displace both vectors (2/3)*(half-widths) eitherside of the centre
    Double_t angDispPhi = ( ( 2.0 / 3.0 ) * bellyPlateSplit ) / ( fAVOuterRadius );
    pointOnEquator1.SetPhi( pointOnEquator1.Phi() + angDispPhi );
    pointOnEquator2.SetPhi( pointOnEquator2.Phi() - angDispPhi );

    // Put these into the equator point storage vector
    fSupportRopeEquatorPoints.push_back( pointOnEquator1 );
    fSupportRopeEquatorPoints.push_back( pointOnEquator2 );

    // Now iterate over 100 points in the z-direction up each of these
    // equator locations to trace out the geometry of the support ropes themselves
    Double_t height = fSupportRopeHeight;
    height /= 100.0;
    for ( Int_t kVal = 1; kVal < 100; kVal++ ){
      TVector3 pointAboveEquator1 = pointOnEquator1 + ( kVal * height * zAxis );
      TVector3 pointAboveEquator2 = pointOnEquator2 + ( kVal * height * zAxis );
      fSupportRopeUpperAVPoints.push_back( pointAboveEquator1 );
      fSupportRopeUpperAVPoints.push_back( pointAboveEquator2 );
    }

  }

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::LoadBellyPlateInfo(){

    DB *db = DB::Get();

  // Obtain the half-height of the belly plate object. We use 'belly_plate_outer'
  // as we are only concerned with the extremities of the belly plate geometry.
  fBellyPlateHalfHeight = db->GetLink( "SOLID", "belly_plate_snoman" )->GetD( "octahedron_half_height" );

  // Obtain the (X,Y,Z) coordinates of all the belly plate locations in the detector geometry
  FillXYZPositions( fBellyPlateCentrePoints, "GEO_LOCATIONS",
                    db->GetLink( "GEO", "inner_av" )->GetS( "belly_plate_locations" ));

  Double_t tol = fBellyPlateTolerance;
  Double_t ztol = fBellyPlateTolerance;
  TVector3 zAxis( 0.0, 0.0, 1.0 );

  for ( Size_t iBP = 0; iBP < fBellyPlateCentrePoints.size(); iBP++ ){

    TVector3 bpPos = fBellyPlateCentrePoints[ iBP ];
    TVector3 bpTop = bpPos + ( fBellyPlateHalfHeight * zAxis );
    bpTop.SetMag( bpPos.Mag() );
    TVector3 bpBottom = bpPos - ( fBellyPlateHalfHeight * zAxis );
    bpBottom.SetMag( bpPos.Mag() );
    Double_t angleBP = bpBottom.Angle( bpTop );
    Double_t startPhi = bpPos.Phi() - ( angleBP / 2.0 );

    for ( Int_t iTol = 0.0; iTol < tol; iTol = iTol + 1.0 ){

      TVector3 tmpPoint = bpBottom;
      tmpPoint.SetPhi( startPhi + ( ( iTol / tol ) * angleBP ) );

      for ( Int_t iZ = 0.0; iZ < ztol; iZ = iZ + 1.0 ){

        TVector3 bpPoint = tmpPoint + ( 2 * ( iZ / ztol ) * fBellyPlateHalfHeight * zAxis );
        bpPoint.SetMag( bpPos.Mag() );
        fBellyPlatePoints.push_back( bpPoint );

      }
    }
  }

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::LoadNCDAnchorInfo(){

  DB *db = DB::Get();

  // Obtain the half-height and radius of the NCD anchor object
  fNCDAnchorHalfHeight = db->GetLink( "SOLID",
                                      db->GetLink( "GEO", "inner_av" )->GetS( "ncd_anchor_definition" ) )->GetD( "half_z" );
  fNCDAnchorRadius = db->GetLink( "SOLID",
                                  db->GetLink( "GEO", "inner_av" )->GetS( "ncd_anchor_definition" ) )->GetD( "r_max" );

  // Obtain the (X,Y,Z) coordinates of all the NCD anchor locations in the detector geometry
  FillXYZPositions( fNCDAnchorPoints, "GEO_LOCATIONS",
                    db->GetLink( "GEO", "inner_av" )->GetS( "ncd_anchor_locations" ) );

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::LoadAVPipeInfo(){

  DB *db = DB::Get();

  // Obtain the two max radii of the two types of AV pipe
  fAVPipeXRadius = db->GetLink( "SOLID", "av_pipe-1" )->GetD( "r_max" );
  fAVPipeLabXRadius = db->GetLink( "SOLID", "av_pipe-labs" )->GetD( "r_max" );

  // Obtain the (X,Y,Z) coordinates of all seven pipes
  FillXYZPositions( fAVPipe1Points, "SOLID", "av_pipe-1" );
  FillXYZPositions( fAVPipe2Points, "SOLID", "av_pipe-2" );
  FillXYZPositions( fAVPipe3Points, "SOLID", "av_pipe-3" );
  FillXYZPositions( fAVPipe4Points, "SOLID", "av_pipe-4" );
  FillXYZPositions( fAVPipeLabSPoints, "SOLID", "av_pipe-labs" );
  FillXYZPositions( fAVPipe6Points, "SOLID", "av_pipe-6" );
  FillXYZPositions( fAVPipeLabRPoints, "SOLID", "av_pipe-labr" );

  // Create a vector which holds all the vector of points defined above
  std::vector< std::vector< TVector3 > > pipeVectors;
  pipeVectors.push_back( fAVPipe1Points );
  pipeVectors.push_back( fAVPipe2Points );
  pipeVectors.push_back( fAVPipe3Points );
  pipeVectors.push_back( fAVPipe4Points );
  pipeVectors.push_back( fAVPipeLabSPoints );
  pipeVectors.push_back( fAVPipe6Points );
  pipeVectors.push_back( fAVPipeLabRPoints );

  // Iterate of each pipe vector, and trace out the points to store
  // in fAVPipePoints through interpolation
  for ( Size_t iPipe = 0; iPipe < pipeVectors.size(); iPipe++ ){

    for ( Size_t iP = 0; iP < pipeVectors[ iPipe ].size() - 1; iP++ ){

      TVector3 p1 = pipeVectors[ iPipe ][ iP ];
      TVector3 p2 = pipeVectors[ iPipe ][ iP + 1 ];
      TVector3 p1Top2 = ( p2 - p1 );

      Double_t kVal = 100.0;
      for ( Double_t k = 0.0; k <= kVal; k += 2.0 ){

        TVector3 pipePoint = p1 + ( ( ( k / kVal ) * p1Top2.Mag() ) * ( p1Top2.Unit() ) );
        fAVPipePoints.push_back( pipePoint );

      }

    }

  }

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::LoadNeckBossInfo()
{

  DB *db = DB::Get();

  // Obtain the z-position for the neck boss, located at the centre of the
  // union of two rings of different heights and radii which form the neckboss
  Double_t nBossZ = db->GetLink( "GEO", "av" )->GetD( "neck_boss_z" );

  // Obtain the half-heights of the aforementioned rings
  Double_t nBossUpperBevelHalfZ = db->GetLink( "SOLID", "neck_boss" )->GetD( "half_height1" );
  Double_t nBossLowerBevelHalfZ = db->GetLink( "SOLID", "neck_boss" )->GetD( "half_height2" );

  // Locate the bottom of the neck boss in the z-coordinates
  Double_t nBossBottomZ = nBossZ - ( ( 2.0 * ( nBossUpperBevelHalfZ + nBossLowerBevelHalfZ ) ) / 2.0 );

  // Now project out from the bottom z location to the edge of the AV.
  TVector3 nBossTmp( 0.0, 0.0, nBossBottomZ );
  TVector3 unitX( 1.0, 0.0, 0.0 );
  TVector3 nBossAVEdge = VectorToSphereEdge( nBossTmp, unitX, fAVInnerRadius, 0 );

  // Use the theta value of this point on the AV as a maximum theta value
  // for which all points with smaller theta values will be shadowed
  fNeckBossTheta = nBossAVEdge.Theta();

}

//////////////////////////////////////
//////////////////////////////////////

void
ShadowingCalculator::FillXYZPositions( std::vector< TVector3 >& positions,
                                       const std::string& tableName,
                                       const std::string& indexName,
                                       const std::string& xField,
                                       const std::string& yField,
                                       const std::string& zField )
{

  DB *db = DB::Get();
  // Obtain the (X,Y,Z) coordinates the geometry represented by 'indexName' located in the 'tblName' table
  std::vector< Double_t > xPos = db->GetLink( tableName, indexName )->GetDArray( xField );
  std::vector< Double_t > yPos = db->GetLink( tableName, indexName )->GetDArray( yField );
  std::vector< Double_t > zPos = db->GetLink( tableName, indexName )->GetDArray( zField );

  // Check that there are an equal number of entries for the X,Y and Z coordinates
  if ( xPos.size() != yPos.size()
       || xPos.size() != zPos.size()
       || yPos.size() != zPos.size() ){
    warn << "ShadowingCalculator::FillXYZPositions: Inequal number of elements across the (X,Y,Z) coordinate arrays!\n";
    warn << "(X,Y,Z) positions not initialised for (" << tableName << ", " << indexName << ").\n";
    return;
  }

  // Now loop over all the locations, putting them into a single std::vector (fNCDAnchorPoints)
  // of TVector3 objects
  positions.clear();
  positions.reserve(xPos.size());
  for ( Int_t iP = 0; iP < (Int_t)xPos.size(); iP++ ){

    // Create the NCD point...
    TVector3 point( xPos[ iP ], yPos[ iP ], zPos[ iP ] );

    // ...and add it to the list of NCD anchor locations
    positions.push_back( point );

  }

}

//////////////////////////////////////
//////////////////////////////////////

TVector3
ShadowingCalculator::VectorToSphereEdge( const TVector3& startPos,
                                         const TVector3& startDir,
                                         const Double_t radiusFromCentre,
                                         const bool outside )
{

  // The a, b, and c coefficients of a typical quadratic equation
  const Double_t aCoeff = 1.0;
  const Double_t bCoeff = 2.0 * startPos.Mag() * TMath::Cos( startPos.Angle( startDir ) );
  const Double_t cCoeff = startPos.Mag2() - ( radiusFromCentre * radiusFromCentre );

  // Parameter (set to 0.0 to begin with)
  Double_t distParam = 0.0;

  // Discriminant for quadratic equation
  const Double_t discrim = TMath::Power( bCoeff, 2 ) - ( 4.0 * aCoeff * cCoeff );

  // Intersection Calculation
  if ( discrim >= 0.0 ){

    const Double_t discrimPlus = ( ( -1.0 * bCoeff ) + TMath::Sqrt(discrim) ) / ( 2.0 * aCoeff );
    const Double_t discrimMinus = ( ( -1.0 * bCoeff ) - TMath::Sqrt(discrim) ) / ( 2.0 * aCoeff );

    if ( discrimPlus > 0.0 && !outside ){ distParam = discrimPlus; }
    else { distParam = std::min( discrimMinus, discrimPlus ); }

  }

  else {
    debug << "ShadowingCalculator::VectorToSphereEdge: Erroneous calculation - negative discriminant" << endl;
    debug << "Discriminant Value: " << discrim << endl;
    debug << "ACoeff: " << aCoeff << ", BCoeff: " << bCoeff << ", CCoeff: " << cCoeff << endl;
    debug << "StartPos.Mag(): " << startPos.Mag() << " and R: " << radiusFromCentre << endl;
  }

  // Implementation of Parameterisation to find intersection point
  TVector3 endPointVec = startPos + ( distParam * startDir );

  return endPointVec;

}


//////////////////////////////////////
//////////////////////////////////////

TVector3
ShadowingCalculator::VectorToCylinderEdge( const TVector3& startPos,
                                           const TVector3& startDir,
                                           const TVector3& cylinderBaseOrigin,
                                           const Double_t cylinderRadius )
{

  TVector3 startDirXY = startDir;
  startDirXY.SetZ( 0.0 );

  TVector3 startPosXY = startPos;
  startPosXY.SetZ( 0.0 );

  TVector3 pointOnCylinderBase = VectorToSphereEdge( startPosXY, startDirXY, cylinderRadius, 0 );

  Double_t angleTmp = ( pointOnCylinderBase.Unit() ).Angle( startDir );
  Double_t heightInCylinder = 0.0;

  if ( angleTmp != halfpi ){
    heightInCylinder = ( ( pointOnCylinderBase - startPosXY ).Mag() ) * TMath::Tan( angleTmp );
  }

  TVector3 pointOnCylinder = pointOnCylinderBase + cylinderBaseOrigin;
  pointOnCylinder.SetZ( pointOnCylinder.Z() + heightInCylinder );

  return pointOnCylinder;

}

//////////////////////////////////////
//////////////////////////////////////

Double_t
ShadowingCalculator::ClosestAngle( const TVector3& eventPos,
                                   const Double_t edgeRadius )
{

  // Calculate one of the angles of the Right-Angled triangle to obtain
  // the other
  Double_t Phi = TMath::ACos( edgeRadius / ( eventPos.Mag() ) );
  Double_t angle = halfpi - Phi;

  return angle;

}

//////////////////////////////////////
//////////////////////////////////////