PmtDriftCalScheme.cxx

Go to the documentation of this file.

// $Id: PmtDriftCalScheme.cxx,v 1.4 2006/04/19 18:50:38 rhatcher Exp $
//
// Calibrator for using even newer PMT-by-PMT drift correction.
// Fast! Compact! Nubile!
//
// n.tagg1@physics.ox.ac.uk
//
//
#include "Calibrator/PmtDriftCalScheme.h"
#include "Calibrator.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexHandle.h"
#include "DatabaseInterface/DbiValidityRec.h"
#include "Validity/VldRange.h"
#include "Conventions/Munits.h"
#include <math.h>

ClassImp(PmtDriftCalScheme)
CVSID("$Id: PmtDriftCalScheme.cxx,v 1.4 2006/04/19 18:50:38 rhatcher Exp $");


//......................................................................
PmtDriftCalScheme::PmtDriftCalScheme()
{
   MSG("Calib",Msg::kVerbose) 
     << "PmtDriftCalScheme::PmtDriftCalScheme" 
     << endl;

   Registry r;
   r.Set("DoPixelCorrection",1);
   r.Set("Mode","PhotonStat");
   r.Set("ReferenceDate",20050101);
   r.Set("ReferenceTime",000010);

   InitializeConfig(r);
}

//......................................................................
void PmtDriftCalScheme::ConfigModified()
{
  const char* str_mode;

  int refdate;
  int reftime;

  bool ok = true;
  ok = ok && GetConfig().Get("DoPixelCorrection",fDoPixelCorrection);
  ok = ok && GetConfig().Get("ReferenceDate",refdate);
  ok = ok && GetConfig().Get("ReferenceTime",reftime);
  ok = ok && GetConfig().Get("Mode",str_mode);
  if(!ok) MSG("Calib",Msg::kWarning) << "PulserDriftCalibrator "
                                     << " Problem when configuring. " <<endl;
  
  // Interpret date.
  fReferenceTime = VldTimeStamp(refdate,reftime,0);

  // Interpret strings.
  if(strcasecmp(str_mode,"PhotonStat")==0) fTask = CalPmtDrift::kPhotonStat;
  else if(strcasecmp(str_mode,"PinDrift")==0) fTask = CalPmtDrift::kPinDrift;
  else {
    fTask = CalPmtDrift::kPhotonStat;
    MSG("Calib",Msg::kError) << "Could not interpret mode " << str_mode 
                             << ". Defaulting to PhotonStat." << std::endl;
  }

  // Ensure the DB has been fixed up.
  Reset(GetContext(),true);
}

//......................................................................
void PmtDriftCalScheme::DoReset( const VldContext& vc )
{

  // Build context for reference point.
  VldContext refContext(vc.GetDetector(),
                        vc.GetSimFlag(),
                        fReferenceTime);


  // Build context for prev, cur, and next drift points.
  // There IS a better way to do this, by querying the DB directly.
  // But that is time consuming and dificult.
  int pointInterval = (int)(3.0 * Munits::hour/Munits::second);;
  if(vc.GetDetector() == Detector::kNear) pointInterval = (int)(1.0*Munits::hour/Munits::second);
  
  VldContext driftVC[3];
  driftVC[0] = VldContext(vc.GetDetector(),
                          vc.GetSimFlag(),
                          vc.GetTimeStamp() - VldTimeStamp(pointInterval,0));
  driftVC[1] = vc;
  driftVC[2] = VldContext(vc.GetDetector(),
                          vc.GetSimFlag(),
                          vc.GetTimeStamp() + VldTimeStamp(pointInterval,0));

  
  // Actually get the tables. Note the order: 
  // it should work well with the L1 cache for
  // forward scanning.
  fRefDrift.NewQuery(refContext,fTask);
  fDrift[0].NewQuery(vc,fTask);
  fDrift[1].NewQuery(vc,fTask);
  fDrift[2].NewQuery(vc,fTask);


  // Error message: no reference rows.
  if(fRefDrift.GetNumRows() == 0) {
    IncrementErrors(kDriftCalibrator,kMissingTable);
    MAXMSG("Calib",Msg::kError,10) 
      << "No rows in CALPMTDRIFT reference table task=" << fTask
      << "  cx=" << refContext.AsString() << std::endl;
  }
  
  // Error message: no current rows.
  if(fDrift[1].GetNumRows() == 0) {
    IncrementErrors(kDriftCalibrator,kMissingTable);
    MAXMSG("Calib",Msg::kWarning,10) 
      << "No rows in CALPMTDRIFT table task=" << fTask
      << "  cx=" << vc.AsString() << std::endl;
  }
  
}


//......................................................................
void PmtDriftCalScheme::PrintConfig( std::ostream& os ) const
{
  os << " PMT Drift Scheme: " << endl;
  os << "  DoPixelCorrection: " << fDoPixelCorrection << endl;
  os << "  ReferenceTime:     " << fReferenceTime.AsString() << endl;
  os << "  Mode:              " << GetConfig().GetCharString("Mode")
     << " (task=" << fTask << ")" << endl;
}


//......................................................................
FloatErr  PmtDriftCalScheme::GetDriftCorrected(FloatErr rawcharge, 
                                           const PlexStripEndId& seid) const
{

  return rawcharge/GetDrift(seid);
}

//......................................................................
FloatErr PmtDriftCalScheme::DecalDrift(FloatErr undrifted, 
                                       const PlexStripEndId& seid) const
{
  return undrifted*GetDrift(seid);
}

//......................................................................
FloatErr  PmtDriftCalScheme::GetDrift(const PlexStripEndId& seid ) const
{

  FloatErr drift1 = GetDrift(seid,1);

  if(drift1 > 0.2 && drift1 < 5.0) // Sanity check!
    return drift1;

  // Ok, now we have a problem. Attempt to ignore this point, and use the surronding ones.

  FloatErr drift0 = GetDrift(seid,0);
  FloatErr drift2 = GetDrift(seid,2);

  FloatErr drift(0,0);
  float n = 0;
  if(drift0 > 0.2 && drift0 < 5.0) { drift+= drift0; n+=1; };
  if(drift2 > 0.2 && drift2 < 5.0) { drift+= drift2; n+=1; };

  if(n==0) {
    IncrementErrors(kDriftCalibrator,kDataInsufficient);
    return FloatErr(1,0.5); // Return 1 with a big error.
  }

  // Use 1 or 2 of the prev/next points, depending on how many are good.
  // Add a systematic error corresponding to the change from d0 to d2.
  drift /= n;
  drift *= FloatErr(1.0, fabs(drift0.GetValue() - drift2.GetValue())*0.5); 

  return drift;
}


//......................................................................
FloatErr  PmtDriftCalScheme::GetDrift(const PlexStripEndId& seid, int whichContext) const
{
  // Find the PMT.
  PlexHandle plex(GetContext());
  PlexPixelSpotId psid = plex.GetPixelSpotId(seid);
  
  UInt_t index = psid.GetUniquePmtEncodedValue();
  
  const CalPmtDrift* drift = fDrift[whichContext].GetRowByIndex(index);
  const CalPmtDrift* ref   = fRefDrift.GetRowByIndex(index);

  if(drift==0) {
    MAXMSG("Calib",Msg::kWarning,10) 
      << "No drift data for strip " << seid.AsString() 
      << " PMT " << psid.AsString("t")
      << " index " << index << endl;
      IncrementErrors(kDriftCalibrator,kMissingRow);
      return FloatErr(0.0,0.0); // Return 0 to trigger error.
  }

  if(ref==0) {
    MAXMSG("Calib",Msg::kError,10) 
      << "No reference drift data for strip " << seid.AsString() << " PMT " << psid.AsString()
      << " index " << index << endl;
      IncrementErrors(kDriftCalibrator,kMissingRow);
      return FloatErr(1.0,0.5); // Add 50% error. Don't attempt interpolation.
  }
  
  FloatErr gain   (drift->GetDrift(),drift->GetStatError());
  gain += FloatErr(0.0,              drift->GetSysError());

  FloatErr gref(ref->GetDrift()); // Set reference point to zero error; let the s2s deal with this.

  if((gain<=0) || (gref<=0)) {
    // One of the gains is faulty, so use the backup plan:
    IncrementErrors(kDriftCalibrator,kDataInsufficient);
    gain = FloatErr(drift->GetCrudeDrift(),drift->GetCrudeDrift()*0.25);
    gref = FloatErr(ref  ->GetCrudeDrift(),  ref->GetCrudeDrift()*0.25);
  }

  FloatErr relativeDrift = gain/gref;

  if(fDoPixelCorrection) relativeDrift *= GetPixelCorrection(psid);
  
  return relativeDrift;
}   


//.........................................................................
FloatErr PmtDriftCalScheme::GetPixelCorrection(const PlexPixelSpotId& /*psid*/) const
{
  //FloatErr peGain;
  //FloatErr width;
  //Calibrator::Instance().DecalGainAndWidth(peGain,width,psid);
  //Apply formula here:
  
  return FloatErr(1.0,0.0); // Default.
}

---