CalVaLinearity.cxx

Go to the documentation of this file.

// $Id: CalVaLinearity.cxx,v 1.9 2007/01/15 19:52:01 rhatcher Exp $
//
// CalVaLinearity
//
// DB Row class corresponding to a single PIN VA liearity correction.
//
// 
//
// Giles Barr g.barr1@physics.ox.ac.uk

#include "Calibrator/CalVaLinearity.h"
#include "MessageService/MsgService.h"
#include "DatabaseInterface/DbiOutRowStream.h"
#include "DatabaseInterface/DbiResultSet.h"
#include "DatabaseInterface/DbiValidityRec.h"
#include <cmath>

ClassImp(CalVaLinearity)

//   Definition of static data members
//   *********************************

CVSID("$Id: CalVaLinearity.cxx,v 1.9 2007/01/15 19:52:01 rhatcher Exp $\n  \
      CVSID_DBIRESULTPTR ");

//  Instantiate associated Result Pointer class.
//  *******************************************

#include "DatabaseInterface/DbiResultPtr.tpl"
template class  DbiResultPtr<CalVaLinearity>;

#include "DatabaseInterface/DbiWriter.tpl"
template class  DbiWriter<CalVaLinearity>;

Double_t CalVaLinearity::Linearize(const Double_t adc) const 
{
  // From raw VA adc value to linearaized VA adc value. For data correction
  Double_t arg[1];
  Double_t par[8];
  
  arg[0] = adc;
  
  for (Int_t i=0; i<6; i++) par[i] = fParam[i];
  par[6] = 0.;   // Internal Pedestal
  par[7] = 1.;   // Internal Gain
  
  // par[0] is the function select (normally 4 or 24)
  
  return CalVaLinearity::FitFunction(arg,par);
}

Double_t CalVaLinearity::UnLinearize(const Double_t linearAdc) const 
{
  // From linearized VA adc to raw adc (i.e. add the saturation in). For MC
  Double_t arg[1];
  Double_t par[8];

  arg[0] = linearAdc;

  for (Int_t i=0; i<6; i++) par[i] = fParam[i];
  par[6] = 0.;   // Internal Pedestal
  par[7] = 1.;   // Internal Gain

  par[0] = par[0] + 1;   // 5 or 25 is the inverse

  return CalVaLinearity::FitFunction(arg,par);
}

void CalVaLinearity::Fill(DbiResultSet& rs, 
                          const DbiValidityRec* /* vrec */) {
//  Purpose:  Fill object from Result Set
//
//  Arguments: 
//    rs           in    Result Set used to fill object
//    vrec         in    Associated validity record (or 0 if filling
//                                                    DbiValidityRec)
//  o Fill object from current row of Result Set.

//  Program Notes
//  Following the template from which this was copied, this method demonstrates 
//  both the "dumb" fill method (just load the data as it comes) and the smart 
//  method (check column name and load according to column order).  The smart method
//  has been 'enhanced' to include a check that all the columns got filled.  It will
//  go kFatal if we don't because we are heading for a divide by zero if some are 
//  incorrect. 

  if ( rs.TableName() == "CALVALINEARITY" ) {
    // Dumb method.
    rs  >> fVAChannel >> fParam[0] 
        >> fParam[1] >> fParam[2] >> fParam[3] >> fParam[4] >> fParam[5];
    fParam[6] = 0.;
    fParam[7] = 1.;
  } else {
    // Smart method
    Int_t numCol = rs.NumCols();
    fVAChannel = 0;

    //  The first column (SeqNo) has already been processed.
    for (Int_t curCol = rs.HasRowCounter() ? 3 : 2; curCol <= numCol; ++curCol) {
      string colName = rs.CurColName();
      if ( colName == "VACHANNEL" )  rs >> fVAChannel; 
      else if( colName == "FSEL"   ) rs >> fParam[0];  
      else if( colName == "PARAM1" ) rs >> fParam[1];  
      else if( colName == "PARAM2" ) rs >> fParam[2];  
      else if( colName == "PARAM3" ) rs >> fParam[3];  
      else if( colName == "PARAM4" ) rs >> fParam[4];  
      else if( colName == "PARAM5" ) rs >> fParam[5];  
      else {
        MSG("Dbi",Msg::kDebug) << "Ignoring column " << curCol 
                               << "(" << colName << ")"
                               << "; not part of CalVaLinearity" 
                               << endl;
        rs.IncrementCurCol();
      } // colname if
    }   // for
  }     // if/else
}
    
//.....................................................................
void CalVaLinearity::Store(DbiOutRowStream& ors,
                               const DbiValidityRec* /* vrec */) const {
  //
  //  Purpose:  Stream object to output row stream
  //
  //  Arguments: 
  //    ors          in     Output row stream.
  //    vrec         in    Associated validity record (or 0 if filling
  
  ors << fVAChannel << fParam[0] 
      << fParam[1] << fParam[2] << fParam[3] << fParam[4] << fParam[5]; 
}

UInt_t CalVaLinearity::Rcid2Index(const RawChannelId& rcid) 
{
  // This static function produces a channel index (used by the database)
  // from a raw channel ID
  return 
    ( rcid.GetDetector() <<28) |
    ( rcid.GetCrate()    <<24) | ( rcid.GetVarcId()   <<20) |
    ( rcid.GetVmm()      <<16) | ( rcid.GetVaAdcSel() <<12) |
    ( rcid.GetVaChip()   << 8) | ( rcid.GetVaChannel()<< 0);
}

RawChannelId CalVaLinearity::Index2Rcid(Int_t index)
{
  // This static function produces a RawChannelId object from a channel index
  // The pedmode, sparsemode, commonmode bits are set to zero (by default)
  UInt_t crate = (index>>24)&0xf;
  UInt_t varc  = (index>>20)&0xf;
  UInt_t vmm   = (index>>16)&0xf;
  UInt_t vfb   = (index>>12)&0xf;
  UInt_t va    = (index>> 8)&0xf;
  UInt_t ch    = (index>> 0)&0xff;

  return RawChannelId(Detector::kFar
                      ,ElecType::kVA
                      ,crate,varc,vmm,vfb,va,ch);
}

Double_t CalVaLinearity::FitFunction(Double_t *arg, Double_t *par) 
{
  // Static function to fit linearity response of PIN channels
  // This is used by the VALinearityConstantMaker package
  // as well as the Calibrator package.

  double x      = *arg;
  double dtype  = par[0];
  int    type   = static_cast<int>(dtype+0.1);
  double k1     = par[1];
  double k2     = par[2];
  double k3     = par[3];
  double k4     = par[4];
  //  double k5     = par[5];
  double ped    = par[6];
  double gain   = par[7];

  // **************************************************************************
  // **** Pole with adjustable power and one variable high end pole kicker ****
  // **************************************************************************

  if (x < 0.1) return x;    // Do not linearise negative values

  if (type == 24) {
    double x1 = x;
    x1 = ped + x1 * gain ;
    x1 = Pole(x1,k1,k2,-1.);
    if ((x1-k3) > 0) {      // High end kicker function
      x1 = (x1 - k3)/(1. - (x1-k3)/k4) + k3;
    }
    return x1;
  }
  else if (type == 25) {    // Non-inverted function
    double x1 = x;
    if ((x1-k3) > 0) {      // High end kicker function
      x1 = (x1 - k3)/(1. + (x1-k3)/k4) + k3;
    }
    x1 = Pole(x1,k1,k2,+1.);
    x1 = (x1-ped)/gain ;    
    return x1;
  }

  // **************************************************************************
  // Pole with adjustable power, no kicker
  // **************************************************************************

  else if (type == 4) {
    double x1 = x;
    x1 = ped + x1 * gain ;
    x1 = Pole(x1,k1,k2,-1.);
    return x1;
  }
  else if (type == 5) {     // Non-inverted function
    double x1 = x;
    x1 = Pole(x1,k1,k2,+1.);
    x1 = (x1-ped)/gain ;    
    return x1;
    //    printf("f %g %g %g %g\n",result,x,k1,k2);
  }

  // **************************************************************************
  // Straight line conversion (used for fitting, not for linearization)
  // **************************************************************************

  else if (type == 0) {
    return ped + x * gain ;
  } 
  else if (type == 1) {
    return (x-ped)/gain;
  }

  // **************************************************************************
  // Remedial functions
  // **************************************************************************

  else if (type == 10) {
    return x;    // Straight line y = x function
  } 
  else if (type == 12) {
    return 0.;   // Dead channel function
  }
  else if (type == 14) {
    return k1;   // Saturated on main Pole
  } 
  else if (type == 16) {
    return 1.e7; // Infinity
  } 

  return x;   // Default (last resort) is to return y=x function
}

Double_t CalVaLinearity::Pole(Double_t x, Double_t k1, 
                              Double_t k2, Double_t sign) 
{
  // Local function to do the pole calculation.  sign is either +1 or -1
  // depending on wheter the function is in use for the inverse or not.
  
  // The function implemented is
  //      x1 = x / (1 + (x1/k1)**k2 ) ** (1./k2)
  
  Double_t x1 = x;
  if (k1 <= 500.) {   // Function is flatline [was if (k1<=0.)]
    x1 = 0.;
  } else if (sign <= 0. && k1 <= x1) { // Above pole: make function 'flatline'
    x1 = 0.;
  } else if (x1 < 0.1) {   // Don't modify (x=y when x is negative)
  } else if (k2 <= 0.1) {  // don't modify with negative power [was if (k2<=0.)]
  } else {
    Double_t a = k2*log(x1/k1);   // Not negative log, both x1 and k1 checked
    if (a > 228.) {
      if (sign < 0.) {   // Only happens when sign positive (print if not) 
        //printf("pole:s x1 %g k1 %g k2 %g sign %g\n",x1,k1,k2,sign);
        x1 = k1;
      }
    } else {  // Should not overflow exp()
      Double_t b = pow(1. + sign * exp(a), 1./k2);
      if (b > 1.0e-20) {   // OK
        x1 = x1 / b;
      } else {  // Very close to pole 
        x1 = x1*1.0e20;
      }
    }
  }
  return x1;
}

---