////////////////////////////////////////////////////////////////////////
///
/// \class
///
/// \brief
///
/// \author David Auty auty@ualberta.ca
///
/// REVISION HISTORY: Jie Hu jhu9@ualberta.ca - new strategy and modified
/// to fit J Tseng's MultiPathFunction
///
/// \detail
///       To fit data in a partial scint-water fill geometry
////
////////////////////////////////////////////////////////////////////////

#include <RAT/MultiPathScintWater.hh>
#include <CLHEP/Random/RandFlat.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <RAT/DB.hh>
#include <RAT/DU/Utility.hh>
#include "TFile.h"
#include "TH1F.h"
#include "TMath.h"
#include <TF1.h>
//#define DEBUG

using namespace RAT;
using namespace ROOT;

MultiPathScintWater::MultiPathScintWater() {
  // parameter definitions
  fPars.push_back(Parameter("x",1678,-8390.0,8390.0));
  fPars.push_back(Parameter("y",1678,-8390.0,8390.0));
  fPars.push_back(Parameter("z",1678,-8390.0,8390.0));
  fPars.push_back(Parameter("t",400,-100.0,300.0));
}

void MultiPathScintWater::BeginOfRun(DS::Run&) {
#ifdef MULTIPATH_FIXED_SEED
  //step = 0.0001;
  step = 0.5001;
#endif
  DB* db = DB::Get();
  DBLinkPtr dbLink = db->GetLink( "FIT_MULTIPATH" );
  fPDFDataCh = dbLink->GetDArray( "sPDF_multipathfit" );
  fPDFDataScint = dbLink->GetDArray( "sPDF_scint" );
  fScintRI = dbLink->GetD( "scint_RI" ); //for grVelocity
  fWaterRI = dbLink->GetD( "water_RI" ); // effective water refraction index; n = 1.33 for normal case
  double fChBinWidth = dbLink->GetD( "time_bin_width" ); // time resolution in ns, set to 0.25 ns
  double fChOffset = dbLink->GetD( "time_offset" );
  fMaxPMTtRes = dbLink->GetD( "max_PMT_tRes" );
  fMaxReflects = dbLink->GetI( "maximum_reflections");
  fBoundary = dbLink->GetD( "boundary_tolerance" );
  fNeckpathEnable = dbLink->GetZ( "neckpath_enable" );
  rpsup = dbLink->GetD( "psup_rad" );
  fPSUPRadius2 = rpsup * rpsup;

  fZoff = dbLink->GetD( "av_offset_z" ); //z offset of AV to PSUP center
  fAVRadius = db->GetLink( "SOLID", "acrylic_vessel_inner" )->GetD( "r_sphere" ); // 6005.0, AV inner radius
  fAVRadiusOuter = db->GetLink( "SOLID", "acrylic_vessel_outer" )->GetD( "r_sphere" ); // 6060, AV outer radius
  fNeckRadius = db->GetLink( "SOLID", "acrylic_vessel_inner" )->GetD( "r_neck" ); // inner radius of neck = 730 mm
  fNeckRadiusOuter = db->GetLink( "SOLID", "acrylic_vessel_outer" )->GetD( "r_neck" ); // outer radius of neck = 785 mm
  fZneckLo = sqrt( fAVRadiusOuter*fAVRadiusOuter - fNeckRadiusOuter*fNeckRadiusOuter ) + fZoff; // bottome of the neck

  double fSpeedOfLight = CLHEP::c_light;
  fSpeedOfLightWater = fSpeedOfLight / fWaterRI;
  fSpeedOfLightScint = fSpeedOfLight / fScintRI;

  fMaxNumStartPos = dbLink->GetI( "max_start_positions" );
  fRhoCut = dbLink->GetD( "rho_sampling_cut" ); // sampling only in AV, for future test
  //fGoodFitThreshold = dbLink->GetI( "low_DeltaLogL_threshold" );// for future test

  DBLinkPtr innerAvDb = db->GetLink( "GEO", "inner_av" );
  fWaterLevel = innerAvDb->GetD( "split_z" );

  fEntriesTimeCh = fPDFDataCh.size();// 1600
  double xtemp = fPDFDataCh[0];
  int j;
  for( j = fEntriesTimeCh - 1; j >= 0; j-- ) {
    if( fPDFDataCh[j] == fPDFDataCh[0] ) {
      fPDFDataCh[j] = xtemp;
      xtemp -= 0.001;//0.1 in FitterAW
    }
  }

  xtemp = fPDFDataCh[ fEntriesTimeCh-1 ];

  for( j = fEntriesTimeCh - 1; fPDFDataCh[j] == xtemp; j-- );

  for( j++; j < fEntriesTimeCh; j++ ) {
    fPDFDataCh[j] = xtemp;
    xtemp -= 0.001;//0.1 in FitterAW
  }

  for( j = 0; j < fEntriesTimeCh; j++ ) {
    fPDFxCh.push_back(fChBinWidth * static_cast<double>(j + 1) - fChOffset);
    fPDFCh.push_back(fPDFDataCh[j]);
  }

  for( j = 0; j < fEntriesTimeCh - 1; j++ )
    fDerivativeCh.push_back( ( fPDFCh[j + 1] - fPDFCh[j] ) / fChBinWidth );

  /// fill scint pdf
  double xtemp1 = fPDFDataScint[0];
  for( j = fEntriesTimeCh-1; j >= 0; j-- ) {
    if( fPDFDataScint[j] == fPDFDataScint[0] ) {
      fPDFDataScint[j] = xtemp1;
      xtemp1 -= 0.001;
    }
  }

  xtemp1 = fPDFDataScint[ fEntriesTimeCh-1 ];

  for( j = fEntriesTimeCh-1; fPDFDataScint[j] == xtemp1; j-- );

  for( j++; j < fEntriesTimeCh; j++ ) {
    fPDFDataScint[j] = xtemp1;
    xtemp1 -= 0.001;
  }

  for( j = 0; j < fEntriesTimeCh; j++ ) {
    fPDFxScint.push_back(fChBinWidth * static_cast<double>(j + 1) - fChOffset);
    fPDFScint.push_back(fPDFDataScint[j]);
  }

  for( j = 0; j < fEntriesTimeCh - 1; j++ )
    fDerivativeScint.push_back( ( fPDFScint[j + 1] - fPDFScint[j] ) / fChBinWidth );

  /// calculate cosTheta, reflection and transimission probabilities (Fresnel equations)
  double step_ct = 1./fMaxReflects;
  for ( double cosTi = 0; cosTi < 1.; cosTi += step_ct ) { // cosThetaIncident
    double Rp = 0., Rs = 0.; // reflection prob for p-wave and s-wave
    double sinTi = sqrt( 1 - cosTi*cosTi ); // sinThetaIncident square
    // incident in water, transmit in scint
    double cosTt = 1 - (fWaterRI/fScintRI*sinTi)*(fWaterRI/fScintRI*sinTi);
    if( cosTt < 0 ) { // total internal reflection
      Rp = 1;
      Rs = 1;
      fRWater.push_back( ( Rp + Rs )/2. );
    } else { // transmission (refraction) and reflection
      cosTt = sqrt( cosTt ); // cosThetaTransmit
      Rs = pow( ( fWaterRI*cosTi - fScintRI*cosTt )/( fWaterRI*cosTi + fScintRI*cosTt ), 2 );
      Rp = pow( ( fWaterRI*cosTt - fScintRI*cosTi )/( fWaterRI*cosTt + fScintRI*cosTi ), 2 );
      fRWater.push_back( ( Rp + Rs )/2. );
    }
    // incident in scint, transmit in water
    Rs = 0, Rp = 0, cosTt = 0;
    cosTt = 1 - (fScintRI/fWaterRI*sinTi)*(fScintRI/fWaterRI*sinTi);
    if( cosTt < 0 ) { // total internal reflection
      Rp = 1;
      Rs = 1;
      fRScint.push_back( ( Rp + Rs )/2. );
      //fTScint.push_back( ( 1-Rp + 1-Rs )/2. );
    } else {
    Rs = pow( ( fScintRI*cosTi - fWaterRI*cosTt )/( fScintRI*cosTi + fWaterRI*cosTt ), 2 );
    Rp = pow( ( fScintRI*cosTt - fWaterRI*cosTi )/( fScintRI*cosTt + fWaterRI*cosTi ), 2 );
    fRScint.push_back( ( Rp + Rs )/2. );
   }
  } // end of for loop
  fRScint.push_back( 0 );
  fRWater.push_back( 0 );

#ifdef DEBUG
  //Check the fitting pdfs and their derivatives for any event
  TFile f("pdfDerivatives_MultiPath.root","RECREATE");
  TH1F h1("h1","pdfCherenkov",fEntriesTimeCh, -100, 300);
  TH1F h1x("h1x","pdfXCh",fEntriesTimeCh, -100, 300);
  TH1F dh1("dh1","DerivativeCherenkov",fEntriesTimeCh,-100, 300);

  TH1F h2("h2","pdfScint",fEntriesTimeCh, -100, 300);
  TH1F h2x("h2x","pdfXScint",fEntriesTimeCh, -100, 300);
  TH1F dh2("dh2","DerivativeScint",fEntriesTimeCh, -100, 300);

  TH1F hRWater("hRWater","Unpolarized light, W-Scint", fMaxReflects, 0, 1);
  TH1F hRScint("hRScint","Unpolarized light, Scint-W", fMaxReflects, 0, 1);


  for( j = 0; j < fEntriesTimeCh; j++ ) {
    h1.SetBinContent(j+1, fPDFCh[j]);
    h1x.SetBinContent(j+1, fPDFxCh[j]);
    dh1.SetBinContent(j+1, fDerivativeCh[j]);
  }
  for( j = 0; j < fEntriesTimeCh; j++ ) {
    h2.SetBinContent(j+1, fPDFScint[j]);
    h2x.SetBinContent(j+1, fPDFxScint[j]);
    dh2.SetBinContent(j+1, fDerivativeScint[j]);
  }

  double ct = 0;
  for( j = 0; j< fMaxReflects; j++) {
    hRWater.Fill(ct,fRWater[j]);
    hRScint.Fill(ct,fRScint[j]);
    ct = ct + 1./double(fMaxReflects);
  }

  f.cd();
  h1.Write();h1x.Write();dh1.Write();
  h2.Write();h2x.Write();dh2.Write();
  hRScint.Write(); hRWater.Write();
  f.Close();
#endif // DEBUG
}

void MultiPathScintWater::Initialise(const std::string&) {
}

void MultiPathScintWater::SetSeed(const DS::FitResult&) {
}

DS::FitResult MultiPathScintWater::MakeResult(const std::vector<double>& pars) {
  TVector3 fittedVertex(pars[0], pars[1], pars[2]);
  DS::FitVertex vertex;
  vertex.SetPosition(fittedVertex);
  vertex.SetTime(pars[3]);
  DS::FitResult result;
  result.SetVertex(0, vertex);
  return result;
}

bool MultiPathScintWater::ValidResult(const std::vector<double>& pars) {
  double r2 = pars[0]*pars[0] + pars[1]*pars[1] + pars[2]*pars[2];
  return r2 < fPSUPRadius2;
}

std::vector<double> MultiPathScintWater::GetStart() {
  // Could set up a seed here, but I'm not using it for now (from aTime and aVertex)
#ifdef MULTIPATH_FIXED_SEED
  double ran0 = step;
  double ran1 = 2.0*step - 1.0;
  double ran2Pi = step*2.0*CLHEP::pi;
  double t = 100.0 + 200.0*step;
  step += 0.0001;
  if (step > 1.0) step = 0.0001;
#else
  double ran0 = CLHEP::RandFlat::shoot( 0.0, 1.0 );
  double ran1 = CLHEP::RandFlat::shoot( -1.0, 1.0 );
  double ran2Pi = CLHEP::RandFlat::shoot( 0.0, 2.0*CLHEP::pi );
  double t = CLHEP::RandFlat::shoot( 100.0, 300.0 );
#endif

  double r = pow(ran0, 1.0/3.0) * 8390; // mm
  double costheta = ran1;
  double sintheta = sqrt(1.0 - costheta*costheta);
  std::vector<double> v;
  v.push_back(r*sintheta*cos(ran2Pi));
  v.push_back(r*sintheta*sin(ran2Pi));
  v.push_back(r*costheta);
  // in principle, MC and data could have different PMT hit time offset,
  // but in latest incarnation of MultiPathFitter code, this doesn't seem to be the case
  v.push_back(t);
  return v;
}

/// returns Likelihood and derivative for a given time of flight t.
double MultiPathScintWater::LAnddLdt( double tDiff, double tof, double &dLdt ) {
  double L = 0.0;
  fRes = tDiff - tof;
  if( fRes > fMaxPMTtRes ) fRes = fMaxPMTtRes; // longer than valid time (late light, crosstalk?)
  if( fRes < -fMaxPMTtRes ) fRes = -fMaxPMTtRes; // eariler than valid time

  int ibin = (int) floor( ( fRes + 100.0 ) * 4.0 );
  if(ibin > 0 && ibin < fEntriesTimeCh - 2) {
    double dt = fRes-fPDFxCh[ibin];
    L = (fPDFCh[ibin] + fDerivativeCh[ibin] * dt);
    if (dt < 0.125) {
      double u = 0.5 + dt * 4.0;
      dLdt = -u*fDerivativeCh[ibin] - (1.0-u)*fDerivativeCh[ibin-1];
    } else {
      double u = 4.0 * dt - 0.5;
      dLdt = -(1.0-u)*fDerivativeCh[ibin] - u*fDerivativeCh[ibin+1];
    }// end of the if statement if dt<0.125
  } else if (ibin <= 0) {
    L = fPDFCh[0];
    dLdt = -fDerivativeCh[0];
  } else {
    L = fPDFCh[fEntriesTimeCh - 2];
    dLdt = -fDerivativeCh[fEntriesTimeCh - 2];
  }
  return L;
}

double MultiPathScintWater::LAnddLdtScint( double tDiff, double tof, double &dLdt ) {
  double L = 0.0;
  fRes = tDiff - tof;
  if( fRes > fMaxPMTtRes ) fRes = fMaxPMTtRes; // longer than valid time (late light, crosstalk?)
  if( fRes < -fMaxPMTtRes ) fRes = -fMaxPMTtRes; // eariler than valid time

  int ibin = (int) floor( ( fRes + 100.0 ) * 4.0 );
  if(ibin > 0 && ibin < fEntriesTimeCh - 2) {
    double dt = fRes-fPDFxScint[ibin];
    L = (fPDFScint[ibin] + fDerivativeScint[ibin] * dt);
    if (dt < 0.125) {
      double u = 0.5 + dt * 4.0;
      dLdt = -u*fDerivativeScint[ibin] - (1.0-u)*fDerivativeScint[ibin-1];
    } else {
      double u = 4.0 * dt - 0.5;
      dLdt = -(1.0-u)*fDerivativeScint[ibin] - u*fDerivativeScint[ibin+1];
    }// end of the if statement if dt<0.125
  } else if (ibin <= 0) {
    L = fPDFScint[0];
    dLdt = -fDerivativeScint[0];
  } else {
    L = fPDFScint[fEntriesTimeCh- 2];
    dLdt = -fDerivativeScint[fEntriesTimeCh - 2];
  }
  return L;
}

double MultiPathScintWater::Calculate(const std::vector<double>& pmt,
                                      const std::vector<double>& par,
                                      std::vector<double>& diff) {
  // look up position each time - is this efficient enough?
  // const TVector3& pos = DU::Utility::Get()->GetPMTInfo().GetPosition(pmt.GetID());
  // double tpmt = pmt.GetTime();

  double dx = pmt[0] - par[0];
  double dy = pmt[1] - par[1];
  double dz = pmt[2] - par[2];
  double dr = sqrt(dx*dx + dy*dy + dz*dz); // | pmtpos - X0|

  TVector3 fVertex( par[0], par[1], par[2] ); // X0 = (x0,y0,z0), trial vertex
  TVector3 fPMTpos( pmt[0], pmt[1], pmt[2] );
  /// following is to calculate time residual by: PMTtime - tof - trial time;
  // fRes = pmt[3] - tof - par[3];
  // in partial fitter, tof is calculated separately in different ways in different situations;
  // instead of fRes, (tDiff, tof) is passed to the likelihood calculation
  double tDiff = pmt[3] - par[3];// tDiff = hitTime - trialTime; then fRes = tDiff - tof

  double L = 0.0;
  double dLdt = 0.0, dLdx = 0.0, dLdy = 0.0, dLdz = 0.0, dLdt0 = 0.0;
  double fx = par[0]; // trial x0
  double fy = par[1]; // trial y0
  double fz = par[2]; // trial z0
  fDepth = fWaterLevel - fz; // L - z0
  fHeight = pmt[2] - fWaterLevel; // Zpmt - L

  if( fWaterLevel >= -fAVRadius + fZoff && fWaterLevel <= fAVRadius + fZoff ) { // valid water level
    double tof = 0.0, alpha = 0.0, reflected = 0.0;

    int layout = 0; // check PMT and vertex layout, 0 is by default: except for geo 1, when scintpath == 0, always direct path in water
    if( fDepth>0 ) {
      if( fHeight<0 ) layout = 1; // geo 1: both PMT and vertex are below interface: find direct and reflected light paths
      else layout = 2; // geo 2: vertex below interface, PMT above
    }
    else {
      layout = 3; // geo 3, 4: vertex above interface, PMT above or below
    }

    /// evaluate the distanceInScint = scintpathInNeck + scintpathInAV
    double distInScint = 0;

    /// evaluate line-AV sphere intersection first
    double scintpathInAV = 0;
    double udDotVtx = fx*dx/dr + fy*dy/dr + (fz-fZoff)*dz/dr; // (P-X0)/|P-X0| * (X0 - oAV)
    double vtx2 = fx*fx + fy*fy + (fz-fZoff)*(fz-fZoff); // (X0 - oAV)*(X0 - oAV)
    double rVertex = sqrt( fx*fx + fy*fy + (fz-fZoff)*(fz-fZoff) );
    double sqrVal = udDotVtx*udDotVtx - vtx2 + fAVRadiusOuter*fAVRadiusOuter;
    double a1 = 0, a2 = 0, aplus = 0, abig = 0, asmall = 0;
    int passAVcheck = 0;
    if( sqrVal>0 ) { // line passes AV sphere
      passAVcheck = 1; // use to distinguish geo 0 and geo 1; for geo 1, line passes AV and count reflection
      /// find the line-sphere intersect points; a1, a2 are the path lengths btw vertex and intersection points
      a1 = -udDotVtx + sqrt(sqrVal);
      a2 = -udDotVtx - sqrt(sqrVal);
      /// find the line-interface plane intersection; a3 is the path length btw vertex and intersection point
      double a3 = -9999; // a3 is not defined if lightPath is parallel to interface
      if( fabs(dz)>fBoundary ) a3 = ( fWaterLevel-fz )*dr/dz; // well-defined a3

      if( rVertex<fAVRadiusOuter ) { // vertex inside the AV, equivalent to a1*a2<0
        if ( a1*a2<0 ) {
          aplus = ( a1>0 || a2>0 )? a1:a2;
          if( a3>0 ) { // a3>0, hit interface plane
            if( fz<fWaterLevel ) { // vertex below
              if( a3<aplus ) scintpathInAV = aplus - a3;
            }
            else { // vertex above
              if( a3<aplus ) scintpathInAV = a3;
              else scintpathInAV = aplus;
            }
          }
          else { // a3<=0, not hit interface plane
            if( fz>=fWaterLevel ) scintpathInAV = aplus; // vertex must above
          }
        }
      } // vertex in AV

      else { // rVertex>=fAVRadiusOuter, vertex in external
        if( a1>0 && a2>0) {
          abig = a1>a2? a1:a2;
          asmall = a1<a2? a1:a2;
          double zbig = fz + abig*dz/dr; // z position of the intersection points
          double zsmall = fz + asmall*dz/dr;
          if( zsmall>=fWaterLevel && zbig>=fWaterLevel )
            scintpathInAV = abig - asmall;
          if( zsmall<fWaterLevel && zbig>fWaterLevel && a3>0 )
            scintpathInAV = abig - a3;
          if( zsmall>fWaterLevel && zbig<fWaterLevel && a3>0 )
            scintpathInAV = a3 - asmall;
        } // ensure a1 and a2 are positive
      } // vertex in external
    } // pass through AV

    /// evaluate line-neck cylinder intersection
    double scintpathInNeck = 0;
    if( fNeckpathEnable ) { // turn off to ignorng high z/neck events

      double fZneckHi = rpsup; // top Z of the neck = Zpsup
      double bneck = (par[0]*dx+par[1]*dy);
      double sqrValneck = bneck*bneck - ( par[0]*par[0]+par[1]*par[1]-fNeckRadiusOuter*fNeckRadiusOuter )*( dx*dx + dy*dy );

      if( sqrValneck > 0 ) { // check the line passes through the neck
        double aneck1 = dr*( -bneck + sqrt(sqrValneck) )/( dx*dx+dy*dy );
        double aneck2 = dr*( -bneck - sqrt(sqrValneck) )/( dx*dx+dy*dy );

        if( aneck1*aneck2<0 ) { // vertex inside the neck
          double aneckplus = ( aneck1>0 || aneck2>0 )? aneck1:aneck2;
          double zplus = fz + aneckplus*dz/dr;
          if( zplus < fZneckHi && zplus > fZneckLo ) {
            scintpathInNeck = aneckplus;
          }
          if( zplus < fZneckLo ) {
            if( a1>0 && a2>0 ) { // must hits the AV at 2 points and outside the AV
              scintpathInNeck = asmall;
            }
          }
        } // vertex inside the neck

        if( aneck1>0 && aneck2>0 ) {
          double aneckbig = aneck1>aneck2? aneck1:aneck2;
          double anecksmall = aneck1<aneck2? aneck1:aneck2;
          double zbig = fz + aneckbig*dz/dr;
          double zsmall = fz + anecksmall*dz/dr;
          if( zbig<fZneckHi && zbig>fZneckLo && zsmall<fZneckHi && zsmall>fZneckLo )
          {
            scintpathInNeck = aneckbig - anecksmall;
            // an extreme condition: pass thru neck as well as the "virtual" sphere inside neck
            if(a1>0 && a2>0) {
              double zAVbig = fz + abig*dz/dr, zAVsmall = fz + asmall*dz/dr;
              if(zAVbig>=fZneckLo && zAVsmall>=fZneckLo) scintpathInNeck = aneckbig-anecksmall-(abig-asmall);
            }
          }

          if( zbig<fZneckLo && zsmall>fZneckLo && zsmall<fZneckHi ) {
            if(a1>0 && a2>0) { // must pass the AV at 2 points and outside the AV
              scintpathInNeck = asmall - anecksmall;
            }
          }

          if( zsmall<fZneckLo && zbig>fZneckLo && zbig<fZneckHi) {
            if(a1*a2<0) {// inside AV
              scintpathInNeck = aneckbig - aplus;
            }
            if(a1>0 && a2>0) { // outside AV
              scintpathInNeck = aneckbig - abig;
            }
          }

        } // vertex outside the neck, aneck1>0 and aneck2>0
      } // line hits the neck
    } // enalbe neck path calculation

   distInScint = scintpathInAV + scintpathInNeck;
    tof = ( dr - distInScint)/fSpeedOfLightWater + distInScint/fSpeedOfLightScint;
    if( distInScint == 0 ) { // if light path is always in water; except for geo 1, the other layout go to default
      if( layout != 1 ) layout = 0;
      else {
        if (!passAVcheck) layout = 0;// for geo1, if lightpath never passes AV, go to default
      }
    }

    switch (layout) {
      case 0: // geo 0: default, MPW like: direct light path in water only and never pass AV
      {
        tof = dr/fSpeedOfLightWater;
        L = LAnddLdt(tDiff, tof, dLdt);
        dLdt0 = dLdt/L;
        dLdx = -dLdt0*dx/(dr*fSpeedOfLightWater);
        dLdy = -dLdt0*dy/(dr*fSpeedOfLightWater);
        dLdz = -dLdt0*dz/(dr*fSpeedOfLightWater);
        L = log(L);
        break;
      }

      case 1: // geo 1: direct light + reflection light, photons pass the detector but always in water
      {
        tof = dr/fSpeedOfLightWater; // direct time
        alpha = fDepth/(fDepth-fHeight);  // find fractional transverse distance between vertex and water intersection
        TVector3 s = (1-alpha)*fVertex + alpha*fPMTpos;  // find a point directly below the intersection with water
        s.SetZ(fWaterLevel);  // move intersection to water surface.
        TVector3 incident = s-fVertex;
        TVector3 reflectedRay = fPMTpos-s;
        double treflected = (incident.Mag()+reflectedRay.Mag())/fSpeedOfLightWater;  // time of flight for reflected ray
        incident = incident.Unit();
        double ct = incident.Z()*fMaxReflects;  // find cos(theta) between vertical and photon trajectory, multiply by fMaxReflects to look up in table
        int ix=(int)ct; double u=ct-ix;
        reflected = (1-u)*fRWater[ix]+u*fRWater[ix+1];  // find reflection probability for this cos(theta), interpolated
        double dLdt1 = 0.0, dLdt2 = 0.0;
        L = ( LAnddLdt(tDiff, tof, dLdt1) + reflected*LAnddLdt(tDiff, treflected, dLdt2) ); // first term is for direct light that never hits the water surface; second term for water-water reflections
        dLdt1 = dLdt1/L;
        dLdt2 = reflected*dLdt2/L;
        L = log(L/(1+reflected));
        //components are proportional to dt/dx,dt/dy,dt/dz; dt/dt0=-1; incident is already set to a unit vector
        dLdx = -( dLdt1*dx/dr + dLdt2*incident.X() )/fSpeedOfLightWater;
        dLdy = -( dLdt1*dy/dr + dLdt2*incident.Y() )/fSpeedOfLightWater;
        dLdz = -( dLdt1*dz/dr + dLdt2*incident.Z() )/fSpeedOfLightWater;
        dLdt0 = dLdt1 + dLdt2;
        break;
      }

      case 2:
      {
        /// geo 2: PMT is above water, vertex below
        L = LAnddLdtScint(tDiff, tof, dLdt);
        dLdt0 = dLdt/L;
        L = log(L);
        /// derivatives are simplified as the vertex moves a little bit in the opposite direction to the light path in water
        dLdx = -dLdt0*dx/( dr*fSpeedOfLightWater );
        dLdy = -dLdt0*dy/( dr*fSpeedOfLightWater );
        dLdz = -dLdt0*dz/( dr*fSpeedOfLightWater );
        break;
      }

      case 3:
      {
        /// geo 3: vertex above: PMT below or above
        L = LAnddLdtScint(tDiff, tof, dLdt);
        dLdt0 = dLdt/L;
        L = log( L );
        /// derivatives are simplified as the vertex moves a little bit in the opposite direction to the light path in water or scint
        double fSpeedOfLightEff = fVertex.Mag()<fAVRadiusOuter? fSpeedOfLightScint:fSpeedOfLightWater;
        dLdx = -dLdt0*dx/( dr*fSpeedOfLightEff );
        dLdy = -dLdt0*dy/( dr*fSpeedOfLightEff );
        dLdz = -dLdt0*dz/( dr*fSpeedOfLightEff );
        break;
      }
    } // switch layout

    diff.assign(4, 0.0);// sBeta
    diff[0] = dLdx;
    diff[1] = dLdy;
    diff[2] = dLdz;
    diff[3] = dLdt0;
   } // valid water level
   else {
    fHeight = pmt[2];
    Log::Die( "Not a valid water height ");
   }

#ifdef DEBUG
//  debug << "AWP::CL pmt=" << pmt[0] << "," << pmt[1] << "," << pmt[2] << ")\n";
//  debug << "        fRes=" << fRes << " ," << dLdx << " "<< dLdy << " " <<dLdz << "\n";
//  debug << "        L=" << L << " diff=(" << diff[0] << "," << diff[1] << "," << diff[2] << "," << diff[3] << ")\n";
#endif

  return L;
}