PdbMutCalibStrip.cc
//-----------------------------------------------------------------------------
//
// The pdbcal package
// Copyright (C) PHENIX collaboration, 2002
//
// Declaration of class PdbMutCalibStrip
//
// Purpose: Store MuTr calibration information; one object, one strip
//
// Description:
//
// Author: silvermy (silvermy@lanl.gov)
//-----------------------------------------------------------------------------
#include <iomanip>
#include <cstdio>
#include "PdbMutCalibStrip.hh"
PdbMutCalibStrip::PdbMutCalibStrip()
{
set_indices(0, 0, 0, 0, 0, 0, 0); // init to dummy values
zero();
}
PdbMutCalibStrip::PdbMutCalibStrip(int iarm, int istation, int ioctant,
int ihalfoctant, int igap, int iplane,
int istrip)
{
set_indices(iarm, istation, ioctant, ihalfoctant, igap, iplane, istrip);
zero();
}
void
PdbMutCalibStrip::zero()
{
pedestal = 0;
gain = 0;
rms = 0;
for (int i=0; i<NMUTCALIBPAR; i++)
{
calibpar[i] = 0;
}
Status = 0;
Step = 0;
Saturation = 0;
}
PdbMutCalibStrip::~PdbMutCalibStrip()
{
}
PdbMutCalibStrip&
PdbMutCalibStrip::operator=(const PdbMutCalibStrip& rhs)
{
ArmNum = rhs.ArmNum;
StationNum = rhs.StationNum;
OctantNum = rhs.OctantNum;
HalfOctantNum = rhs.HalfOctantNum;
GapNum = rhs.GapNum;
PlaneNum = rhs.PlaneNum;
StripNum = rhs.StripNum;
pedestal = rhs.pedestal;
gain = rhs.gain;
rms = rhs.rms;
for (int i=0; i<NMUTCALIBPAR; i++)
{
calibpar[i] = rhs.calibpar[i];
}
Status = rhs.Status;
Step = rhs.Step;
Saturation = rhs.Saturation;
return *this;
}
void
PdbMutCalibStrip::print() const
{
write(cout);
}
void
PdbMutCalibStrip::write(ostream& os) const
{
os << setiosflags(ios::fixed)
<< setw(1); // one digit only for the first numbers
os << ArmNum << " "
<< StationNum << " "
<< OctantNum << " "
<< HalfOctantNum << " "
<< GapNum << " "
<< PlaneNum<< " ";
os << setw(3) << StripNum << " ";
os << setw(8) << setprecision(3)
<< pedestal << " "
<< setw(6) << gain << " "
<< setw(5) << rms << " " << setw(7);
for (int i=0; i<NMUTCALIBPAR; i++)
{
os << calibpar[i] << " " << setw(9) << setprecision(6);
}
os << setw(4) << Status << " "
<< setw(1) << Step << " "
<< setw(3) << Saturation;
os << endl;
}
void
PdbMutCalibStrip::read(istream& is)
{
// opposite of write method
is >> ArmNum
>> StationNum
>> OctantNum
>> HalfOctantNum
>> GapNum
>> PlaneNum;
is >> StripNum
>> pedestal
>> gain
>> rms;
for (int i=0; i<NMUTCALIBPAR; i++)
{
is >> calibpar[i];
}
is >> Status
>> Step
>> Saturation;
}
void
PdbMutCalibStrip::set_indices (int iarm, int istation, int ioctant,
int ihalfoctant, int igap, int iplane,
int istrip)
{
ArmNum = iarm;
StationNum = istation;
OctantNum = ioctant;
HalfOctantNum = ihalfoctant;
GapNum = igap;
PlaneNum = iplane;
StripNum = istrip;
}
void
PdbMutCalibStrip::set_values (float fped, float fgain, float frms)
{
pedestal = fped;
gain = fgain;
rms = frms;
}
void
PdbMutCalibStrip::set_calibparAll (float temp[NMUTCALIBPAR])
{
for (int i=0; i<NMUTCALIBPAR; i++)
{
calibpar[i] = temp[i];
}
}
void
PdbMutCalibStrip::get_indices (int & iarm, int & istation, int & ioctant,
int & ihalfoctant, int & igap, int & iplane,
int & istrip) const
{
iarm = ArmNum;
istation = StationNum;
ioctant = OctantNum;
ihalfoctant = HalfOctantNum;
igap = GapNum;
iplane = PlaneNum;
istrip = StripNum;
}
void
PdbMutCalibStrip::get_values (float & fped, float & fgain, float & frms) const
{
fped = pedestal;
fgain = gain;
frms = rms;
}
void
PdbMutCalibStrip::get_calibparAll (float temp[NMUTCALIBPAR]) const
{
for (int i=0; i<NMUTCALIBPAR; i++)
{
temp[i] = calibpar[i];
}
}
int
PdbMutCalibStrip::getUniqueId () const
{
int id = StripNum;
// arm: 0,1; bit 17
// station: 0,1,2; bits 15,16
// octant: 0-7; bits 12-14
// hoct: 0,1; bit 11
// gap: 0,1,2; bit 9,10
// plane: 0,1; bit 8
// strip: 0..160 (max on North); bits 0-7
id += PlaneNum << 8;
id += GapNum << 9;
id += HalfOctantNum << 11;
id += OctantNum << 12;
id += StationNum << 15;
id += ArmNum << 17;
return id;
}
bool
PdbMutCalibStrip::getName(char *name =
"muTr_ArmX_StaY_OctZ_HOctW_GapI_PlaneJ_StripK") const
{
char namestring[50];
sprintf(namestring,"muTr_Arm%d_Sta%d_Oct%d_HOct%d_Gap%d_Plane%d_Strip%d",
ArmNum,StationNum,OctantNum,HalfOctantNum,GapNum,PlaneNum, StripNum);
strcpy(name,namestring);
return true;
}
int
PdbMutCalibStrip::getAdc(int dac) const
{
// check input somewhat
if ( dac<0 || dac>MAXMUTRDAC )
{
cerr << "PdbMutCalibStrip::getAdc - DAC value " << dac
<< " is out of bounds " << endl;
return -999;
}
float convped = MAXMUTRADC - pedestal; // pedestal is a value around MAXMUTRADC
// therefore conv(enient)ped is a value around 0.
// use values that increase with the signal throughout routine
// and return reversed value
float adcval = convped + gain*dac; // default is just a straight line
float center = calibpar[0];
float distCenter, distSat;
// Make an 'integration' of how much the non-linear deviation has cost you
// up to this point
int nsteps = dac/Step;
float sumdev = 0;
float dacval, olddiff = 0;
for (int i = 0; i<nsteps; i++)
{
dacval = i*Step;
distCenter = dac - center;
if (distCenter<0) // low side
{
sumdev = sumdev + calibpar[1]*distCenter*distCenter;
}
else
{
if ( (dacval - Saturation) < 0)
{ // not yet saturated: normal
olddiff = calibpar[2]*distCenter*distCenter;
sumdev += olddiff;
}
else
{// we have passed saturation point, new rules apply here
// olddiff is kept from previous step
distSat = dacval - Saturation;
sumdev = sumdev + olddiff + calibpar[3]*distSat;
}
}
}
adcval -= sumdev;
// we are done. Now reverse it back to the real scale
adcval = MAXMUTRADC - adcval;
int adc = (int) (adcval+0.5); // return nearest integer
return adc;
}
int
PdbMutCalibStrip::getDac(int adc) const
{
const int MAXITERATIONS = 10;
// just do the things done in getAdc backwards..
float adcval = MAXMUTRADC - (float) adc; // reversed to convenient scale
float dacval = (adcval - pedestal)/gain; // initial guess
int dac = (int) (dacval+0.5);
if (dac > MAXMUTRDAC)
{
return MAXMUTRDAC; // this is the maximum DAC value
}
float center = calibpar[0];
float distCenter, distSat;
// Make an 'integration' of how much the non-linear deviation has cost you
// up to this point
int nsteps = dac/Step;
float sumdev = 0;
float dacstepval, olddiff = 0;
for (int i = 0; i<nsteps; i++)
{
dacstepval = i*Step;
distCenter = dacstepval - center;
if (distCenter<0) // low side
{
sumdev = sumdev + calibpar[1]*distCenter*distCenter;
}
else
{
if ( (dacstepval - Saturation) < 0)
{ // not yet saturated: normal
olddiff = calibpar[2]*distCenter*distCenter;
sumdev += olddiff;
}
else
{// we have passed saturation point, new rules apply here
// olddiff is kept from previous step
distSat = dacstepval - Saturation;
sumdev = sumdev + olddiff + calibpar[3]*distSat;
}
}
}
// Some minor adjustments can be needed: maybe off-by-one or a few steps
dacval += sumdev/gain;
int extrastep = (int) (dacval/Step) - nsteps;
int diffstep, sign, iteration=0;
while ( extrastep!=0 && iteration<MAXITERATIONS )
{ // we didn't go far enough or too far, avoid an infinite loop also
sumdev = 0;
if (extrastep>0)
{
sign = 1; // go forward
}
else
{
extrastep = - extrastep;
sign = -1; // back up
}
for (int i = 0; i<extrastep; i++)
{
dacstepval = (nsteps+i*sign)*Step;
distCenter = dacstepval - center;
if (distCenter<0) // low side
{
sumdev = sumdev + sign*calibpar[1]*distCenter*distCenter;
}
else
{
if ( (dacstepval - Saturation) < 0)
{ // not yet saturated: normal
olddiff = calibpar[2]*distCenter*distCenter;
sumdev += olddiff;
}
else
{// we have passed saturation point, new rules apply here
// olddiff is kept from previous step
distSat = dacstepval - Saturation;
sumdev = sumdev + olddiff + calibpar[3]*distSat;
}
}
}
dacval += sumdev/gain;
diffstep = sign*Step;
nsteps += diffstep;
extrastep = (int) (dac/Step) - nsteps;
if ( extrastep!=0 )
{ // continued iteration: maybe not necessary?
if ( (extrastep+diffstep) == 0) //
{ // You just want to go back where you just came from.
extrastep = 0;
}
if (extrastep>2)
{ // something wrong - trying too many steps
// possibly a bad channel or bad fit
cerr << " PdbMutCalibStrip::getDac - stopped iteration - too many steps" << endl;
iteration = MAXITERATIONS - 1; // abort, back to initial guess
}
}
iteration++;
}
if (iteration == MAXITERATIONS)
{ // did not converge! we'll have to settle for the initial guess
dacval = (adcval - pedestal)/gain;
}
dac = (int) (dacval+0.5);
return dac;
}