CMS 3D CMS Logo

Public Member Functions | Private Member Functions

ElectronMomentumCorrector Class Reference

#include <ElectronMomentumCorrector.h>

List of all members.

Public Member Functions

void correct (reco::GsfElectron &, TrajectoryStateOnSurface &)
 ElectronMomentumCorrector ()
 ~ElectronMomentumCorrector ()

Private Member Functions

float energyError (float E, float *par) const

Detailed Description

Definition at line 19 of file ElectronMomentumCorrector.h.


Constructor & Destructor Documentation

ElectronMomentumCorrector::ElectronMomentumCorrector ( ) [inline]

Definition at line 25 of file ElectronMomentumCorrector.h.

{}
ElectronMomentumCorrector::~ElectronMomentumCorrector ( ) [inline]

Definition at line 27 of file ElectronMomentumCorrector.h.

{}

Member Function Documentation

void ElectronMomentumCorrector::correct ( reco::GsfElectron electron,
TrajectoryStateOnSurface vtxTsos 
)

Definition at line 33 of file ElectronMomentumCorrector.cc.

References reco::GsfElectron::classification(), reco::GsfElectron::correctMomentum(), reco::GsfElectron::ecalEnergy(), reco::GsfElectron::ecalEnergyError(), reco::GsfElectron::GAP, reco::GsfElectron::GOLDEN, reco::GsfElectron::isEB(), reco::GsfElectron::isEE(), reco::GsfElectron::isGap(), GaussianSumUtilities1D::mode(), MultiGaussianStateTransform::multiState1D(), reco::GsfElectron::p4(), reco::GsfElectron::P4_COMBINATION, reco::GsfElectron::p4Error(), reco::GsfElectron::SHOWERING, mathSSE::sqrt(), reco::btau::trackMomentum, reco::GsfElectron::trackMomentumAtVtx(), reco::GsfElectron::UNKNOWN, and SingleGaussianState1D::variance().

Referenced by GsfElectronAlgo::createElectron().

 {
  if (electron.p4Error(reco::GsfElectron::P4_COMBINATION)!= 999.)
   {
    edm::LogWarning("ElectronMomentumCorrector::correct")<<"already done" ;
    return ;
   }

  math::XYZTLorentzVector newMomentum = electron.p4() ; // default
  int elClass = electron.classification() ;

  // irrelevant classification
  if ( (elClass <= reco::GsfElectron::UNKNOWN) ||
           (elClass>reco::GsfElectron::GAP) )
   {
    edm::LogWarning("ElectronMomentumCorrector::correct")<<"unexpected classification" ;
    return ;
   }

  float scEnergy = electron.ecalEnergy() ;
  float errorEnergy_ = electron.ecalEnergyError() ;

  float trackMomentum  = electron.trackMomentumAtVtx().R() ;
  float errorTrackMomentum = 999. ;

  // retreive momentum error
  MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(vtxTsos,0));
  GaussianSumUtilities1D qpUtils(qpState);
  errorTrackMomentum = trackMomentum*trackMomentum*sqrt(qpUtils.mode().variance());

  float finalMomentum = electron.p4().t(); // initial
  float finalMomentumError = 999.;


  // first check for large errors

  if (errorTrackMomentum/trackMomentum > 0.5 && errorEnergy_/scEnergy <= 0.5) {
    finalMomentum = scEnergy;    finalMomentumError = errorEnergy_;
   }
  else if (errorTrackMomentum/trackMomentum <= 0.5 && errorEnergy_/scEnergy > 0.5){
    finalMomentum = trackMomentum;  finalMomentumError = errorTrackMomentum;
   }
  else if (errorTrackMomentum/trackMomentum > 0.5 && errorEnergy_/scEnergy > 0.5){
    if (errorTrackMomentum/trackMomentum < errorEnergy_/scEnergy) {
      finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum;
     }
    else{
      finalMomentum = scEnergy; finalMomentumError = errorEnergy_;
     }
  }

  // then apply the combination algorithm
  else {

     // calculate E/p and corresponding error
    float eOverP = scEnergy / trackMomentum;
    float errorEOverP = sqrt(
                             (errorEnergy_/trackMomentum)*(errorEnergy_/trackMomentum) +
                             (scEnergy*errorTrackMomentum/trackMomentum/trackMomentum)*
                             (scEnergy*errorTrackMomentum/trackMomentum/trackMomentum));

    if ( eOverP  > 1 + 2.5*errorEOverP )
     {
      finalMomentum = scEnergy; finalMomentumError = errorEnergy_;
      if ((elClass==reco::GsfElectron::GOLDEN) && electron.isEB() && (eOverP<1.15))
        {
          if (scEnergy<15) {finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum;}
        }
     }
    else if ( eOverP < 1 - 2.5*errorEOverP )
     {
      finalMomentum = scEnergy; finalMomentumError = errorEnergy_;
      if (elClass==reco::GsfElectron::SHOWERING)
       {
              if (electron.isEB())
               {
                      if (scEnergy<18)
                       { finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum; }
               }
              else if (electron.isEE())
               {
          if (scEnergy<13)
           {finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum;}
               }
              else
               { edm::LogWarning("ElectronMomentumCorrector::correct")<<"nor barrel neither endcap electron ?!" ; }
             }
            else if (electron.isGap())
             {
              if (scEnergy<60)
               { finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum ; }
             }
      }
    else
     {
      // combination
      finalMomentum = (scEnergy/errorEnergy_/errorEnergy_ + trackMomentum/errorTrackMomentum/errorTrackMomentum) /
        (1/errorEnergy_/errorEnergy_ + 1/errorTrackMomentum/errorTrackMomentum);
      float finalMomentumVariance = 1 / (1/errorEnergy_/errorEnergy_ + 1/errorTrackMomentum/errorTrackMomentum);
      finalMomentumError = sqrt(finalMomentumVariance);
     }

  }

  math::XYZTLorentzVector oldMomentum = electron.p4() ;
  newMomentum = math::XYZTLorentzVector
   ( oldMomentum.x()*finalMomentum/oldMomentum.t(),
     oldMomentum.y()*finalMomentum/oldMomentum.t(),
     oldMomentum.z()*finalMomentum/oldMomentum.t(),
     finalMomentum ) ;


  // final set
  electron.correctMomentum(newMomentum,errorTrackMomentum,finalMomentumError);

 }
float ElectronMomentumCorrector::energyError ( float  E,
float *  par 
) const [private]