PhotonLookupTable.cxx

Go to the documentation of this file.

#include "PhotonLookupTable.h"
#include <TGraph.h>
#include <TH1F.h>
#include <TRandom.h>
#include <iostream>

using namespace std;

ClassImp(PhotonLookupTable)

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

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

void PhotonLookupTable::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 PhotonLookupTable::Rebuild()
{
  //
  // Rebuilds arrays if the DB result changes.
  //
  fx.clear();
  fy.clear();
  fint.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)->GetX() ] =  fResPtr.GetRow(i)->GetY();
  }
  
  double integral = 0;
  for( sortMapItr = sortMap.begin(); sortMapItr != sortMap.end(); sortMapItr++) {
    fx.push_back(sortMapItr->first);
    fy.push_back(sortMapItr->second);
    integral += sortMapItr->second;
    fint.push_back(integral);
  }
  fN = fx.size();

}

Int_t PhotonLookupTable::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 PhotonLookupTable::Interpolate( Double_t value )
{
  //
  // Find the y corresponding to x. Interpolate.
  //
  int i = BinarySearch(fx,value,true);
  //cout << "Binary search: " << i << "\t" << fx[i] << endl;
  return (value-fx[i])/(fx[i+1]-fx[i])*(fy[i+1]-fy[i]) + fy[i];
}

Double_t PhotonLookupTable::Pick(Double_t r)
{
  //
  // Pick a random x from the integral of y
  // using random number r.
  //
  double pick = r* fint[fN-1];
  int i = BinarySearch(fint,pick);

  double r2 = (pick-fint[i])/(fint[i+1]-fint[i]); // is 0 to 1
  //double yy1 = fy[i]*fy[i];
  //double yy2 = fy[i+1]*fy[i+1];
  //double dy = fy[i+1]-fy[i];
  double fraction;
  //if(dy!=0) 
  //fraction = (-fy[i] + sqrt(yy1+(yy2-yy1)*r2))/(fy[i+1]-fy[i]);
    //else
  fraction = r2;
  double val = fx[i] + (fx[i+1]-fx[i])*fraction;
  
  //cout << "Pick: " 
  //    << fx[i] << "\t" 
  //   << (fraction) << "\t" 
  //   << fx[i+1]-fx[i] << "\t" 
  //   << val << endl;

  return val;
}

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

TGraph* PhotonLookupTable::GetGraph()
{
  // 
  // Draw a TGraph of the table.
  //
  if(!fGraph) {
    double* x = new double[fN];
    double* y = new double[fN];
    
    for(UInt_t i=0;i<fN;i++) {
      x[i] = fx[i];
      y[i] = fy[i];
    }

    fGraph = new TGraph(fN,x,y);
    
    delete [] x;
    delete [] y;
  }
      
  return fGraph;
}

void PhotonLookupTable::Draw(Option_t* option)
{
  GetGraph()->Draw(option);
}

TH1* PhotonLookupTable::MakeRandomHisto(Int_t n)
{
  TH1* h = new TH1F("tempTableHisto","tempTableHisto",fN*10,fx[0],fx[fN-1]);
  TRandom rand;
  for(int i=0; i<n; i++) {
    h->Fill(Pick(rand.Rndm()));
  }
  h->Draw();
  return h;
}


void PhotonLookupTable::Print(Option_t*) const
{
  // 
  // Print out the table.
  //
  cout << fTableName << endl;
  for(UInt_t i=0; i<fN; i++) {
    cout << "\t" << fx[i] << "\t" << fy[i] << endl;
  }
}

---