AlgFarDetEvent.cxx

Go to the documentation of this file.

#include "AlgFarDetEvent.h"
#include "FarDetEventHandle.h"

#include "Algorithm/AlgConfig.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Candidate/CandContext.h"

#include "RecoBase/CandFitTrackHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandStripHandle.h"

#include "Calibrator/Calibrator.h"

#include "MuELoss/BetheBlochModel.h"
#include "MuELoss/StandardMaterial.h"
#include "MuELoss/MuEnergyLoss.h"

#include "TString.h"
#include <cmath>

ClassImp(AlgFarDetEvent)    

CVSID("$Id: AlgFarDetEvent.cxx,v 1.11 2009/02/13 15:03:01 blake Exp $");

AlgFarDetEvent::AlgFarDetEvent()
{

}   

AlgFarDetEvent::~AlgFarDetEvent()
{

}

// Event Building Algorithm
// ========================
//
// (0) Set primary track/shower properties; basic event properties.
// (1) Identify shower strips and calculate MIPs for each strip.
// (2) Make correction to MIPs to account for track energy deposits.
// (3) Add up the energy deposits linearly and calibrate to Ehad,
//     use this value to decide which deweighting should be used.
// (4) Add up the deweighted energy and calibrate to Ehad.
// (5) Set different energy values
//

void AlgFarDetEvent::RunAlg(AlgConfig& ac, CandHandle& ch, CandContext& cx)
{
  MSG("FarDetEvent",Msg::kDebug) << " AlgFarDetEvent::RunAlg(....) " << endl;

  FarDetEventHandle& myevent = dynamic_cast<FarDetEventHandle&>(ch);

  Calibrator& cal = Calibrator::Instance();
  cal.ReInitialise(*(cx.GetCandRecord()->GetVldContext()));

  // Unpack CandContext
  // ==================
  const TObjArray* eventArray = dynamic_cast<const TObjArray*>(cx.GetDataIn());

  CandFitTrackHandle* fit = (CandFitTrackHandle*)(eventArray->At(0));
  CandTrackHandle* trk = (CandTrackHandle*)(eventArray->At(1));
  CandShowerHandle* shw = (CandShowerHandle*)(eventArray->At(2));
  // TObjArray* HadronicShowers = (TObjArray*)(eventArray->At(3));
  // TObjArray* BremsAndDeltas = (TObjArray*)(eventArray->At(4));
  // TObjArray* AssociatedStrips = (TObjArray*)(eventArray->At(5));
  TObjArray* AllStrips = (TObjArray*)(eventArray->At(6));
  
  // Set primary track and shower properties
  // =======================================
  MSG("FarDetEvent",Msg::kDebug) << " Setting track properties " << endl;
  this->SetTrackProperties(&myevent, fit, trk);

  MSG("FarDetEvent",Msg::kDebug) << " Setting shower properties " << endl;
  this->SetShowerProperties(&myevent, shw);

  if( myevent.GetPrimaryTrack()==0 && myevent.GetPrimaryShower()==0 )
    {
      MSG("FarDetEvent",Msg::kWarning) << " AlgFarDetEvent: Event has no primary track or shower!" << endl; 
    }

  // Calculate other event properties
  // ================================
  MSG("FarDetEvent",Msg::kDebug) << " Setting event properties " << endl;

  int minplane(-999), maxplane(-999);
  double minz(-999.9), maxz(-999.9);

  for(Int_t i=0; i<AllStrips->GetEntries(); i++)
    {
      CandStripHandle* strip = (CandStripHandle*)(AllStrips->At(i));
    
      if( minplane<0 || strip->GetPlane()<minplane ){
        minplane = strip->GetPlane(); 
        minz = strip->GetZPos();
      }
 
      if( maxplane<0 || strip->GetPlane()>maxplane ){
        maxplane = strip->GetPlane();  
        maxz = strip->GetZPos();
      }

    }

  myevent.SetMinPlane(minplane); 
  myevent.SetMaxPlane(maxplane);
  myevent.SetMinZ(minz); 
  myevent.SetMaxZ(maxz);

  // Calculate event kinematics
  // ==========================
  MSG("FarDetEvent",Msg::kDebug) << " Calculating event kinematics " << endl;

  double EhadLinear(0.), EhadDeweight(0.);

  // (1) get MIP values for each shower strip
  // ========================================
  // input: primary shower, primary track
  // output: list of MIPs, track correction, shared planes 

  MSG("FarDetEvent",Msg::kDebug) << " (1) calculate MIP values for each strip " << endl;

  vector<double> TempShw;
  vector<double> TempTrk;

  this->GetShower( myevent.GetPrimaryShower(), myevent.GetPrimaryTrack(), TempShw, TempTrk );

  MSG("FarDetEvent",Msg::kDebug) << "    ... shower strips = " << TempShw.size() << endl;

  // (2) calculate pulse height and energy
  // =====================================
  // input: list of MIPs, track correction, shared planes, deweight factor
  // output: corrected pulse height, then shower energy

  MSG("FarDetEvent",Msg::kDebug) << " (2) calculate pulse height and energy " << endl;

  // linear pulse height and energy
  double linearfactor(1.0), linearph(0.0);

  this->GetPulseHeight( TempShw, TempTrk, linearph, linearfactor );

  EhadLinear = this->CalibrateShowerEnergy( linearph, linearfactor, ac );

  MSG("FarDetEvent",Msg::kDebug) << "    ... linear pulse height = " << linearph << endl;
  MSG("FarDetEvent",Msg::kDebug) << "        linear energy = " << EhadLinear << endl;

  // deweighted pulse height and energy
  double deweightfactor(1.0), deweightph(0.0);

  deweightfactor = this->GetDeWeightFactor( EhadLinear, ac );

  this->GetPulseHeight( TempShw, TempTrk, deweightph, deweightfactor );

  EhadDeweight = this->CalibrateShowerEnergy( deweightph, deweightfactor, ac );

  MSG("FarDetEvent",Msg::kDebug) << "    ... deweighted pulse height = " << deweightph << endl;
  MSG("FarDetEvent",Msg::kDebug) << "        deweighted energy = " << EhadDeweight << endl;

  // set reconstructed shower energy
  myevent.SetShwMipsLinear(linearph);            
  myevent.SetShwMipsDeweighted(deweightph);        
  myevent.SetShwEnergyGeVLinear(EhadLinear); 
  myevent.SetShwEnergyGeVDeweighted(EhadDeweight); 
  myevent.SetShwEnergyGeV(EhadLinear); 

  // (3) calculate overall energy
  // ============================
  double Enu(0.0);
  
  myevent.SetEnergy(Enu);

  return;
}

// ==========================================================
//  SET PRIMARY TRACK PROPERTIES
// ==========================================================
void AlgFarDetEvent::SetTrackProperties(FarDetEventHandle* evt, CandFitTrackHandle* fit, CandTrackHandle* trk)
{
  MSG("FarDetEvent",Msg::kVerbose) << " AlgFarDetEvent::SetTrackProperties(....) " << endl;

  CandTrackHandle* track = 0; 

  double erange,ecurve,emcharge;
  // Use the fitted track if we have one...
  if(fit) 
    {
      track = fit;
      erange = fit->GetMomentumRange();
      ecurve = fabs(fit->GetMomentumCurve());
      emcharge = fit->GetEMCharge();
    }  
  // ...or else fall back to found track
  else if( trk ) 
    {
      track = trk;
      erange = trk->GetMomentum();
      ecurve = 0.0;
      emcharge = 0.0;
    }
  
  if( track )
    {
      // set primary track properties
      MSG("FarDetEvent",Msg::kVerbose) << "   (found primary track) " << endl;
      evt->AddTrack(track);
      evt->SetPrimaryTrack(track);
      evt->SetMuReco(1);
      evt->SetMuMomentumRange(erange);
      evt->SetMuMomentumCurve(ecurve);
      evt->SetEMCharge(emcharge);
      evt->SetMuDirCosU(track->GetDirCosU()); 
      evt->SetMuDirCosV(track->GetDirCosV()); 
      evt->SetMuDirCosZ(track->GetDirCosZ());
      evt->SetVtxTime(track->GetVtxT());
      evt->SetVtxU(track->GetVtxU()); 
      evt->SetVtxV(track->GetVtxV()); 
      evt->SetVtxZ(track->GetVtxZ());
      evt->SetVtxPlane(track->GetVtxPlane()); 
    }

  MSG("FarDetEvent",Msg::kVerbose) << endl
    << " Track Properties: " << endl
    << "  MuReco = " << evt->GetMuReco() << endl
    << "  MuMomentumRange = " << evt->GetMuMomentumRange() << endl
    << "  MuMomentumCurve = " << evt->GetMuMomentumCurve() << endl
    << "  EMCharge = " << evt->GetEMCharge() << endl
    << "  MuDirCosU = " << evt->GetMuDirCosU() << endl
    << "  MuDirCosV = " << evt->GetMuDirCosV() << endl
    << "  MuDirCosZ = " << evt->GetMuDirCosZ() << endl   
    << "  (VtxTime = " << evt->GetVtxTime() << ")" << endl
    << "  (VtxU = " << evt->GetVtxU() << ")" << endl
    << "  (VtxV = " << evt->GetVtxV() << ")" << endl
    << "  (VtxZ = " << evt->GetVtxZ() << ")" << endl
    << "  (VtxPlane = " << evt->GetVtxPlane() << ")" << endl;

  return;
}

// ==========================================================
//  SET PRIMARY SHOWER PROPERTIES
// ==========================================================

void AlgFarDetEvent::SetShowerProperties(FarDetEventHandle* evt, CandShowerHandle* shw)
{
  MSG("FarDetEvent",Msg::kVerbose) << " AlgFarDetEvent::SetShowerProperties(....) " << endl;

  if( shw )
    {
      // set primary shower properties
      MSG("FarDetEvent",Msg::kVerbose) << "   (found primary shower) " << endl;
      evt->SetShwReco(1);
      evt->AddShower(shw);
      evt->SetPrimaryShower(shw);
      evt->SetShwDirCosU(shw->GetDirCosU()); 
      evt->SetShwDirCosV(shw->GetDirCosV()); 
      evt->SetShwDirCosZ(shw->GetDirCosZ());
      if( evt->GetMuReco()==0 )
        {
          evt->SetVtxTime(shw->GetVtxT());
          evt->SetVtxU(shw->GetVtxU()); 
          evt->SetVtxV(shw->GetVtxV()); 
          evt->SetVtxZ(shw->GetVtxZ());
          evt->SetVtxPlane(shw->GetVtxPlane()); 
        }
    }

  MSG("FarDetEvent",Msg::kVerbose) << endl
    << " Shower properties: " << endl
    << "  ShwReco = " << evt->GetShwReco() << endl
    << "  ShwDirCosU = " << evt->GetShwDirCosU() << endl
    << "  ShwDirCosV = " << evt->GetShwDirCosV() << endl
    << "  ShwDirCosZ = " << evt->GetShwDirCosZ() << endl   
    << "  (VtxTime = " << evt->GetVtxTime() << ")" << endl
    << "  (VtxU = " << evt->GetVtxU() << ")" << endl
    << "  (VtxV = " << evt->GetVtxV() << ")" << endl
    << "  (VtxZ = " << evt->GetVtxZ() << ")" << endl
    << "  (VtxPlane = " << evt->GetVtxPlane() << ")" << endl;

  return;
}

// ==========================================================
// METHODS FOR SHOWER ENERGY CALCULATION
// ==========================================================

bool AlgFarDetEvent::Contained(const CandTrackHandle* track)
{
  // vertex containment
  double VtxU = track->GetVtxU();
  double VtxV = track->GetVtxV();
  double VtxR = sqrt(VtxU*VtxU+VtxV*VtxV);
  bool VtxContained = (VtxR<3.75 && VtxR>0.4) && ((track->GetVtxPlane()>5 && track->GetVtxPlane()<243) || (track->GetVtxPlane()>254 && track->GetVtxPlane()<480));

  // end containment
  double EndU = track->GetEndU();
  double EndV = track->GetEndV();
  double EndR = sqrt(EndU*EndU+EndV*EndV);
  bool EndContained = (EndR<3.75 && EndR>0.4) && ((track->GetEndPlane()>5 && track->GetEndPlane()<243) || (track->GetEndPlane()>254 && track->GetEndPlane()<480));

  return (VtxContained && EndContained);
}

double AlgFarDetEvent::DeDx(Float_t emu)
{
  double dedx = 0.;
  double mm = 105.658389;
  double mm_2 = mm*mm;
  
  if(emu*emu-mm_2<0.){return 0.;}
  double a_2     = 1./(137.036*137.036);
  double pi = 3.141;
  double Na = 6.023;
  double lamda_2 = 3.8616*3.8616;
  double Z_A = 0.5377;
  double me  = 0.51099906;
  double me_2 = me*me; 
  double beta = sqrt(emu*emu-mm_2)/emu;
  double beta_2 = beta*beta;
  double gamma = emu/105.658389;
  double gamma_2 = gamma*gamma;
  double I = 68.7e-6;
  double I_2 = I*I;
  double p = sqrt(emu*emu-mm_2);
  double p_2 = p*p;
  double Em = 2 * me * p_2 / ( me_2 + mm_2 + 2*me*emu );
  double Em_2= Em*Em;
  double emu_2 = emu*emu;
  double X0 = 0.165;
  double X1 = 2.503;
  double X = log10(beta*gamma);
  double a = 0.165;
  double C = -3.3;
  double m = 3.222;
  double d = 0;
  if(X0 < X && X < X1){d = 4.6052 * X + a * pow(X1-X,m) + C;}
  if(X > X1){d = 4.6052 * X + C;}

  dedx =  a_2 * 2*pi * Na * lamda_2 * Z_A * (me / beta_2) *( log( 2*me*beta_2*gamma_2*Em/I_2 ) - 2*beta_2 + 0.25*(Em_2/emu_2) - d );

  return 10*dedx;
}

 // GetShower 
 // =========
 // This stores the MIP values of strips contained in the shower in 
 // TempShw and calculates the number of planes where the shower and 
 // track overlap as well as the amount of pulseheight to remove per 
 // strip to account for the energy deposited by the muon

void AlgFarDetEvent::GetShower(const CandShowerHandle* shower, 
                               const CandTrackHandle* track, 
                               vector<double>& TempShw,
                               vector<double>& TempTrk)
{

  Calibrator& cal = Calibrator::Instance();
  int indx(-1),pln(-1),strp(-1);


  // make map of shower strips
  // ========================
  map<int,int> isShowerStrip;
  if( shower ) 
    {
      TIter shwitr(shower->GetDaughterIterator());
      while(CandStripHandle* strip = (CandStripHandle*)(shwitr()))
        { 
          if( strip )
            {
              pln = strip->GetPlane();
              strp = strip->GetStrip();
              indx = strp + 192*pln;
              isShowerStrip[indx] = 1;
            }
        }
    }


  // make map of track strips
  // ========================
  map<int,int> isTrackStrip;
  
  double sharedplanes(0.0);
  if( track ) 
    {
      int minpln(-999),maxpln(-999);
      TIter trkitr(track->GetDaughterIterator());
      while(CandStripHandle* strip = (CandStripHandle*)(trkitr()))
        { 
          if( strip )
            {
              pln = strip->GetPlane();
              strp = strip->GetStrip();
              indx = strp + 192*pln;
              isTrackStrip[indx] = 1;

              if( isShowerStrip[indx]==1 )
                {
                  if(minpln<0 || pln<minpln) minpln=pln;
                  if(maxpln<0 || pln>maxpln) maxpln=pln;
                }
            }
        }
    
      if( maxpln>=0 && minpln>=0 && maxpln-minpln>=0 )
        {
          sharedplanes = 1+maxpln-minpln;
        }
    }

  MSG("FarDetEvent",Msg::kVerbose) << "    ... shared planes = " << sharedplanes << endl;    


  // Calculate track correction (MIPs per strip) to remove from shower
  // =================================================================
  double trackcorrection = 0.0;
  if(track)
    { 
       double pmu = track->GetMomentum();
       double trackdircosZ = fabs(1.0/(track->GetVtxDirCosZ())); 
       MSG("FarDetEvent",Msg::kVerbose) << "    ... track direction = " << trackdircosZ << endl;

       const CandFitTrackHandle* fit = dynamic_cast<const CandFitTrackHandle*>(track);
       if( fit )
         {
           if( this->Contained(track) ) pmu = fit->GetMomentumRange(); 
           else pmu = fit->GetMomentumCurve();
         }

       // (Note: alternative dEdX calculation)
       // Material  poly(MinosMaterial::ePolystyreneMinos);
       // BetheBlochModel bbmodel(poly);
       // double dedx = bbmodel.dE_dx(emu*1000.)/1.985; 

       double dedx = DeDx(1000.0*pmu)/1.985;
       double dedx_end = dedx;
       double pmuend = pmu - sharedplanes/(30.*trackdircosZ);
       MSG("FarDetEvent",Msg::kVerbose) << "    ... muon momentum: vtx = " << pmu << " end = " << pmuend << endl;
       if( pmuend>0.3 ) // beta*gamma~3 => p~300MeV
         {
           dedx_end = DeDx(1000.0*pmuend)/1.985;
         }
       
       double avg_dedx = 0.5*(dedx+dedx_end);
       MSG("FarDetEvent",Msg::kVerbose) << "    ... average dedx per strip (MIPs) = " << avg_dedx << endl;

       double temp_trackdircosZ = fabs(trackdircosZ);
       if( temp_trackdircosZ<0.5 ) temp_trackdircosZ = 0.5;
       trackcorrection = avg_dedx/temp_trackdircosZ;
    }
  
  MSG("FarDetEvent",Msg::kVerbose) << "    ... track correction = " << trackcorrection << endl;

  // Calculate the MIP values for each strip in the shower
  // =====================================================
  
  if(shower)
    {
      double opos(0.0);
      double sigmapN(0.0),sigmapP(0.0);
      double shwMIP(0.0),trkMIP(0.0); 
      TIter shwitr(shower->GetDaughterIterator());
      while(CandStripHandle* strip = (CandStripHandle*)(shwitr()))
        { 
          if( strip )
            {
              pln = strip->GetPlane();
              strp = strip->GetStrip();
              indx = strp + 192*pln;

              shwMIP = 0.0; 
              trkMIP = 0.0;
              opos = 0.0;

              if(strip->GetPlaneView()==PlaneView::kU) opos = -1.0*shower->GetV(pln); 
              if(strip->GetPlaneView()==PlaneView::kV) opos =  1.0*shower->GetU(pln); 
              if( fabs(opos)<999.9 )
                {
                  sigmapN=cal.GetAttenCorrected(strip->GetCharge(StripEnd::kNegative,CalDigitType::kSigCorr), opos, strip->GetStripEndId(StripEnd::kNegative) );
                  shwMIP+=cal.GetMIP(sigmapN, strip->GetStripEndId(StripEnd::kNegative) );
                  sigmapP=cal.GetAttenCorrected(strip->GetCharge(StripEnd::kPositive,CalDigitType::kSigCorr), opos, strip->GetStripEndId(StripEnd::kPositive) );
                  shwMIP+=cal.GetMIP(sigmapP, strip->GetStripEndId(StripEnd::kPositive) );
                }
              
              if( isTrackStrip[indx]==1 )
                {
                  trkMIP = trackcorrection;
                }
              TempShw.push_back(shwMIP);
              TempTrk.push_back(trkMIP);
            }   
        }  
    }
  
  return;
}

 // GetDeweightFactor(...)
 // ======================
 // Calculate deweight factor as function of linear energy

double AlgFarDetEvent::GetDeWeightFactor(const double& LinearEhad, AlgConfig & ac)
{
  double deweightfactorCC(1.0);

  double fdeweightLowScale(0.0);
  double fdeweightC0(0.0);
  double fdeweightC1(0.0);
  double fdeweightC2(0.0);
  double fdeweightC3(0.0);

  bool warn =false;
  TString errstring="AlgFarDetEvent::GetDeWeightFactor - Failed to get AlgConfig info for the following keys: ";
  
  if( ac.KeyExists("deweightLowScale") ) fdeweightLowScale = ac.GetDouble("deweightLowScale");
  else {errstring+="deweightLowScale, "; warn =true; }
  if( ac.KeyExists("deweightC0") )fdeweightC0 = ac.GetDouble("deweightC0");
  else {errstring+="deweightC0, "; warn =true; }
  if( ac.KeyExists("deweightC1") )fdeweightC1 = ac.GetDouble("deweightC1");
  else {errstring+="deweightC1, "; warn =true; }
  if( ac.KeyExists("deweightC2") )fdeweightC2 = ac.GetDouble("deweightC2");
  else {errstring+="deweightC2, "; warn =true; }
  if( ac.KeyExists("deweightC3") )fdeweightC3 = ac.GetDouble("deweightC3");
  else {errstring+="deweightC3"; warn =true; }
  if(warn) MSG("FarDetEvent",Msg::kWarning) << errstring.Data() <<endl;

  if(LinearEhad < fdeweightLowScale)
    {
      deweightfactorCC = fdeweightC0 + 
                         fdeweightC1*LinearEhad +
                         fdeweightC2*LinearEhad*LinearEhad +
                         fdeweightC3*LinearEhad*LinearEhad*LinearEhad;
    }
  
  MSG("FarDetEvent",Msg::kVerbose) << " LinearEnergy=" << LinearEhad << " DeweightFactor=" << deweightfactorCC << endl;

  return deweightfactorCC;
}

 // GetPulseHeight(...)
 // ===================
 // Calculate the summed MIPs for the current shower

void AlgFarDetEvent::GetPulseHeight(vector<double>& TempShw, vector<double>& TempTrk, double& totph, const double& deweight)
{

  // zero pulse height
  totph=0.0;
  
  double ph;
  double pulseheight_check_1(0);
  double pulseheight_check_2(0);
  double pulseheight_check_3(0);

  unsigned int nstrips = TempShw.size();
  if( nstrips<1 ) return;

  if( deweight>=1.0 )
    {
      for(unsigned int i=0; i<nstrips; ++i)
        {
          ph = TempShw[i]-TempTrk[i];
          pulseheight_check_1 += TempShw[i];
          pulseheight_check_2 += TempTrk[i];
          totph += ph;
        }
      MSG("FarDetEvent",Msg::kVerbose) << "    ... linear check: ShwPH=" << pulseheight_check_1 << ", TrkPH=" << pulseheight_check_2 << endl;
    }
  else
    {
      for(unsigned int i=0; i<nstrips; ++i)
        {
          if( TempShw[i]-1.4*TempTrk[i]>0.0 )
            {
              ph = TempShw[i]-1.4*TempTrk[i];
              pulseheight_check_3 += ph;
              totph += pow(ph,deweight);
            }
        }
        MSG("FarDetEvent",Msg::kVerbose) << "    ... deweight check: ShwPH-TrkPH=" << pulseheight_check_3 << endl;   
    }

  return;
}

 // CalibrateShowerEnergy(...)
 // ==========================
 // Calculate GeV from MIPs (using polynomial paramaterizations)

double AlgFarDetEvent::CalibrateShowerEnergy(const double& totph, const double& deweight, AlgConfig & ac)
{

  // Initialize shower energy
  double cal_e(0.0);

  if( totph<=0 ) 
    {
      MSG("FarDetEvent",Msg::kVerbose) << " 0 MIPs -> 0 GeV" << endl;
      return cal_e;
    }

  // Get Calibration Constants from AlgConfig.
  if( deweight>=1.0 )
    {
      // USE LINEAR CONSTANTS
      // ====================
      double fCCLinLowScale(0.0); 
      double fCCLinLowC0(0.0); 
      double fCCLinLowC1(0.0); 
      double fCCLinLowC2(0.0); 
      double fCCLinHiC0(0.0); 
      double fCCLinHiC1(0.0);
 
      bool warn =false;
      TString errstring="AlgFarDetEvent::CalibrateShowerEnergy - Linear - Failed to get AlgConfig info for the following keys: ";
      
      if( ac.KeyExists("CCLinLowScale") ) fCCLinLowScale = ac.GetDouble("CCLinLowScale"); 
      else {errstring+="CCLinLowScale, "; warn =true; }
      if( ac.KeyExists("CCLinLowC0") )    fCCLinLowC0 = ac.GetDouble("CCLinLowC0");
      else {errstring+="CCLinLowC0, "; warn =true; }
      if( ac.KeyExists("CCLinLowC1") )    fCCLinLowC1 = ac.GetDouble("CCLinLowC1");
      else {errstring+="CCLinLowC1, "; warn =true; }
      if( ac.KeyExists("CCLinLowC2") )    fCCLinLowC2 = ac.GetDouble("CCLinLowC2");
      else {errstring+="CCLinLowC2, "; warn =true; }
      if( ac.KeyExists("CCLinHiC0") )     fCCLinHiC0 = ac.GetDouble("CCLinHiC0");
      else {errstring+="CCLinHiC0, "; warn =true; }
      if( ac.KeyExists("CCLinHiC1") )     fCCLinHiC1 = ac.GetDouble("CCLinHiC1"); 
      else {errstring+="CCLinHiC1"; warn =true; }

      if(warn) MSG("FarDetEvent",Msg::kWarning) << errstring.Data() <<endl;

      if( totph<fCCLinLowScale )
        {
          MSG("FarDetEvent",Msg::kVerbose) << "  using linear constants, low scale... " << endl;
          cal_e = fCCLinLowC0 + fCCLinLowC1*totph + fCCLinLowC2*totph*totph
                   + 0.21*sqrt(totph)*(1.0-totph/25.0)*exp(-totph/25.0);
        }
      else
        {
          MSG("FarDetEvent",Msg::kVerbose) << "  using linear constants, high scale... " << endl;
          cal_e = fCCLinHiC0 + fCCLinHiC1*totph;
        }
    }
  else
    {
      // USE DEWEIGHTED CONSTANTS
      // ========================
      double fCCWtLowScale(0.0);
      double fCCWtMidScale(0.0);
      double fCCWtLowC0(0.0);
      double fCCWtLowC1(0.0);
      double fCCWtLowC2(0.0);
      double fCCWtLowC3(0.0);
      double fCCWtLowC4(0.0);
      double fCCWtMidC0(0.0);
      double fCCWtMidC1(0.0);
      double fCCWtHiC0(0.0);
      double fCCWtHiC1(0.0);

      bool warn =false;
      TString errstring="AlgFarDetEvent::CalibrateShowerEnergy - Deweight - Failed to get AlgConfig info for the following keys: ";
      
      if( ac.KeyExists("CCWtLowScale") ) fCCWtLowScale = ac.GetDouble("CCWtLowScale");  
      else {errstring+="CCWtLowScale, "; warn =true; }
      if( ac.KeyExists("CCWtMidScale") ) fCCWtMidScale = ac.GetDouble("CCWtMidScale");  
      else {errstring+="CCWtMidScale, "; warn =true; } 
      if( ac.KeyExists("CCWtLowC0") )    fCCWtLowC0 = ac.GetDouble("CCWtLowC0");
      else {errstring+="CCWtLowC0, "; warn =true; }
      if( ac.KeyExists("CCWtLowC1") )    fCCWtLowC1 = ac.GetDouble("CCWtLowC1");
      else {errstring+="CCWtLowC1, "; warn =true; }
      if( ac.KeyExists("CCWtLowC2") )    fCCWtLowC2 = ac.GetDouble("CCWtLowC2");
      else {errstring+="CCWtLowC2, "; warn =true; }
      if( ac.KeyExists("CCWtLowC3") )    fCCWtLowC3 = ac.GetDouble("CCWtLowC3");
      else {errstring+="CCWtLowC3, "; warn =true; }
      if( ac.KeyExists("CCWtLowC4") )    fCCWtLowC4 = ac.GetDouble("CCWtLowC4");
      else {errstring+="CCWtLowC4, "; warn =true; }
      if( ac.KeyExists("CCWtMidC0") )    fCCWtMidC0 = ac.GetDouble("CCWtMidC0");
      else {errstring+="CCWtMidC0, "; warn =true; }
      if( ac.KeyExists("CCWtMidC1") )    fCCWtMidC1 = ac.GetDouble("CCWtMidC1"); 
      else {errstring+="CCWtMidC1"; warn =true; }
      if( ac.KeyExists("CCWtHiC0") )     fCCWtHiC0 = ac.GetDouble("CCWtHiC0");
      else {errstring+="CCWtHiC0, "; warn =true; }
      if( ac.KeyExists("CCWtHiC1") )     fCCWtHiC1 = ac.GetDouble("CCWtHiC1"); 
      else {errstring+="CCWtHiC1"; warn =true; }

      if(warn) MSG("FarDetEvent",Msg::kWarning) << errstring.Data() <<endl;

      if( totph<fCCWtLowScale )
        {
          MSG("FarDetEvent",Msg::kVerbose) << "  using deweighted constants, low scale... " << endl;
          cal_e = fCCWtLowC0 +
                  fCCWtLowC1*totph + 
                  fCCWtLowC2*totph*totph +
                  fCCWtLowC3*totph*totph*totph +
                  fCCWtLowC4*totph*totph*totph*totph
                   + 0.18*sqrt(totph)*(1.0-totph/5.0)*exp(-totph/10.0);
        }
      else if( totph<fCCWtMidScale )
        {
          MSG("FarDetEvent",Msg::kVerbose) << "  using deweighted constants, middle scale... " << endl;
          cal_e = fCCWtMidC0 + fCCWtMidC1*totph;
        }
      else 
        {
          MSG("FarDetEvent",Msg::kVerbose) << "  using deweighted constants, high scale... " << endl;
          cal_e = fCCWtHiC0 + fCCWtHiC1*totph;
        }

    }

  MSG("FarDetEvent",Msg::kVerbose) << "   " << totph << " MIPs -> " << cal_e << " GeV" << endl;

  return cal_e;
}

void AlgFarDetEvent::Trace(const char* /*c*/ ) const
{

}

        

---