EnergyResolution.cxx

Go to the documentation of this file.

#include "EnergyResolution.h"  

#include <cmath>

// neutrino energy resolution
// ==========================
Double_t EnergyResolution::GetResolution(Double_t trkmomentum, Double_t trkmomentumerr, Double_t shwenergy, Double_t shwenergyerr, Bool_t rangecurve, Bool_t containedvertex, Detector::Detector_t detector, ReleaseType::Release_t release)
{
  Double_t muon_resolution = MuonResolution(trkmomentum, trkmomentumerr,
                                            rangecurve, containedvertex, detector, release );

  Double_t shower_resolution = ShowerResolution(shwenergy, shwenergyerr,
                                                detector, release );

  return sqrt( muon_resolution*muon_resolution + shower_resolution*shower_resolution );
}

// shower energy resolution
// ========================
Double_t EnergyResolution::ShowerResolution(Double_t shwenergy, Double_t shwenergyerr, Detector::Detector_t detector, ReleaseType::Release_t release)
{
  // dogwood resolution
  if( (release&ReleaseType::kDogwood)==(ReleaseType::kDogwood) ){
    if( detector == Detector::kNear ) return ShowerResolutionDogwoodNear(shwenergy, shwenergyerr);
    if( detector == Detector::kFar ) return ShowerResolutionDogwoodFar(shwenergy, shwenergyerr);
  }

  // cedar resolution
  else if( (release&ReleaseType::kCedar)==(ReleaseType::kCedar) ){
    if( detector == Detector::kNear ) return ShowerResolutionCedar(shwenergy, shwenergyerr);
    if( detector == Detector::kFar ) return ShowerResolutionCedar(shwenergy, shwenergyerr);
  }
  
  // default to zero
  return 0.0;
}

// muon momentum resolution
// ========================
Double_t EnergyResolution::MuonResolution(Double_t trkmomentum, Double_t trkmomentumerr, Bool_t rangecurve, Bool_t containedvertex, Detector::Detector_t detector, ReleaseType::Release_t release)
{
  // dogwood resolution
  if( (release&ReleaseType::kDogwood)==(ReleaseType::kDogwood) ){
    if( detector == Detector::kNear ) return MuonResolutionDogwoodNear(trkmomentum,trkmomentumerr,
                                                                       rangecurve,containedvertex);
    if( detector == Detector::kFar ) return MuonResolutionDogwoodFar(trkmomentum,trkmomentumerr,
                                                                     rangecurve,containedvertex);
  }

  // cedar resolution
  else if( (release&ReleaseType::kCedar)==(ReleaseType::kCedar) ){
    if( detector == Detector::kNear ) return MuonResolutionCedar(trkmomentum,trkmomentumerr,
                                                                 rangecurve,containedvertex);
    if( detector == Detector::kFar ) return MuonResolutionCedar(trkmomentum,trkmomentumerr,
                                                                rangecurve,containedvertex);
  }

  // default to zero
  return 0.0;
}
  

// muon momentum resolution (cedar/dogwood, near/far)
// ==================================================

Double_t EnergyResolution::MuonResolutionDogwoodNear( Double_t trkmomentum, Double_t trkmomentumerr, Bool_t rangecurve, Bool_t containedvertex )
{
  return MuonResolutionCedar(trkmomentum, trkmomentumerr, rangecurve, containedvertex);
}
  
Double_t EnergyResolution::MuonResolutionDogwoodFar( Double_t trkmomentum, Double_t trkmomentumerr, Bool_t rangecurve, Bool_t containedvertex )
{
  return MuonResolutionCedar(trkmomentum, trkmomentumerr, rangecurve, containedvertex);
}
  
Double_t EnergyResolution::MuonResolutionCedar( Double_t trkmomentum, Double_t trkmomentumerr, Bool_t rangecurve, Bool_t containedvertex )
{
  if( containedvertex ){

    // momentum from range
    if( rangecurve ){
      Double_t trkmomentum_input = trkmomentum;
      if( trkmomentum_input<0.0 ) trkmomentum_input = 0.0;

      Double_t fRange0 = 0.0;    // Emu = 0  : sigma   -> 0 GeV
      Double_t fRange1 = 0.0467;
      Double_t fRange2 = 0.0700; // Emu >> 0 : sigma/E -> 7%

      Double_t rangesigma0 = fRange0;
      Double_t rangesigma1 = fRange1*sqrt(trkmomentum_input);
      Double_t rangesigma2 = fRange2*trkmomentum_input;

      return sqrt( rangesigma0*rangesigma0 + rangesigma1*rangesigma1 + rangesigma2*rangesigma2 );
    }

    // momentum from curvature
    else{
      Double_t trkmomentumerr_input = trkmomentumerr;
      if( trkmomentumerr_input<0.0 ) trkmomentumerr_input = 0.0;

      Double_t fCurve0 = 0.0;
      Double_t fCurve1 = 1.237;
      Double_t fCurve2 = 1.198;

      Double_t curvesigma0 = fCurve0;
      Double_t curvesigma1 = fCurve1*sqrt(trkmomentumerr_input);
      Double_t curvesigma2 = fCurve2*trkmomentumerr_input;

      return sqrt( curvesigma0*curvesigma0 + curvesigma1*curvesigma1 + curvesigma2*curvesigma2 );
    }
  }

  // default to zero
  return 0.0;
}


// shower energy resolution (cedar/dogwood, near/far)
// ==================================================

Double_t EnergyResolution::ShowerResolutionDogwoodNear( Double_t shwenergy, Double_t shwenergyerr )
{
  return ShowerResolutionCedar(shwenergy, shwenergyerr);
}

Double_t EnergyResolution::ShowerResolutionDogwoodFar( Double_t shwenergy, Double_t shwenergyerr )
{
  return ShowerResolutionCedar(shwenergy, shwenergyerr);
}

Double_t EnergyResolution::ShowerResolutionCedar( Double_t shwenergy, Double_t shwenergyerr )
{
  Double_t shwenergy_input = shwenergy;
  if( shwenergy_input<0.0 ) shwenergy_input = 0.0;

  Double_t shwenergyerr_input = shwenergyerr;
  if( shwenergyerr_input<0.0 ) shwenergyerr_input = 0.0;

  Double_t fShw0 = 0.190;  // Eshw = 0  : sigma   -> 190 MeV
  Double_t fShw1 = 0.509;  // Eshw >> 0 : sigma/E -> 51%/sqrt(E) 
  Double_t fShw2 = 0.0;

  Double_t shwsigma0 = fShw0;
  Double_t shwsigma1 = fShw1*sqrt(shwenergy_input);
  Double_t shwsigma2 = fShw2*shwenergy_input;

  return sqrt( shwsigma0*shwsigma0 + shwsigma1*shwsigma1 + shwsigma2*shwsigma2 );
}  

---