QuadLinearityCalScheme.cxx

Go to the documentation of this file.

#include "QuadLinearityCalScheme.h"
#include "MessageService/MsgService.h"
#include "PulserCalibration/PulserConventions.h"
#include <limits>

ClassImp(QuadLinearityCalScheme)
CVSID("$Id: QuadLinearityCalScheme.cxx,v 1.10 2008/10/14 21:04:46 tagg Exp $");
//......................................................................

const int kFitType = 44;

QuadLinearityCalScheme::QuadLinearityCalScheme() 
{
 MSG("Calib",Msg::kVerbose) << "QuadLinearityCalScheme::QuadLinearityCalScheme" 
                            << endl;
 Registry r;
 r.Set("PulseToBucketCorrection",0.44);
 r.Set("BucketMode",1);

 InitializeConfig(r);
}

//......................................................................
QuadLinearityCalScheme::~QuadLinearityCalScheme()
{
}


//......................................................................

FloatErr QuadLinearityCalScheme::GetLinearized(FloatErr rawadc, 
                                               const PlexStripEndId& seid) const
{
  //Purpose: Convert from raw ADC to linearized ADC
  //It would find the Pin value in the gain curve of the requested channel for the given raw ADC value either through interpolation or extrapolation. Then multiply that pin value by the slope from a linear fit of the same gain curve in the linear region to take out the PMT saturation and get the linearized ADC.
  //In: raw ADC
  //Out: linearized ADC (siglin)

  // Get fit:
  const CalPulserFits* row = fFitTable.GetRowByIndex(seid.BuildPlnStripEndKey());
  
  Float_t beta = 0; // The nonlinearity parameter, in ADC^-1. Should be -ve
  Float_t dbeta = 0.; // The error on beta.
  Float_t maxadc = 0; // The max ADC value pulsed by the system
  Float_t meanres = rawadc; // The mean residual, in ADC counts.
  if(row) {
    beta = row->GetPar1();
    dbeta = row->GetPar2();
    maxadc = row->GetAdcMax();
    meanres = row->GetMeanRes();
  } else {
    if(fFitTable.GetNumRows()>0) {
      MAXMSG("Calib",Msg::kWarning,10) << "No CALPULSERFITS quadratic fit for " << seid.AsString() << endl;
      IncrementErrors(kLinCalibrator,kMissingRow,seid);
    }
  }

  if(fBucketMode) {
    beta /= fPulseToBucketCorrection;
    dbeta /= fPulseToBucketCorrection;
  }

  // Apply the nonlinearity correction.

  Float_t err2 = rawadc.GetError2(); // error squared
  Float_t cal = rawadc;
  
  // Total saturation limit:
  float s = 1.+4.*rawadc*beta;
  if(s<0) {
    IncrementErrors(kLinCalibrator,kBadInput,seid);
    s=0;
    err2 += pow(rawadc - 1.0/(-4.*beta),2); // Increase error by distance to saturation.
  } else {
    // The correct way to calibrate:
    if(beta<0) cal = 1.0/(2.*beta) * (sqrt(s) - 1.0);     
  }

  const float kMaxMeanRes = sqrt(std::numeric_limits<float>::max())*1e-6;
  if(meanres>kMaxMeanRes) meanres = kMaxMeanRes;
  if(cal>2500)  
    err2 += meanres*meanres; // Add on mean residual as error

  if(rawadc>maxadc) {
    // We're out of bounds. Apply an additional error
    // as the error on beta times the correction.
    IncrementErrors(kLinCalibrator,kDataInsufficient,seid);
    Float_t temp = (cal-rawadc)*(dbeta);
    if(beta!=0) temp = temp/beta;
    err2 += temp*temp;
  }
  FloatErr retval;
  retval.SetValueError2(cal,err2);
  return retval;
}

//......................................................................

FloatErr QuadLinearityCalScheme::DecalLinearity(FloatErr lin, 
                                                const PlexStripEndId& seid) const
{

  const CalPulserFits* row = fFitTable.GetRowByIndex(seid.BuildPlnStripEndKey());
  
  Float_t beta = 0; // The nonlinearity parameter, in ADC^-1. Should be -ve
  Float_t dbeta = 0.; // The error on beta.
  Float_t maxadc = 0; // The max ADC value pulsed by the system
  Float_t meanres = lin; // Mean residual, in ADC counts.
  if(row) {
    beta = row->GetPar1();
    dbeta = row->GetPar2();
    maxadc = row->GetAdcMax();
    meanres = row->GetMeanRes();
  } else {
    if(fFitTable.GetNumRows()>0) {
      MAXMSG("Calib",Msg::kWarning,10) << "No CALPULSERFITS quadratic fit for " << seid.AsString() << endl;
      IncrementErrors(kLinCalibrator,kMissingRow,seid);
    }
  }

  if(fBucketMode) {
    beta /= fPulseToBucketCorrection;
    dbeta /= fPulseToBucketCorrection;
  }

  // Apply the nonlinearity correction.

  Float_t x   = lin.GetValue();
  Float_t err2 = lin.GetError2(); // error squared
  Float_t nonlin = x * (1.0 + x*beta);
  
  const float kMaxMeanRes = sqrt(std::numeric_limits<float>::max())*1e-6;
  if(meanres>kMaxMeanRes) meanres = kMaxMeanRes;
  if(x>2500) 
    err2 += meanres*meanres;

  if(x>maxadc) {
    // We're out of bounds. Apply an additional error
    // as the error on beta times the correction.
    IncrementErrors(kLinCalibrator,kBadInput,seid);
    Float_t temp = (lin-nonlin)*(dbeta);
    if(beta!=0) temp = temp/beta;
    err2 += temp*temp;
  }


  FloatErr retval;
  retval.SetValueError2(nonlin,err2);
  return retval;
}

//......................................................................

void QuadLinearityCalScheme::ConfigModified() 
{
  Bool_t ok = true;
  ok = ok& GetConfig().Get("PulseToBucketCorrection",fPulseToBucketCorrection);
  ok = ok& GetConfig().Get("BucketMode",fBucketMode);
  if(!ok) MSG("Calib",Msg::kError) << "QuadLinearityCalScheme has a configration problem!" << endl;
}

//......................................................................
//
void QuadLinearityCalScheme::PrintConfig( std::ostream& os ) const
{
  os << "  Calibration mode:           ";
  if(fBucketMode) os << "Bucket-by-bucket" << endl;
  else            os << "Total Pulseheight" << endl;
  os << "  Pulse-to-Bucket correction: " << fPulseToBucketCorrection << endl; 
}

//......................................................................

void QuadLinearityCalScheme::DoReset(const VldContext& vc)
{
  // Check validity.
  if(vc.GetDetector() != Detector::kNear) {
    MAXMSG("Calib",Msg::kWarning,10) << "** WARNING: You are using the QuadLinearityCalScheme " << endl;
    MAXMSG("Calib",Msg::kWarning,10) << "** on non-Near data. This will not correctly calibrate your data." << endl;
    MAXMSG("Calib",Msg::kWarning,10) << "** Try: Calibrator::Instance().Set(\"LinCalibrator=PulserLinearityCalScheme\");" << endl;
  }
  // Requery table.
  fFitTable.NewQuery(vc,kFitType);

  if(fFitTable.GetNumRows()==0) {
    MAXMSG("Calib",Msg::kWarning,10)
      << "No data in CALPULSERFITS task " << kFitType << " vc = " << vc.AsString() << endl;
    IncrementErrors(kLinCalibrator,kMissingTable);
  }

   
}
//......................................................................

---