// ----------------------------------------------------------------------------
// $Id: LCDLookUp2d.cxx,v 1.1 2001/05/10 17:56:02 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDLookUp2d.cxx,v $
// Revision 1.1  2001/05/10 17:56:02  toshi
// Changes C++ name convention from *.C to *.cxx to allow Windows usage.
//
//
// two dimensional lookup table
// Version 1.0  7-Jan-1999 Richard creation
// Version 1.1 20-Apr-2000 Toshi modified to privent negative
//                         diagonal error matrix element...

#include "LCDLookUp2d.h"

//______________________________________________________________________
//
// LookUp2d
//
// utility class which reads a 2-d table and uses it for interpolation

ClassImp(LCDLookUp2d)
 LCDLookUp2d::LCDLookUp2d(TString name, 
			 Double_t* key1, Double_t* key2, Double_t* mat, 
			 Int_t numBins1, Int_t numBins2) :
  m_name(name), m_key1(key1), m_key2(key2), m_matrix(mat), 
  m_numBins1(numBins1), m_numBins2(numBins2) {
};

 LCDLookUp2d::~LCDLookUp2d() {};


 Double_t LCDLookUp2d::interpolateVal(Double_t val1, Double_t val2, Int_t* rc) {
  // interpolates table given input val1, val2

  //  const Int_t badVal1LowerBound = 1;
  const Int_t badVal1UpperBound = 2;
  //  const Int_t badVal2LowerBound = 4;
  const Int_t badVal2UpperBound = 8;
  const Int_t goodBounds = 0;
    
  *rc = goodBounds;
  Int_t index1;
  Int_t index2;
  Double_t t;
  Double_t u;

  if (val1 > m_key1[m_numBins1-1] ) {
    *rc   += badVal1UpperBound;
    index1 = m_numBins1 - 1;
    t      = 0;
  } else {
    index1 = TMath::BinarySearch(m_numBins1,m_key1,val1);
    index1 = ( index1 < 0 ) ? 0 : index1;
    t = (val1 - m_key1[index1])/(m_key1[index1+1] - m_key1[index1]);
  } 
  if (val2 > m_key2[m_numBins2-1] ) {
    *rc   += badVal2UpperBound;
    index2 = m_numBins2 - 1;
    u      = 0;
  } else {
    index2 = TMath::BinarySearch(m_numBins2,m_key2,val2);
    index2 = ( index2 < 0 ) ? 0 : index2;
    u = (val2 - m_key2[index2])/(m_key2[index2+1] - m_key2[index2]);
  }

  t = ( t < 0 ) ? 0 : t ;
  u = ( u < 0 ) ? 0 : u ;

  Double_t y1 = 0;
  Double_t y2 = 0;
  Double_t y3 = 0;
  Double_t y4 = 0;
  y1 = m_matrix[index1*m_numBins2+index2];
  if (index1+1<m_numBins1) {
    y2 = m_matrix[(index1+1)*m_numBins2+index2];
  }
  if (index1+1<m_numBins1 && index2+1<m_numBins2) {
    y3 = m_matrix[(index1+1)*m_numBins2+index2+1];
  }
  if (index2+1<m_numBins2) {
    y4 = m_matrix[index1*m_numBins2+index2+1];
  }
  
  Double_t yNew = (1.-t)*(1.-u)*y1 + t*(1.-u)*y2 + t*u*y3 + (1.-t)*u*y4;
  
  return yNew;
};