PulserTimeCalScheme.cxx

Go to the documentation of this file.

// $Id: PulserTimeCalScheme.cxx,v 1.23 2008/10/31 14:36:12 tagg Exp $
//
// Calibrator scheme for doing both LI hardware correction
// and muon correction.
//
// Inherits from the plain-old TimeCalScheme.
//
// n.tagg1@physics.ox.ac.uk

#include "PulserTimeCalScheme.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexHandle.h"
#include "Plex/PlexSEIdAltL.h"
#include "Plex/PlexStripEndId.h"
#include "DatabaseInterface/DbiValidityRec.h"
#include "Validity/VldRange.h"
#include "Calibrator.h"

CVSID("$Id: PulserTimeCalScheme.cxx,v 1.23 2008/10/31 14:36:12 tagg Exp $");
ClassImp(PulserTimeCalScheme)


const RawChannelId kVaRcid(Detector::kFar,ElecType::kVA,0,0);

PulserTimeCalScheme::PulserTimeCalScheme()
{
  Registry r;
  // Set configurable options.
  r.Set("UseMuonCalibration",1);
  r.Set("UsePulserCalibration",0);
  r.Set("UseShieldCalibration",1);
  r.Set("UseJumps",1);
  r.Set("JumpTask",1); // (Obsolete - now tied to MuonTask)
  r.Set("DoWalkCorrection",1);
  r.Set("PulserLed1",2);
  r.Set("PulserLed2",3);
  r.Set("PulserLed3",4);
  r.Set("DefaultTimingCardSetting",260.);
  //r.Set("OverallOffset",30.*Munits::ns);
  r.Set("ReferenceDate",20040930);
  r.Set("ReferenceTime",0);
  r.Set("MuonTask","Dogwood");
  InitializeConfig(r); 
}

//.....................................................................
void PulserTimeCalScheme::ConfigModified()
{                                       
  Int_t refDate = 0;
  Int_t refTime = 0;;
  
  const char* taskname;

  bool ok = true;
  ok = ok && GetConfig().Get("UseMuonCalibration",fUseMuonCalibration);
  ok = ok && GetConfig().Get("UseShieldCalibration",fUseShieldCalibration);
  ok = ok && GetConfig().Get("UsePulserCalibration",fUsePulserCalibration);
  ok = ok && GetConfig().Get("UseJumps",fUseJumps);
  //ok = ok && GetConfig().Get("JumpTask",fJumpTask);
  ok = ok && GetConfig().Get("DoWalkCorrection",fDoWalkCorrection);
  ok = ok && GetConfig().Get("DefaultTimingCardSetting",fDefaultTimingCardSetting);
  //ok = ok && GetConfig().Get("OverallOffset",fOverallOffset);
  ok = ok && GetConfig().Get("ReferenceDate",refDate);
  ok = ok && GetConfig().Get("ReferenceTime",refTime);
  for(int i=0;i<kMaxLeds; i++) {
    ok = ok && GetConfig().Get(Form("PulserLed%d",i+1),fPulserLed[i]);
  }
  ok = ok && GetConfig().Get("MuonTask",taskname);

  if(!ok){
    MSG("Calib",Msg::kWarning) << "PulserTimeCalScheme couldn't configure properly!" << endl;
    GetConfig().PrettyPrint(std::cerr);
  }

  if     (0==strcasecmp(taskname,"Dogwood")){ fMuonTask = 2; }
  else if(0==strcasecmp(taskname,"Andys"  )){ fMuonTask = 1; }
  else if(0==strcasecmp(taskname,"Brians" )){ fMuonTask = 0; }
  else{
    // Feh?
    MSG("Calib",Msg::kWarning) << "PulserTimeCalScheme doesn't recognise task " << taskname << endl;
  }  

   fReferenceTime = VldTimeStamp(refDate,refTime,0);

  // Ensure that the DB has been changed for this event.
  Reset(GetContext(),true); 
}

//.....................................................................
void PulserTimeCalScheme::PrintConfig( std::ostream& os ) const
{
  char taskname[100];
  switch(fMuonTask) {
  case 0: sprintf(taskname,"Brians"); break;
  case 1: sprintf(taskname,"Andys"); break;
  case 2: sprintf(taskname,"Dogwood"); break;
  default: sprintf(taskname,"Task %d",fMuonTask); break;
  }

  os << "  UseMuonCalibration:   " << ((fUseMuonCalibration)?("On"):("Off")) << endl;
  os << "  MuonTask:             " << taskname << endl;
  os << "  UseShieldCalibration: " << ((fUseShieldCalibration)?("On"):("Off")) << endl;
  os << "  UseJumps:             " << ((fUseJumps)?("On"):("Off")) << endl;
  // os << "  JumpTask:             " << fJumpTask << endl;
  os << "  UsePulserCalibration: " << ((fUsePulserCalibration)?("On"):("Off")) << endl;  
  os << "  DoWalkCorrection:     " << ((fDoWalkCorrection)?("On"):("Off")) << endl;  
  os << "  Reference date:       " << fReferenceTime.AsString() << endl;
  os << "  PulserLeds to use:    ";
  for(int i=0;i<kMaxLeds;i++) {
    os << fPulserLed[i];
    if(i!=kMaxLeds-1) os << ",";
  }
  os << endl;
    
}

 //.....................................................................
void PulserTimeCalScheme::DoReset( const VldContext& context )
{
  MSG("Calib",Msg::kVerbose) << "PulserTimeCalScheme::DoReset()" << endl;
  
  if(fUseMuonCalibration) {
    fMuonResPtr.NewQuery(context,fMuonTask);
    
    if(fMuonResPtr.GetNumRows()==0) {
      MAXMSG("Calib",Msg::kWarning,10) 
        << "No rows in CALTIMECALIBRATION with validity context "
        << context.AsString() << "  No muon timing calibration will be applied." << endl;
      IncrementErrors(kTimeCalibrator,kMissingTable);
    }
  }
  
  if(fUseShieldCalibration && context.GetDetector()==Detector::kFar) {
    fShieldResPtr.NewQuery(context,99); // Task 99 is for shield.
    if(fShieldResPtr.GetNumRows()==0) {
        MAXMSG("Calib",Msg::kWarning,10) 
          << "No rows for FD shield timing calibration CALTIMECALIBRATION task " 
          << 99 << " with validity context "
          << context.AsString() << endl;
        IncrementErrors(kTimeCalibrator,kMissingTable);
    }
  }

  if(fUseJumps && fMuonTask>0) {
    fJumpResPtr.NewQuery(context,fMuonTask);
    // No errors if no data.
  }

  VldContext refContext(context.GetDetector(),
                        context.GetSimFlag(),
                        fReferenceTime);

  if(fUsePulserCalibration) {
    // Note: we use the LED number to set the Task in the DB call.
    // This allows us to use different LEDs for the calibration if we desire.

    for(int iled=0;iled<kMaxLeds;iled++) {
      fDriftPtr[iled].NewQuery(context,fPulserLed[iled]);
      if(fDriftPtr[iled].GetNumRows()==0) {
        MAXMSG("Calib",Msg::kWarning,10) 
          << "No rows in PULSERTIMEDRIFT task " << fPulserLed[iled] << " with validity context "
          << context.AsString() << endl;
        IncrementErrors(kTimeCalibrator,kMissingTable);
      }

      
      fDriftRefPtr[iled].NewQuery(refContext,fPulserLed[iled]);
      if(fDriftRefPtr[iled].GetNumRows()==0) {
        MAXMSG("Calib",Msg::kWarning,10) 
          << "No rows in PULSERTIMEDRIFT task " << fPulserLed[iled] << " with reference context "
          << refContext.AsString() << endl;
        IncrementErrors(kTimeCalibrator,kMissingTable);
      }
      
    }

    // Time card and reference.

    fTimingCardResPtr.NewQuery(context,0);
    if(fTimingCardResPtr.GetNumRows()==0) {
      MAXMSG("Calib",Msg::kWarning,10) 
        << "No rows in PULSERTIMINGCARDSETTING with validity context "
        << context.AsString() << endl;
      IncrementErrors(kTimeCalibrator,kMissingTable);
    }

    fTimingCardReferencePtr.NewQuery(refContext,0);
    if(fTimingCardReferencePtr.GetNumRows()==0) {
      MAXMSG("Calib",Msg::kWarning,10) 
        << "No rows in PULSERTIMINGCARDSETTING with reference context "
        << refContext.AsString() << endl;
      IncrementErrors(kTimeCalibrator,kMissingTable);
    }

  }
}

//.....................................................................
DoubleErr PulserTimeCalScheme::GetCalibratedTime(DoubleErr rawtime,
                                                FloatErr rawcharge,
                                                const PlexStripEndId& seid) const
{

  // Find raw channel.
  DoubleErr time = rawtime;

  if(seid.IsVetoShield()) {

    // Shield hits:
    if(fUseShieldCalibration) time = CalibrateShield(time,rawcharge,seid,true);     
  } else {

    // Non-shield hits:

    if(GetContext().GetDetector()==Detector::kFar) {
      PlexHandle plex(GetContext());
      RawChannelId rcid = plex.GetRawChannelId(seid);
      if(fUseJumps)             time = CalibrateByJumps(time,rcid,true);
      if(fUsePulserCalibration) time = CalibrateByPulser(time,rcid,true);
    }
    if(fUseMuonCalibration)   time = CalibrateByMuon(time,seid,true);
    if(fDoWalkCorrection)     time-= WalkCorrection(rawcharge,seid);
  }

  return time;
}

//.....................................................................
DoubleErr PulserTimeCalScheme::DecalTime(DoubleErr caltime, 
                                        FloatErr rawcharge,
                                        const PlexStripEndId& seid) const
{

  // Find raw channel.
  DoubleErr time = caltime;

  if(seid.IsVetoShield()) {

    // Shield hits:
    if(fUseShieldCalibration) time = CalibrateShield(time,rawcharge,seid,false); 

  } else {

    // Non-shield hits:

    if(GetContext().GetDetector()==Detector::kFar) {
      PlexHandle plex(GetContext());
      RawChannelId rcid = plex.GetRawChannelId(seid);
      if(fUseJumps)             time = CalibrateByJumps(time,rcid,false);
      if(fUsePulserCalibration) time = CalibrateByPulser(time,rcid,false);
    }
    if(fUseMuonCalibration)   time = CalibrateByMuon(time,seid,false);
    if(fDoWalkCorrection)     time+= WalkCorrection(rawcharge,seid);
  }

  return time;
}


//.....................................................................
DoubleErr PulserTimeCalScheme::CalibrateByJumps(DoubleErr rawtime,               
                                                const RawChannelId& rcid,
                                                Bool_t calibrate) const
{
  // Find index number into table.
  UInt_t index = CalTimeJump::GetIndex(rcid);

  // If a jump entry exists, apply it.

  const CalTimeJump* jump = fJumpResPtr.GetRowByIndex(index);
  if(jump) {
    if(calibrate)
      return rawtime - (double)jump->GetJump();
    else
      return rawtime + (double)jump->GetJump();
  }

  // Otherwise, do nothing.

  return rawtime;
}

//.....................................................................
DoubleErr PulserTimeCalScheme::CalibrateByMuon(DoubleErr rawtime,               
                                           const PlexStripEndId& seid,
                                           Bool_t calibrate) const
{
  // Now need to get the row which corresponds to the stripendnum.
  const CalTimeCalibration* dpgc = fMuonResPtr.GetRowByIndex(seid.BuildPlnStripEndKey());

  if(dpgc==0) {
    if(fMuonResPtr.GetNumRows()>0) {
      MAXMSG("Calib",Msg::kWarning,10) 
        << "PulserTimeCalScheme: No database row for StripEnd " << seid.AsString() << "\n";
      
      IncrementErrors(kTimeCalibrator,kMissingRow,seid);
    }
    return rawtime + DoubleErr(0,20.*Munits::ns); // Show there's no calib constant.
  }
  if(calibrate)
    return (rawtime/(double)dpgc->GetScale()) -  (double)dpgc->GetOffset();

  // else decalibrate:
  // No inverse function exists in the row function, so I'll kludge one here.
  // NJT 7/04
  return (rawtime +  (double)dpgc->GetOffset())*(double)dpgc->GetScale();
}



//.....................................................................
DoubleErr PulserTimeCalScheme::CalibrateByPulser(DoubleErr rawtime,               
                                             const RawChannelId& rcid,
                                             Bool_t calibrate) const
{
  // Find index number into table.
  UInt_t index = PulserTimeDrift::GetIndex(rcid);
  
  const PulserTimeDrift* drift    = 0;
  const PulserTimeDrift* driftref = 0;

  // Get drift point and drift reference point
  Int_t whichled = 0;
  bool driftok = GetPulserCalibration(whichled,index,drift,driftref);
  
  // If not good, look through the list of backups.
  while(!driftok){
    whichled++;
    if(whichled<kMaxLeds) {
      driftok = GetPulserCalibration(whichled,index,drift,driftref);
      MSG("Calib",Msg::kDebug) << "PulserTiming: " << rcid.AsString() << " Fell back to LED " << fPulserLed[whichled] << endl;
    } else {
      MAXMSG("Calib",Msg::kWarning,10) 
        << "PulserTimeCalScheme: Failed to find calibration for channel " 
        << rcid.AsString() << " index " << index << " in primary or fallback tables.\n";
      IncrementErrors(kTimeCalibrator,kMissingRow,rcid);
      return rawtime + DoubleErr(0,50.*Munits::ns); // no calibration. Add 50ns error to show.
    }
  }
  
  
  // Get Settings row, if it exists.
  const PulserTimingCardSetting* cardSetting = 
    fTimingCardResPtr.GetRowByIndex(rcid.GetCrate());
  const PulserTimingCardSetting* refCardSetting = 
    fTimingCardReferencePtr.GetRowByIndex(rcid.GetCrate());

  Double_t tpmt_delay_ns = fDefaultTimingCardSetting;
  Double_t ref_delay_ns  = fDefaultTimingCardSetting;
  
  // Check for errors getting setting.
  if(cardSetting==0) {
    if(fTimingCardResPtr.GetNumRows()>0) {
      MAXMSG("Calib",Msg::kWarning,10) << "Incomplete PulserTimingCardSettings table!" << endl;      
    }
    IncrementErrors(kTimeCalibrator,kMissingRow,rcid);    
  } else {
    tpmt_delay_ns = cardSetting->GetDelayNanoSecs();
  }
  
  
  // Check for errors getting reference setting.
  if(refCardSetting==0) {
    if(fTimingCardReferencePtr.GetNumRows()>0) {
      MAXMSG("Calib",Msg::kWarning,10) << "Incomplete PulserTimingCardSettings table!" << endl;      
    }
    IncrementErrors(kTimeCalibrator,kMissingRow,rcid);    
  } else {
    ref_delay_ns = cardSetting->GetDelayNanoSecs();
  }
  
  
  // Finally, do the correction:
  
  // Crate offset due to changed timing card setting:
  double crate_offset_t = (tpmt_delay_ns - ref_delay_ns)*Munits::ns;
  
  // Offset as calibrated by LI. Longer timing card settings
  // mean that the drift point is too big a number.
  double drift_tdcs = (drift->GetMean() - driftref->GetMean());
  double drift_err  = (drift->GetRms() / sqrt(drift->GetEntries()));
  double convert = Calibrator::Instance().GetTDCConvert(kVaRcid); // s per tick
  DoubleErr drift_t(drift_tdcs * convert,
                    drift_err  * convert);
  
  // And the time walk.
  
  
  DoubleErr offset = drift_t - crate_offset_t;
  
  if(calibrate) 
    return rawtime - offset;
  
  // else:
  return rawtime + offset;
}

Bool_t 
PulserTimeCalScheme::GetPulserCalibration(Int_t iled, UInt_t index,
                                          const PulserTimeDrift* &drift, 
                                          const PulserTimeDrift* &ref) const
{
  // Retrieve rows from the DB. Return true if the rows exist and
  // are good, false if they are bad.

  drift = fDriftPtr[iled].GetRowByIndex(index);
  ref   = fDriftRefPtr[iled].GetRowByIndex(index);

  // Do they exist?
  if(drift==0) return false;
  if(ref==0) return false;

  // Are the means valid to a good precision?
  // Error limit: 0.2 ticks.
  // So, want to fail if RMS/Entries > 0.2 ticks
  // Rearrange to get below:
  // (The clever bit here is that this automatically checks for entries==0)

  if( drift->GetRms() >= sqrt(drift->GetEntries())*0.2 ) return false;
  if(   ref->GetRms() >= sqrt(  ref->GetEntries())*0.2 ) return false;


  return true;
}



//.....................................................................
DoubleErr PulserTimeCalScheme::WalkCorrection(FloatErr rawcharge, 
                                              const PlexStripEndId&) const
{

  // These are Andy's paramaterisations for the FD, hard-coded for now.
  if(GetContext().GetDetector()==Detector::kFar) {

    const float kParData[4] = {  5.44066e+00,
                              -5.94887e-01,
                              -7.36609e-02,
                              7.66637e-03 };

    const float kParMC[4] = { 4.18608e+00,
                              -3.75033e-01,
                              -9.50218e-02,
                              9.52560e-03 };
    
    const float* par = kParData;
    if(GetContext().GetSimFlag()==SimFlag::kMC) par = kParMC;
      
    if(rawcharge<=0) rawcharge = 1; // Be careful.
    double cx = log(rawcharge()/2.3);
    double walk = (
                   par[0]
                   + (par[1] * cx)
                   + (par[2] * cx *cx)
                   + (par[3] * cx * cx *cx) 
                   );
    // These parameters are in units of meters, so convert to time:
    return walk / Munits::c_light; // Fixme: does not paramaterize error.
      
  }

  return 0;
}

//.....................................................................
DoubleErr PulserTimeCalScheme::CalibrateShield(DoubleErr rawtime,
                                               FloatErr rawcharge,
                                               const PlexStripEndId& seid, 
                                               Bool_t calibrate) const
{
   // Now need to get the row which corresponds to the stripendnum.
  const CalTimeCalibration* dpgc = fShieldResPtr.GetRowByIndex(seid.BuildPlnStripEndKey());

  if(dpgc==0) {
    if(fShieldResPtr.GetNumRows()>0) {
      MAXMSG("Calib",Msg::kWarning,10) 
        << "PulserTimeCalScheme: No database row for shield StripEnd " << seid.AsString() << "\n";      
      IncrementErrors(kTimeCalibrator,kMissingRow,seid);
    }
    return rawtime - DoubleErr(-30.0*Munits::ns, // -30ns is the usual average
                               20.*Munits::ns); // Show there's no calib constant.
  }

  // apply the 'extra' walk correction for shield hits.
  if(rawcharge<1) rawcharge=1; // Protect against -ve numbers
  Double_t logq = log(rawcharge);

  Double_t walk = 0;
  if(fDoWalkCorrection) {
    const float kgainadjust=60./90.; // I think this compensates for the extra gain in the shield
    walk = WalkCorrection(rawcharge*kgainadjust,seid);
    walk += dpgc->GetSlewC1()
    + dpgc->GetSlewC2()*logq
    + dpgc->GetSlewC3()*logq*logq
    + dpgc->GetSlewC4()*logq*logq*logq;
  }


  if(calibrate)
    return rawtime - walk + (double)dpgc->GetOffset();
  
  // else decalibrate:
  return  rawtime +  walk - (double)dpgc->GetScale();

}

---