MuonDriftCalScheme.cxx

Go to the documentation of this file.

// $Id: MuonDriftCalScheme.cxx,v 1.11 2007/04/01 14:00:20 evans 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.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.11 2007/04/01 14:00:20 evans Exp $");


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

   Registry r;
   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 = "";

  //Make all drift numbers consistent with R1_18_2.
  Double_t r1_18Factor = 1.0; //R1_18_2 = r1_18Factor*R1_18;
  Double_t r1_18_2Factor = 1.0;
  Double_t r1_18_4Factor = 1.0;
  Double_t cedarFactor = 1.0;
  if (Detector::kNear == vc.GetDetector()){
    r1_18Factor = 1.00129; //R1_18_2 = r1_18Factor*R1_18;
    r1_18_2Factor = 1.0;
    r1_18_4Factor = 1.0;
    cedarFactor = 1.0007; //R1_18_2 = cedarFactor*cedar;
  }
  if (Detector::kFar == vc.GetDetector()){
    r1_18Factor = 1.0;
    r1_18_2Factor = 1.0;
    r1_18_4Factor = 1.0;
    cedarFactor = 0.96688; //R1_18_2 = cedarFactor*cedar;
  }

  // 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";
  } else {
    const CalDrift* refRow = fRefTable.GetRow(0); // We know at least one row exists.
    reference = refRow->GetMedian();
    //    if ("R1_18" == refRow->GetRecoVersion()){reference *= corrFact;}
    if ("R1_18" == refRow->GetRecoVersion()){reference *= r1_18Factor;}
    if ("R1_18_2" == refRow->GetRecoVersion()){reference *= r1_18_2Factor;}
    if ("R1_18_4" == refRow->GetRecoVersion()){reference *= r1_18_4Factor;}
    if ("Cedar" == refRow->GetRecoVersion()){reference *= cedarFactor;}
    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 prevFact = 1.0;
      Float_t nextFact = 1.0;
      if ("R1_18" == nextRow->GetRecoVersion()){nextFact = r1_18Factor;}
      if ("R1_18_2" == nextRow->GetRecoVersion()){nextFact = r1_18_2Factor;}
      if ("R1_18_4" == nextRow->GetRecoVersion()){nextFact = r1_18_4Factor;}
      if ("Cedar" == nextRow->GetRecoVersion()){nextFact = cedarFactor;}
      if ("R1_18" == prevRow->GetRecoVersion()){prevFact = r1_18Factor;}
      if ("R1_18_2" == prevRow->GetRecoVersion()){prevFact = r1_18_2Factor;}
      if ("R1_18_4" == prevRow->GetRecoVersion()){prevFact = r1_18_4Factor;}
      if ("Cedar" == prevRow->GetRecoVersion()){prevFact = cedarFactor;}
      FloatErr y1(prevRow->GetMedian(), prevRow->GetMedianErr());
      FloatErr y2(nextRow->GetMedian(), nextRow->GetMedianErr());
      FloatErr y = (y2*nextFact-y1*nextFact)*(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. ";
      FloatErr y(nextRow->GetMedian(), nextRow->GetMedianErr());
      if ("R1_18" == nextRow->GetRecoVersion()){y *= r1_18Factor;}
      if ("R1_18_2" == nextRow->GetRecoVersion()){y *= r1_18_2Factor;}
      if ("R1_18_4" == nextRow->GetRecoVersion()){y *= r1_18_4Factor;}
      if ("Cedar" == nextRow->GetRecoVersion()){y *= cedarFactor;}
      // 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. ";
      FloatErr y(prevRow->GetMedian(), prevRow->GetMedianErr());
      if ("R1_18" == prevRow->GetRecoVersion()){y *= r1_18Factor;}
      if ("R1_18_2" == prevRow->GetRecoVersion()){y *= r1_18_2Factor;}
      if ("R1_18_4" == prevRow->GetRecoVersion()){y *= r1_18_4Factor;}
      if ("Cedar" == prevRow->GetRecoVersion()){y *= cedarFactor;}
      // 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);
    current.Set(nowRow->GetMedian(), nowRow->GetMedianErr());
    if ("R1_18" == nowRow->GetRecoVersion()){current *= r1_18Factor;}
    if ("R1_18_2" == nowRow->GetRecoVersion()){current *= r1_18_2Factor;}
    if ("R1_18_4" == nowRow->GetRecoVersion()){current *= r1_18_4Factor;}
    if ("Cedar" == nowRow->GetRecoVersion()){current *= cedarFactor;}
    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;
}

---