RangeLookupTable.cxx

Go to the documentation of this file.

#include "RangeLookupTable.h"
#include <TGraph.h>
#include <TH1F.h>
#include <TRandom.h>
#include <iostream>
#include <iomanip>

using namespace std;

ClassImp(RangeLookupTable)

VldContext RangeLookupTable::fgDefaultContext(Detector::kNear,
                                              SimFlag::kData,
                                              VldTimeStamp());

RangeLookupTable::RangeLookupTable(const char* tablename, 
                                     const VldContext& vc)
  : fTableName(tablename),
    fResPtr(string(tablename), vc),
    fGraph(0)
{
  Rebuild();
  fVldRange =  fResPtr.GetValidityRec()->GetVldRange();
}

void RangeLookupTable::Reset(const VldContext& vc)
{
  //
  // Switches to a new validity context.
  // 
  fResPtr.NewQuery(vc,0);
  VldRange newrange = fResPtr.GetValidityRec()->GetVldRange();

  // See if the table has somehow changed.
  if( (newrange.GetTimeStart()!=fVldRange.GetTimeStart()) ||
      (newrange.GetTimeEnd()!=fVldRange.GetTimeEnd()) ||
      ( newrange.GetDetectorMask()!=fVldRange.GetDetectorMask()) ||
      ( newrange.GetSimMask()!=fVldRange.GetSimMask()) )  {

    // Table changed. Rebuild the data.
    fVldRange = newrange;
    Rebuild(); 
  }
}

void RangeLookupTable::Print(Option_t* /* option */) const
{
  //
  // Print the values
  //
  std::cout << "RangeLookupTable " << fTableName 
            << " has " << fN << " points." << std::endl;
  std::cout << "  indx      Range in Fe     Mu Momentum        Integral"
            << std::endl;
  std::cout << "            (g/cm^2)        (GeV/c)" 
            << std::endl;
  for (UInt_t i=0; i < fN ; ++i) {
    std::cout << " [ " << std::setw(3) << i << "] "
              << setw(15) << setprecision(7) << frange[i] << " " 
              << setw(15) << setprecision(7) << fmomentum[i] << " " 
              << std::endl;
  }
  std::cout << resetiosflags(ios::floatfield); // reset to "%g" format
  std::cout << setprecision(6);
}


void RangeLookupTable::Rebuild()
{
  //
  // Rebuilds arrays if the DB result changes.
  //
  frange.clear();
  fmomentum.clear();
 
  // Use a map to sort the data as it comes in.
  std::map<double,double> sortMap;
  std::map<double,double>::iterator sortMapItr;

  for(UInt_t i=0; i<fResPtr.GetNumRows(); i++) {
    sortMap[ fResPtr.GetRow(i)->GetRange() ] =  fResPtr.GetRow(i)->GetMomentum();
  }

   for( sortMapItr = sortMap.begin(); sortMapItr != sortMap.end(); sortMapItr++) 
{
    frange.push_back(sortMapItr->first);
    fmomentum.push_back(sortMapItr->second);
  }
  fN = frange.size();

}

Int_t RangeLookupTable::BinarySearch( std::vector<Double_t> &array,
                                       Double_t value, 
                                       Bool_t doGuess )
{
  // Find the nearest entry below value.

  // Binary search. Lifted mercilessly from TMath.cxx
  // fx is supposed  to be sorted prior to this call.
  // If match is found, function returns position of element.
  // If no match found, function gives nearest element smaller than value.

  Int_t nabove, nbelow, middle;
  UInt_t n = array.size();
  if(value<=array[0]) return 0;
  if(value>=array[n-2]) return n-2;
  nabove = n+1;
  nbelow = 0;
  
  // Make the guess.
  if(doGuess) {
    Double_t guess = (value-array[0])/(array[n-1]-array[0]);
    if(guess<0.5) middle = (Int_t)(guess*n) +2;
    else          middle = (Int_t)(guess*n) -2;

    // Ensure the guess paramater isn't silly.
    if(middle>(Int_t)(n-1)) middle = n-1; 
    if(middle<1) middle = 1;
      
  } else {
    middle = (nabove-nbelow)/2;
  }
  
  while(nabove-nbelow > 1) {
    if (value == array[middle-1]) return middle-1;
    if (value  < array[middle-1]) nabove = middle;
    else                          nbelow = middle;
    middle = (nabove+nbelow)/2;
  }
  return nbelow-1;
}


Double_t RangeLookupTable::Interpolate( Double_t value )
{
  //
  // Find the y corresponding to x. Interpolate.
  //
  int i = BinarySearch(frange,value,true);
  //cout << "Binary search: " << i << "\t" << frange[i] << endl;
  return (value-frange[i])/(frange[i+1]-frange[i])*(fmomentum[i+1]-fmomentum[i]) + fmomentum[i];
}

Double_t RangeLookupTable::GetMaximum()
{
  double max = -1e20;
  for(UInt_t i=0;i<fN;i++) {
    if(fmomentum[i]>max) max = fmomentum[i];
  }
  return max;
}

---