MCNNAnalysis/EnergyCorrections.cxx

Go to the documentation of this file.

//
// EnergyCorrections.h
//
// Corrections to energy based quantities
//
// Created:  M. Kordosky -- December, 2005
//
// $Author: rhatcher $ 
//
// $Revision: 1.2 $
// 
// $Name:  $
//
// $Id: EnergyCorrections.cxx,v 1.2 2009/05/26 22:04:51 rhatcher Exp $
//
// $Log: EnergyCorrections.cxx,v $
// Revision 1.2  2009/05/26 22:04:51  rhatcher
// Use Detector:: rather than DetectorType:: (synonym) everywhere, so we can
// someday eliminate the deprecated DetectorType::.
//
// Add #include <cmath> in enough places sufficient for resolving random
// use of fabs() when compiled w/ gcc 4.x
//
// Revision 1.1  2009/03/23 09:28:27  rbpatter
// initial import of new singleton MCNN package
//
// Revision 1.4  2006/02/15 15:57:29  kordosky
// add correction from for momentum via curvature as advocated in 1430-v2. Only modifies momenta of positively charged tracks.
//
// Revision 1.3  2006/02/13 00:20:02  kordosky
// correction of +1.8% to FD data shower energy
//
// Revision 1.2  2006/02/12 22:09:58  kordosky
// Update momentum correction to account for new density and thickness measurements.
//
// Revision 1.1  2005/12/16 03:38:51  kordosky
// place implementation in a cxx file
//
// Revision 1.4  2005/12/16 02:52:43  kordosky
// Forgot the include guards.  Thansk Niki.
//
// Revision 1.3  2005/12/15 22:44:31  kordosky
// retweaked A Culling formulae in reaction to his presentation this morning
//
// Revision 1.2  2005/12/15 22:26:19  kordosky
// E=sqrt(p*p+m*m)  != p   (almost negligible)
//
// Revision 1.1  2005/12/15 21:57:09  kordosky
// moved from Mad
//
// Revision 1.2  2005/12/14 20:35:59  kordosky
// Updated energy corrections.
//
// Revision 1.1  2005/12/13 19:41:35  kordosky
// Standalone routines to correct the momentum via range to to groom et al range tables. Also, correct data to 7.755 g/cc which is our current best estimate of the detector density.
//

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

float CorrectMomentumFromRange(float p, bool isdata, Detector::Detector_t det){
  static const float c[4]={1.01334,0.05563,-0.05346,0.01205};
  
  // correction for difference in data mc steel density
  if(isdata){
    // inital correction, pre-Oxford 2006
    //static const float dcor=7.755/7.87;// data/mc density
    float dcor=1;
    if (det==Detector::kNear) dcor=(7.85*2.563)/(7.87*2.54); 
    else if(det==Detector::kFar) dcor=(7.85*2.558)/(7.87*2.54);

    p*=dcor;
  }
  // 
  float pcor=p/(c[0] + c[1]*log(p) + c[2]*pow(log(p), 2) + c[3]*pow(log(p),3));
  return pcor;
}

float CorrectSignedMomentumFromCurvature(float p, bool /* isdata */, Detector::Detector_t /* det */){
  // input is the signed momentum (e.g. p/q)
  // isdata and det are not used... but maybe in the future
  float pcor=p;
  if(pcor!=0) {
    // correction advertised in 1430-v2, J. Musser
    // note: for 1/p < 0   C=1, so correction only matters for mu+
    float C = (1.01+0.1*fabs(1/p))/(1.01-0.1*(1/p));
    pcor*=(1.0/C);
  }
  return pcor;
}


float CorrectEnergyFromRange(float E, bool isdata, Detector::Detector_t det){
   const float m=0.1057;// mon mass
   float p = sqrt(E*E -m*m);
   float pcor = CorrectMomentumFromRange(p,isdata,det);
   return sqrt(pcor*pcor +m*m);
} 


float CorrectShowerEnergyNear(float E, const CandShowerHandle::ShowerType_t& st, int mode, bool /* isdata */){

   float ecor=E;
   if(st==CandShowerHandle::kCC){
      if(mode==1){
         // Niki Correction
         ecor=E/1.18;
      }
      else if(mode==2){
         // Andy Correction
         ecor=((E/1.05)*(1.-0.35*exp(-0.18*E/1.06)));
      }
   }
   else if(st==CandShowerHandle::kWtCC){
      if(mode==1){
         // Niki Correction
         ecor=((E)*(1.+0.50*exp(-1.00*E)));
      }
      else if(mode==2){
         // Andy Correction
         ecor=E/1.03;
      }
   }
   return ecor;
}

float CorrectShowerEnergyFar(float E, const CandShowerHandle::ShowerType_t& st, int mode, bool isdata){

   float ecor=E;
   if(st==CandShowerHandle::kCC){
     if(isdata) {
       // a correction for the FD MIP scale
       // for the beam data, one measured MIP 
       // actually corresponds to 1.018 MIPs
       // so we must correct shower energy up
       const float mip_scale_correction=1.018;
       E*=mip_scale_correction;
     }

     if(mode==1){
       // Niki Correction
       ecor=((E)*(1.-0.12*exp(-0.12*E)));
     }
     else if(mode==2){
       // Andy Correction
       //        ecor=(E)*(1.-0.2*exp(-0.2*E));
       ecor=E*(0.99-0.035*E*exp(-0.25*E));
     }
   }
   else if(st==CandShowerHandle::kWtCC){
      if(mode==1){
         // Niki Correction
         ecor=((E)*(1.+0.18*exp(-0.35*E)));
      }
      else if(mode==2){
         // Andy Correction
         E=ecor;
      }
   }
   return ecor;
}
float CorrectShowerEnergy(float E, const Detector::Detector_t& det, 
                          const CandShowerHandle::ShowerType_t& st, 
                          int mode, bool isdata){
  float ecor=E;
   if(det==Detector::kNear){
      ecor = CorrectShowerEnergyNear(E,st,mode,isdata);
   }
   else if(det==Detector::kFar){
      ecor = CorrectShowerEnergyFar(E,st,mode,isdata);
   }

   return ecor;
   
}

---