StripAttenCalScheme.cxx

Go to the documentation of this file.

// $Id: StripAttenCalScheme.cxx,v 1.20 2007/03/08 21:19:18 tagg Exp $
//
// StripAtten Calibrator class
//
// Responsible for correcting light output from center of strip
// to position along strip.
//
// Nathaniel Tagg
// n.tagg1@physics.ox.ac.uk
#include "StripAttenCalScheme.h"
#include "CalibrationSmearer.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "MessageService/MsgService.h"

#include <cmath>

ClassImp(StripAttenCalScheme)

CVSID( " $Id: StripAttenCalScheme.cxx,v 1.20 2007/03/08 21:19:18 tagg Exp $ ");

StripAttenCalScheme::StripAttenCalScheme()
  : fResPtr(),
    fDefaultAtten(0)
{
  Registry r;
  r.Set("defF12",0.333);
  r.Set("defL1",1.0); 
  r.Set("defL2",7.0);
  r.Set("Task",0); 
  InitializeConfig(r);
}


//......................................................................
StripAttenCalScheme::~StripAttenCalScheme()
{
  if(fDefaultAtten) delete fDefaultAtten;
}


//......................................................................
void StripAttenCalScheme::DoReset( const VldContext& vc )
{
  MSG("Calib",Msg::kVerbose) << "StripAttenCalScheme::DoReset(VldContext) Task: " << fTask << endl;

  fResPtr.NewQuery(vc,fTask);
  if(fResPtr.GetNumRows()==0) {
    MAXMSG("Calib",Msg::kWarning,10) 
      << "StripAtten Scheme: No rows found for CALSTRIPATTEN database table"
      << vc.AsString() << ".\n";         
    IncrementErrors(kAttenCalibrator,kMissingTable);
  }
}

//......................................................................
void StripAttenCalScheme::ConfigModified()
{

  double f12,l1,l2;
  bool ok = true;
  ok = ok && GetConfig().Get("Task",fTask);
  ok = ok && GetConfig().Get("defF12",f12);
  ok = ok && GetConfig().Get("defL1",l1);
  ok = ok && GetConfig().Get("defL2",l2);
  if(!ok) MSG("Calib",Msg::kError) << "Problem setting up StripAttenCalScheme config." << endl;

  if(fDefaultAtten) delete fDefaultAtten;
  
  // Create a dummy seid. This is because the row constructor calls
  // seid.BuildPlnStripEndKey(), which prints out warning messages we don't want.
  // Doesn't matter, we don't use this seid for anything.
  PlexStripEndId dummy_seid(Detector::kNear,1,1,StripEnd::kPositive);
    
  // Build the 'default' row.  Uses the default values harvested above to create a
  // single DB row, which can do our calibration!
  fDefaultAtten = new CalStripAtten(dummy_seid,
                                    l1, l2, f12,  // L1, L2, N1/(N1+N2)
                                    l1/2., l2/2., f12/2.  // 50% errors
                                    );
  
  // Ensure that the DB has been changed for this event.
  Reset(GetContext(),true); 
}

//......................................................................
void StripAttenCalScheme::PrintConfig( std::ostream& os ) const
{
  os << " (NEW attenuation calibration scheme) " << endl;
  os << "  Task = " << fTask << endl;
  os << "  Default attenuation is " 
     << "F12 = " << fDefaultAtten->GetFrac1() << "\t"
     << "Lambda1 = " << fDefaultAtten->GetLambda1() << "\t"
     << "Lambda2 = " << fDefaultAtten->GetLambda2() << endl;
}


//......................................................................
FloatErr StripAttenCalScheme::GetAttenCorrected(FloatErr sigcorr, 
                                                FloatErr stripX, 
                                                const PlexStripEndId& seid) const
{
  // This bit ripped from the old StripCalibrator:
  // Change from 'stripX' to 'length of green fibre'

  UgliGeomHandle ugli(GetContext());
  UgliStripHandle ustrip = ugli.GetStripHandle(seid);

  float outside = fabs(stripX) - ustrip.GetHalfLength();
  if(outside>0) {
    if(outside > 10.0*Munits::cm) {
      MSG("Calib",Msg::kDebug) 
        << "Reconstruction asked for non-physical attenuation correction: "
        << seid.AsString()
        << endl
        << "halflen=" << ustrip.GetHalfLength()
        << "  requested=" << stripX 
        << "  ... trimming to physical bounds now..."
        << endl;
      IncrementErrors(kAttenCalibrator,kBadInput);
    }
    // Trim it.
    stripX = TMath::Sign(ustrip.GetHalfLength(),stripX);
  }
  
  double halfgreen = ustrip.GetHalfLength() + ustrip.WlsPigtail(seid.GetEnd());
  double wlsLen = halfgreen;
  if(seid.GetEnd()==StripEnd::kPositive) wlsLen -= stripX;
  if(seid.GetEnd()==StripEnd::kNegative) wlsLen += stripX;
  // If not +ve or -ve, do nothing.

  const CalStripAtten* atten = fResPtr.GetRowByIndex(seid.BuildPlnStripEndKey());
  if(atten==0) {
    if((fResPtr.GetNumRows()>0)&&(!seid.IsVetoShield())) {
      MAXMSG("Calib",Msg::kWarning,10) << "StripAtten Scheme: Cannot find row for stripend "
                                       << seid.AsString() << ". Using default." << endl;
      IncrementErrors(kAttenCalibrator,kMissingRow,seid);
    }
    atten = fDefaultAtten;
  }

  // The correction is for the center of the strip to the new position
  // on the strip.
  
  float    a = atten->GetAttenuation(halfgreen);
  FloatErr b = atten->GetAttenuation(wlsLen);
  b.SetError( fabs( b.GetValue() - atten->GetAttenuation(wlsLen - stripX.GetError()) ) );

  if ( b < 4.0e-10 ) {  
    // The value of "b" gets raised to the 4-th power in
    // the error calculation of a/b. We cut below 4.0e-10 because it
    // is ~(1.175e-38)^(1/4) where that is the smallest representable
    // value for a Float_t; smaller values underflow to zero and are
    // then used as a denominator.
    MAXMSG("Calib",Msg::kWarning,10) 
      << "StripAtten Scheme: caught a divide-by-zero error on " 
      << seid.AsString() << endl;
    IncrementErrors(kAttenCalibrator,kFPE,seid);
    return sigcorr * FloatErr(1,0.5);
  }

  // New, March 2007. NJT, JL
  // Clamp the correction at the near end to be less that a factor of 10
  FloatErr ratio = a/b;
  if(ratio < 0.1) ratio = FloatErr(0.1,fabs(0.1-ratio.GetValue()));
  
  return sigcorr * ratio;
}
  

//......................................................................
FloatErr StripAttenCalScheme::DecalAttenCorrected(FloatErr sigmap,
                                                FloatErr stripX,
                                                const PlexStripEndId& seid) const
{

  UgliGeomHandle ugli(GetContext());
  UgliStripHandle ustrip = ugli.GetStripHandle(seid);
  
  // Don't clip the MC.. trust that it's checked strip bounds.
  // I sometimes want to inject light into the pigtails, so I need the
  // flexibility. --Nathaniel, Jan 2006
//   if(fabs(stripX) > ustrip.GetHalfLength()) {
//     MSG("Calib",Msg::kDebug) << "Reconstruction asked for non-physical attenuation correction: "
//                           << seid.AsString() << endl;
//     MSG("Calib",Msg::kDebug) << "halflen=" << ustrip.GetHalfLength()
//                           << "  requested=" << stripX 
//                           << "  ... trimming to physical bounds now..."
//                           << endl;
    
//     // Increment error if it's more than 10 cm out.
//     if(fabs(stripX) > ustrip.GetHalfLength()+0.1) 
//       IncrementErrors(kAttenCalibrator,kBadInput);

//     // Trim it.
//     stripX = TMath::Sign(ustrip.GetHalfLength(),stripX);
//   }

  double halfgreen = ustrip.GetHalfLength() + ustrip.WlsPigtail(seid.GetEnd());
  double wlsLen = halfgreen;
  if(seid.GetEnd()==StripEnd::kPositive) wlsLen -= stripX;
  if(seid.GetEnd()==StripEnd::kNegative) wlsLen += stripX;
  // If not +ve or -ve, do nothing.

  const CalStripAtten* atten = fResPtr.GetRowByIndex(seid.BuildPlnStripEndKey());
  if(atten==0) {
    if((fResPtr.GetNumRows()>0)&&(!seid.IsVetoShield())) {
      MAXMSG("Calib",Msg::kWarning,10) << "StripAtten Scheme: Cannot find row for stripend "
                                       << seid.AsString() << ". Using default." << endl;
      IncrementErrors(kAttenCalibrator,kMissingRow,seid);
    }
    atten = fDefaultAtten;
  }

  // The correction is for the center of the strip to the new position
  // on the strip.
  
  float a = atten->GetAttenuation(halfgreen);
  FloatErr b = atten->GetAttenuation(wlsLen);
  b.SetError( fabs( b.GetValue() - atten->GetAttenuation(wlsLen - stripX.GetError()) ) );

  if(a < 4e-10) {
    MAXMSG("Calib",Msg::kWarning,10) 
      << "StripAtten Scheme: caught a divide-by-zero error on " << seid.AsString() << endl;
    IncrementErrors(kAttenCalibrator,kFPE,seid);
    return sigmap * FloatErr(1,0.5);
  }
  
  // New, March 2007. NJT, JL
  // Clamp the correction at the near end to be less that a factor of 10
  FloatErr ratio = b/a;
  if(ratio > 10.) ratio = FloatErr(10.,fabs(10.-ratio.GetValue()));
  
  return sigmap * ratio;
}

---