Calibrator.cxx

Go to the documentation of this file.

// $Id: Calibrator.cxx,v 1.7 2009/10/19 17:15:21 jpaley Exp $
//
// Singleton class that handles calibrations for a particular 
// validity context (eg: detector+timestamp).
//

#include "Calibrator/Calibrator.h"
#include "CMap/CMap.h"
#include "JobControl/Exception.h"
#include "Geometry/geo.h"
#include "Database/Table.h"
#include "TMath.h"
#include "TSystem.h"
#include "TFile.h"
#include <iostream>

using namespace std;
using namespace calib;

Calibrator* Calibrator::fInstance = 0;

//------------------------------------------------------------

Calibrator::Calibrator() : 
  fCurrTimeStamp(0),
  fCurrDetector(rawdata::kNear),
  fCurrGeo(0),
  fPEACperMIP(0), fPEACperMIPErr(0),
  fLoadedFineTimeDistros(0)
{
  fADCtoPE.clear();
  fMIPtoGeV.clear();
  fMIPtoGeVErr.clear();
  fFTDistNDX.clear();
  fFTDistNDY.clear();
  fFTDistFDX.clear();
  fFTDistFDY.clear();
}


//------------------------------------------------------------

void Calibrator::SetContext(TTimeStamp timestamp,
                            rawdata::det_id_ detector,
                            const geo::Geometry &geometry) {
  fCurrTimeStamp = timestamp;
  fCurrDetector = detector;
  fCurrGeo = &geometry;
  LoadFineTimeDistros();
  /*
  LoadADCtoPE();
  LoadPEtoSigCor();
  LoadSCtoPEAC();
  LoadPEACtoMIP();
  LoadMIPtoGeV();
  */
}

//------------------------------------------------------------

Calibrator& Calibrator::Instance(TTimeStamp timestamp,
                                 rawdata::det_id_ detector,
                                 const geo::Geometry &geometry) {
  if (!fInstance)
    fInstance = new Calibrator;
  fInstance->SetContext(timestamp, detector, geometry);
  return *fInstance;
}

//------------------------------------------------------------
//   Load ADC to PE database constants
//------------------------------------------------------------

bool Calibrator::LoadADCtoPE()
{
  fADCtoPE.clear();
  
  std::string dbFile = "ADCtoPE.xml";

  db::Table* table = new db::Table(dbFile);
  
  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  // now fill STL map...
  std::string str1 = "Channel";
  std::string str2 = "PEperADC";
  int chan;
  float val;
  for (unsigned int i=0; i<table->size(); ++i) {    
    (*table)[i].Get(str1,chan);
    (*table)[i].Get(str2,val);
    fADCtoPE[chan] = val;
  }
  
  delete table;

  return true;
}

//------------------------------------------------------------
//   Load PE to SigCor database constants
//------------------------------------------------------------

bool Calibrator::LoadPEtoSigCor()
{
  fPEtoSigCor.clear();
  
  std::string dbFile = "PEtoSigCor.xml";

  db::Table* table = new db::Table(dbFile);

  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  // now fill STL map...
  std::string str1 = "Channel";
  std::string str2 = "SCperPE";
  int chan;
  float val;
  for (unsigned int i=0; i<table->size(); ++i) {    
    (*table)[i].Get(str1,chan);
    (*table)[i].Get(str2,val);
    fPEtoSigCor[chan] = val;
  }
  
  delete table;
  
  return true;
}

//------------------------------------------------------------
//   Load SigCor to PEAC database constants
//------------------------------------------------------------

bool Calibrator::LoadSCtoPEAC()
{
  fSCtoPEAC.clear();
  
  std::string dbFile = "SigCortoPEAC.xml";

  db::Table* table = new db::Table(dbFile);

  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  // now fill STL map...
  std::string str1 = "Channel";
  std::string str2 = "A1";
  std::string str3 = "A2";
  std::string str4 = "L1";
  std::string str5 = "L2";
  int chan;
  for (unsigned int i=0; i<table->size(); ++i) {    
    attlen a;
    (*table)[i].Get(str1,chan);
    (*table)[i].Get(str2,a.a1);
    (*table)[i].Get(str2+"Err",a.a1sig);
    (*table)[i].Get(str3,a.a2);
    (*table)[i].Get(str3+"Err",a.a2sig);
    (*table)[i].Get(str4,a.l1);
    (*table)[i].Get(str4+"Err",a.l1sig);
    (*table)[i].Get(str5,a.l2);
    (*table)[i].Get(str5+"Err",a.l2sig);

    fSCtoPEAC[chan] = a;
  }
  
  delete table;
  
  return true;
}

//------------------------------------------------------------
//   Load PEAC to MIP database constant
//------------------------------------------------------------

bool Calibrator::LoadPEACtoMIP()
{
  fPEACperMIP = 1.;
  
  std::string dbFile = "PEACtoMIP.xml";

  db::Table* table = new db::Table(dbFile);

  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  std::string str1 = "PEACperMIP";
  (*table)[0].Get(str1,fPEACperMIP);

  delete table;
  
  return true;
}

//------------------------------------------------------------
//   Load MIP to GeV database constant
//------------------------------------------------------------

bool Calibrator::LoadMIPtoGeV()
{
  fMIPtoGeV.clear();
  
  std::string dbFile = "MIPtoGeV.xml";

  db::Table* table = new db::Table(dbFile);

  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  float val;
  char buf[20];
  std::string str1;
  for (int i=0; i<5; ++i) {    
    sprintf(buf,"par%d",i);
    str1 = buf;
    (*table)[0].Get(str1,val);
    fMIPtoGeV.push_back(val);
  }

  delete table;
  
  if (fMIPtoGeV.empty()) return false;

  return true;
}

//------------------------------------------------------------
//   Load Fine-time Distributions from dB
//------------------------------------------------------------

bool Calibrator::LoadFineTimeDistros() 
{
  fFTDistNDX.clear();
  fFTDistNDY.clear();
  fFTDistFDX.clear();
  fFTDistFDY.clear();

  // get values for ND first
  std::string dbFile = "FTDistroND.xml";

  db::Table* table = new db::Table(dbFile);

  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  float x,y;
  std::string strx = "x";
  std::string stry = "y";
  for (unsigned int i=0; i<table->size(); ++i) {
    (*table)[i].Get(strx,x);
    (*table)[i].Get(stry,y);
    fFTDistNDX.push_back(x);
    fFTDistNDY.push_back(y);
  }

  delete table;
  
  // now get values for FD 
  dbFile = "FTDistroFD.xml";

  table = new db::Table(dbFile);

  if (!table) return false;

  if (table->Connect()) {
    table->Init(fCurrTimeStamp.GetSec());
    
    if (table->empty()) return false;

    table->Disconnect();
  }

  for (unsigned int i=0; i<table->size(); ++i) {
    (*table)[i].Get(strx,x);
    (*table)[i].Get(stry,y);
    fFTDistFDX.push_back(x);
    fFTDistFDY.push_back(y);
  }

  delete table;

  fLoadedFineTimeDistros = 1;

  return true;
}

//------------------------------------------------------------
//   RawDigit-to-CellHit methods
//------------------------------------------------------------

void Calibrator::MakeCellHit(const RawDigit *rawdigit, CellHit *cellhit) {

  // copy inherited data members
  (*((RawDigit*)cellhit)) = (*rawdigit);

  cmap::CMap& cmap = cmap::CMap::Instance(fCurrTimeStamp,fCurrDetector);

  // new info
  cellhit->SetStatus(0); // for now, it always works!
  cellhit->SetCell  (cmap.GetCell(rawdigit));
  cellhit->SetPlane (cmap.GetPlane(rawdigit));
  cellhit->SetView  (fCurrGeo->Plane(cellhit->Plane()).View());
  cellhit->SetTNS   (GetTNS(rawdigit));
  cellhit->SetPE    (GetPE(rawdigit));
  cellhit->SetPECorr(cellhit->PE()*GetPECorrFactor(rawdigit));
}

//-----
float Calibrator::GetPE(const RawDigit *rawdigit) {
  // Look up gain from DB.  For now, must hard code the
  // number.  (This should not be put in any XML file.
  // That'll just confuse the issue.)
  float gain = 1.43;
  return rawdigit->ADC(0)/gain;
}

//-----
float Calibrator::GetPECorrFactor(const RawDigit* /*rawdigit*/) {
  // Look up channel-to-channel correction factors in DB.
  // For now, must hard code the factor.
  return 1.0;
}

//-----
float Calibrator::GetTNS(const RawDigit *rawdigit) {
  // You'd think this would be a simple function... Not anymore!
  // The new readout scheme allows for super time resolution, and
  // this is the code that gets it.

  Double_t t = 0;

  // First, see if we're using the new-style readout.  This
  // condition should be unnecessary in the future.

  if (rawdigit->NADC()<=1) {

    // Old style readout.  Return old style answer:
    t = rawdigit->TDC(0)*62.5;

  }
  else {

    // New style readout.  Use ADC info to determine fine timing.
    // See doc-???? for details.

    // ADC readings before and after peak, scaled to peak,
    // yield "x" and "y":
    Double_t peakadc = rawdigit->ADC(0);
    Double_t x = (rawdigit->ADC(1))/peakadc;
    Double_t y = (rawdigit->ADC(2))/peakadc;
    
    // These constants belong in the the DB.  For now, hard code them.
    // Tuned up by RBP, 2009-10-07
    Double_t Xa = 0.0;
    Double_t Ya = 1.0;
    Double_t Xb = 0.0;
    Double_t Yb = (fCurrDetector==rawdata::kNear?0.90:0.885);
    Double_t Xc = 1.0;
    Double_t Yc = 0.800;

    // These could be evaluated offline and stuck in the DB.
    // For now, I'll evaluate them explicitly here so the math
    // isn't lost:
    Double_t L1 = sqrt((Yb-Ya)*(Yb-Ya)+(Xb-Xa)*(Xb-Xa));
    Double_t L2 = sqrt((Yc-Yb)*(Yc-Yb)+(Xc-Xb)*(Xc-Xb));

    // distance to first line segment:
    Double_t U1 = ((x-Xa)*(Xb-Xa)+(y-Ya)*(Yb-Ya))/L1/L1;
    if (U1<0) U1=0;
    if (U1>0.99999) U1=0.99999;
    Double_t PX1 = Xa + U1*(Xb-Xa);
    Double_t PY1 = Ya + U1*(Yb-Ya);
    Double_t D1 = sqrt((PX1-x)*(PX1-x)+(PY1-y)*(PY1-y));

    // distance to second line segment:
    Double_t U2 = ((x-Xb)*(Xc-Xb)+(y-Yb)*(Yc-Yb))/L2/L2;
    if (U2<0) U2=0;
    if (U2>0.99999) U2=0.99999;
    Double_t PX2 = Xb + U2*(Xc-Xb);
    Double_t PY2 = Yb + U2*(Yc-Yb);
    Double_t D2 = sqrt((PX2-x)*(PX2-x)+(PY2-y)*(PY2-y));

    // get the better of the two
    Double_t U;
    if (D1<D2) U = U1*L1/(L1+L2);
    else       U = (L1+U2*L2)/(L1+L2);

    // map best U value to a fine time:
    vector<float>& ftx = fFTDistNDX;
    vector<float>& fty = fFTDistNDY;
    if (fCurrDetector == rawdata::kFar) {
      ftx = fFTDistFDX;
      fty = fFTDistFDY;
    }

    unsigned int bin=1;
    if (U < ftx[0])
      bin = 0;
    else
      for (; bin<ftx.size()-1; ++bin) {
        if (U >= ftx[bin-1] && U < ftx[bin]) break;
      }
    if (bin == ftx.size()) {
      // should never get in here
      std::cout << "In GetTNS(): bin=" << bin << " u=" << U << std::endl;
      bin=ftx.size()-1;
    }
    Double_t timefrac = 0.5*( fty[bin] + fty[bin-1]);

    // Gap between digitizations
    Double_t gap = (fCurrDetector==rawdata::kNear?125.0:500.0);

    t = rawdigit->TDC(0)*62.5 - timefrac*gap;
  }

  return t;
}


//------------------------------------------------------------
//   CellHit-to-RecoHit methods
//------------------------------------------------------------

void Calibrator::MakeRecoHit(const CellHit *cellhit, RecoHit *recohit, float w) {

  // copy inherited data members
  (*((CellHit*)recohit)) = (*cellhit);

  // new info:
  recohit->SetW(w);

  // Need position info to continue...
  Double_t xyzd[4];
  GetXYZD(cellhit,w,xyzd);

  recohit->SetX    (xyzd[0]);
  recohit->SetY    (xyzd[1]);
  recohit->SetZ    (xyzd[2]);

  recohit->SetT    (GetT    (cellhit,xyzd[3]));
  recohit->SetPEAtt(GetPEAtt(cellhit,xyzd[3]));

  recohit->SetMIP  (recohit->PEAtt() * GetMIPScale(cellhit));
  recohit->SetGeV  (recohit->MIP()   * GetGeVScale(cellhit));
}

//-----
float Calibrator::GetT(const CellHit *cellhit, float dist_to_readout) {
  // This needs to be calibrated and put in a DB table
  Double_t speedOfFiberTransport = 15.3; // cm/ns, "first principles" calib.
                                         // Differs from c/n due to zigzag
                                         // paths down fiber.  But, this
                                         // value may be way off (10%? 30%?).
  return (cellhit->TNS()-dist_to_readout/speedOfFiberTransport);
}

//-----
float Calibrator::GetPEAtt(const CellHit *cellhit, float dist_to_readout) {

  // These constants belong in the DB.  Don't put them in an XML file
  // in the meantime, as this will just confuse the issue...
  Double_t F0 = 5.744; // P. Shanahan tuning from old CalHit.xml
  Double_t L0 = 315.3; // .
  Double_t F1 = 1.546; // .
  Double_t L1 = 3329.; // .

  Double_t attenuation = F0*TMath::Exp(-dist_to_readout/L0) +
                         F1*TMath::Exp(-dist_to_readout/L1);

  return cellhit->PECorr() / attenuation;
}

//-----
float Calibrator::GetMIPScale(const CellHit* /* cellhit */ ) {
  // This needs to be calibrated and put in the DB, obviously!
  return 0.033;
}

//-----
float Calibrator::GetGeVScale(const CellHit* /* cellhit */ ) {
  // This needs to be calibrated and put in the DB, obviously!
  // Super-rough number based on 1 MIP = ~2 MeV/cm.
  return 0.013;
}

//-----
void Calibrator::GetXYZD(const CellHit *cellhit, float w, double *xyzd) {
  // Return (via 'xyzd' argument) position of deposition in world
  // coordinates and the distance to the readout, assuming deposition
  // has cell-local-z position = w.

  UInt_t plane = cellhit->Plane();
  UInt_t cell  = cellhit->Cell();
  if (plane>=fCurrGeo->NPlanes()) {
    throw geo::Exception(__FILE__,__LINE__,geo::Exception::BAD_PLANE_NUMBER);
  }
  const geo::PlaneGeo& p = fCurrGeo->Plane(plane);
  if (cell>=p.Ncells()) {
    throw geo::Exception(__FILE__,__LINE__,geo::Exception::BAD_CELL_NUMBER);
  }
  const geo::CellGeo& c = p.Cell(cell);

  // Should be replaced with alignment tables or whatnot.  Okay for now.

  // w is not localZ, which makes this a lot messier than just
  // c.GetCenter(xyzd,w) and c.DistToReadOut(w)
  c.GetCenter(xyzd);
  if (cellhit->View()==geo::kX)
    xyzd[1] = w;
  else
    xyzd[0] = w;
  Double_t localxyz[3];
  c.WorldToLocal(xyzd,localxyz);
  xyzd[3] = c.DistToReadOut(localxyz[2]);
}

---