MuonDriftCalScheme.cxx

Go to the documentation of this file.

// $Id: MuonDriftCalScheme.cxx,v 1.13 2009/05/22 10:59:42 tinti Exp $
//
// Calibrator for detector-wide drift by muon measurements.
//
// Justin Evans, Nathaniel Tagg
// j.evans2@physics.ox.ac.uk, tagg@minos.phy.tufts.edu
//
// $Log: MuonDriftCalScheme.cxx,v $
// Revision 1.13  2009/05/22 10:59:42  tinti
// Get rid of duplicate code. Now all reco version factors are stored in a single function in CalDrift. To get the right constants it's enough to use CalDrift::GetMeadianAsR1_18_2(..) - the same for the associate error.
//
// Revision 1.12  2009/05/21 18:24:43  tinti
// reco factor for Dogwood. Calibration drift constant is FrontMedian() in the Near detector if it is Dogwood. Change the interpolator to apply reco scale factor consistently
//
// Revision 1.11  2007/04/01 14:00:20  evans
// Adding the Cedar reco version factor for the far detector.
//
// Revision 1.10  2007/03/21 12:07:55  evans
// Ensuring the version shear factors are detector-specific.
//
// Revision 1.9  2007/03/16 17:59:53  evans
// Putting in a version shear correction factor for all reconstruction
// versions. Cedar is not yet corrected for (waiting for a batch job to run).
//
// All versions are corrected to be consistent with R1_18_2.
//
// Revision 1.8  2006/04/19 18:50:38  rhatcher
// Change DetectorType:: to Detector:: everywhere.
//
// Revision 1.7  2006/04/09 16:48:39  evansj
// Removing another cout debugging statement.
//
// Revision 1.6  2006/04/09 16:46:12  evansj
// Removing a cout debugging statement.
//
// Revision 1.5  2006/04/09 16:38:10  evansj
// Updating the DoReset() method to include the version shear factor.
// Now all drift correction factors are corrected to be consistent with R1_18 reconstruction.
//
// Revision 1.4  2006/04/07 22:01:49  tagg
// Add <math.h> to make gcc 4.0 happy.
//
// Revision 1.3  2006/04/07 16:58:33  evansj
// Adding a RecoVersion flag to CalDrift.
// Minor change to a warning message in MuonDriftCalScheme.
// Making CalDrift copy constructor public.
//
// Revision 1.2  2006/02/27 20:43:29  tagg
// Dangit! Erased my last CVS commit message.
// Change the default reference date to Dec 01, 2005.
//
// Revision 1.1  2006/02/27 20:40:50  tagg
// *** empty log message ***
//
//
#include <math.h>
#include <string>

#include "Calibrator/MuonDriftCalScheme.h"
#include "Calibrator.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexHandle.h"
#include "DatabaseInterface/DbiValidityRec.h"
#include "DatabaseInterface/DbiSqlContext.h"
#include "Validity/VldRange.h"
#include "Conventions/Munits.h"

using std::string;

ClassImp(MuonDriftCalScheme)
CVSID("$Id: MuonDriftCalScheme.cxx,v 1.13 2009/05/22 10:59:42 tinti Exp $");


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

   Registry r;
   //task default is still 0, meaning pre Dogwood database
   //from Dogwood1 on the task number needs to be set to 1
   r.Set("Task",0);
   r.Set("ReferenceDate",20051201);
   r.Set("ReferenceTime",000010);

   InitializeConfig(r);
}

//......................................................................
void MuonDriftCalScheme::ConfigModified()
{
  int refdate;
  int reftime;

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

  // Ensure the DB has been fixed up.
  Reset(GetContext(),true);
}
//......................................................................
void MuonDriftCalScheme::DoReset( const VldContext& vc )
{
 

  fCurError = 0;
  fCurDrift.Set(1.0 ,0.5); // Default: No drift, 50% error.
  string refRecoVersion = "";
  string recoVersion = "";

  // Build context for reference point.
  VldContext refContext(vc.GetDetector(),
                        vc.GetSimFlag(),
                        fReferenceTime);
  
  // Get the reference point data.
  fRefTable.NewQuery(refContext,fTask);

  // The reference value:
  FloatErr reference(1,0);

  // Maybe we don't have any data there....
  if(fRefTable.GetNumRows() == 0) {
    IncrementErrors(kDriftCalibrator,kMissingTable);
    MAXMSG("Calib",Msg::kError,10) 
      << "No rows in CALDRIFT reference table task=" << fTask
      << "  cx=" << refContext.AsString() << std::endl;
    MAXMSG("Calib",Msg::kError,10) 
      << "Your database is bad or MuonDriftCalScheme is misconfigured!" << endl;
    // Flag it.
    fCurError=1;
    refRecoVersion = "Null";
    assert(fRefTable.GetNumRows() && "No Table row for reference date"); 
 } else {
    const CalDrift* refRow = fRefTable.GetRow(0); // We know at least one row exists.
    reference = refRow->GetMedianAsR1_18_2(vc.GetDetector());
    reference.SetError(0); // The error on the reference point is not relevant.
    refRecoVersion = refRow->GetRecoVersion();
 }
  
  // The current value:
  FloatErr current(1,0);
  
  fCurTable.NewQuery(vc,fTask);

  if(fCurTable.GetNumRows() == 0) {
    IncrementErrors(kDriftCalibrator,kMissingTable);

    MAXMSG("Calib",Msg::kError,10) 
      << "No rows in CALDRIFT table task=" << fTask
      << "  cx=" << vc.AsString() 
//       << "  cx=" << refContext.AsString() 
      << ". Trying extrapolation." << std::endl;
    // Ok, now we have work to do.

    VldTimeStamp prevTime;
    const CalDrift* prevRow =0;
    VldTimeStamp nextTime;
    const CalDrift* nextRow =0;

    // Attempt to find the previous valid data.
    const char* sqlprev = Form("(TIMEEND<='%s') "
                              "and (TASK=%d) "
                              "and (DETECTORMASK & %d) "
                              "and (SIMMASK & %d) "
                              "order by TIMEEND desc limit 1", 
                               vc.GetTimeStamp().AsString("s"), 
                               fTask,                    
                               vc.GetDetector(),   
                               vc.GetSimFlag() );
    //MSG("Calib",Msg::kDebug) << "PrevRequest: " << sqlprev << endl;        
    DbiSqlContext prevContext(sqlprev);
    DbiResultPtr<CalDrift> prevTable("CALDRIFT",prevContext);
    if(prevTable.GetNumRows()>0) {
      prevRow = prevTable.GetRow(0);
      VldRange r = prevTable.GetValidityRec(prevRow)->GetVldRange();
      prevTime = r.GetTimeStart() + 0.5*(r.GetTimeEnd() - r.GetTimeStart());
    }    
    
    // Attempt to find the next valid data.
    const char* sqlnext = Form("(TIMESTART>='%s') "
                              "and (TASK=%d) "
                              "and (DETECTORMASK & %d) "
                              "and (SIMMASK & %d) "
                              "order by TIMESTART asc limit 1", 
                               vc.GetTimeStamp().AsString("s"), 
                               fTask,                    
                               vc.GetDetector(),   
                               vc.GetSimFlag() );
    //MSG("Calib",Msg::kDebug) << "NextRequest: " << sqlnext << endl;        
    DbiSqlContext nextContext(sqlnext);
    DbiResultPtr<CalDrift> nextTable("CALDRIFT",nextContext);
    if(nextTable.GetNumRows()>0) {
      nextRow = nextTable.GetRow(0);
      VldRange r = nextTable.GetValidityRec(nextRow)->GetVldRange();
      nextTime = r.GetTimeStart() + 0.5*(r.GetTimeEnd() - r.GetTimeStart());
   }    
  
    if((nextRow)&&(prevRow)) {

      // Interpolate between these two straddling points.
      //MSG("Calib",Msg::kDebug) << "Interpolating.      ";
      double dt = (nextTime-prevTime);
      Float_t y1a=prevRow->GetMedianAsR1_18_2(vc.GetDetector());
      Float_t y1b=prevRow->GetMedianErrAsR1_18_2(vc.GetDetector());

      FloatErr y1(y1a, y1b);
      Float_t y2a=nextRow->GetMedianAsR1_18_2(vc.GetDetector());
      Float_t y2b=nextRow->GetMedianErrAsR1_18_2(vc.GetDetector());
 
      FloatErr y2(y2a,y2b);
      FloatErr y = (y2-y1)*(vc.GetTimeStamp()-prevTime)/dt + y1;
      if(fCurError==0)
        fCurDrift = y/reference;
      
    } else if (nextRow) {

      // Extrapolate from the next point.
      //MSG("Calib",Msg::kDebug) << "Extrapolating next. ";
      Float_t ya=nextRow->GetMedianAsR1_18_2(vc.GetDetector());
      Float_t yb= nextRow->GetMedianErrAsR1_18_2(vc.GetDetector());
      FloatErr y(ya,yb);
      
      // Degrade the error by 2% per month it's out of date.
      y *= FloatErr(1.0, 0.02 * fabs(nextTime-vc.GetTimeStamp())/
                    (30*Munits::day));
      if(fCurError==0)
        fCurDrift = y/reference;
      
    } else if (prevRow) {
     
      // Extrapolate from the previous point.
      MSG("Calib",Msg::kDebug) << "Extrapolating prev. ";
      Float_t ya=prevRow->GetMedianAsR1_18_2(vc.GetDetector());
      Float_t yb=prevRow->GetMedianErrAsR1_18_2(vc.GetDetector());
     
      FloatErr y(ya,yb);
      // Degrade the error by 2% per day it's out of date.
      y *= FloatErr(1.0, 0.02 * fabs(prevTime-vc.GetTimeStamp())/
                    (30*Munits::day));
      if(fCurError==0)
        fCurDrift = y/reference;
      
    } else {
      // No data at all!
      fCurError = 1;
      MAXMSG("Calib",Msg::kError,10) 
        << "No rows in CALDRIFT task=" << fTask
        << "for ANY time!  Badness." << std::endl;
    }
  } else {
    // We have current data, so all the above is unneccessary.
    MSG("Calib",Msg::kDebug) << "Have data.          ";
    const CalDrift* nowRow = fCurTable.GetRow(0);
    Float_t ya=nowRow->GetMedianAsR1_18_2(vc.GetDetector());
    Float_t yb=nowRow->GetMedianErrAsR1_18_2(vc.GetDetector());
    current.Set(ya,yb); 
        
    if(fCurError == 0)
      fCurDrift = current/reference;
  }
  MSG("Calib",Msg::kDebug) << "Drift constant = " << fCurDrift.AsString()
                           << " for time " << vc.AsString() << std::endl;
  
  
}


//......................................................................
void MuonDriftCalScheme::PrintConfig( std::ostream& os ) const
{
  os << " Muon Drift Scheme: " << endl;
  os << "  ReferenceTime:     " << fReferenceTime.AsString() << endl;
  os << "  Task:              " << fTask << endl;
}


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

  if(fCurError) {
    IncrementErrors(kDriftCalibrator,kMissingRow);
  }
  return rawcharge/fCurDrift;
}

//......................................................................
FloatErr MuonDriftCalScheme::DecalDrift(FloatErr undrifted, 
                                        const PlexStripEndId& /*seid*/) const
{

  if(fCurError) {
    IncrementErrors(kDriftCalibrator,kMissingRow);
  }
  return undrifted*fCurDrift;
}

---