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 29 of file ElectronMomentumCorrector.cc.

References reco::GsfElectron::BADTRACK, reco::GsfElectron::BIGBREM, reco::GsfElectron::classification(), reco::GsfElectron::correctedEcalEnergy(), reco::GsfElectron::correctedEcalEnergyError(), reco::GsfElectron::correctMomentum(), funct::false, reco::GsfElectron::GAP, reco::GsfElectron::GOLDEN, reco::GsfElectron::isEB(), reco::GsfElectron::isEE(), GaussianSumUtilities1D::mode(), MultiGaussianStateTransform::multiState1D(), reco::GsfElectron::p4(), reco::GsfElectron::P4_COMBINATION, reco::GsfElectron::p4Error(), pileupReCalc_HLTpaths::scale, reco::GsfElectron::SHOWERING, mathSSE::sqrt(), reco::btau::trackMomentum, reco::GsfElectron::trackMomentumAtVtx(), funct::true, 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 ;
   }


  //=======================================================================================
  // cluster energy
  //=======================================================================================

  float scEnergy = electron.correctedEcalEnergy() ;
  float errorEnergy = electron.correctedEcalEnergyError() ;


  //=======================================================================================
  // track  momentum
  //=======================================================================================

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

  // tracker momentum scale corrections (Mykhailo Dalchenko)
  double scale=1.;
  if (electron.isEB()){
    if (elClass==0) {scale = 1./(0.00104*sqrt(trackMomentum)+1);}
    if (elClass==1) {scale = 1./(0.0017*sqrt(trackMomentum)+0.9986);}
    if (elClass==3) {scale = 1./(1.004 - 0.00021*trackMomentum);}
    if (elClass==4) {scale = 0.995;}
  } else if (electron.isEE()){
    if (elClass==3){scale = 1./(1.01432-0.00201872*trackMomentum+0.0000142621*trackMomentum*trackMomentum);}
    if (elClass==4){scale = 1./(0.996859-0.000345347*trackMomentum);}
  }
  if (scale<0.) scale = 1.;  // CC added protection
  trackMomentum = trackMomentum*scale ;

  // error (must be done after trackMomentum rescaling)
  MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(vtxTsos,0)) ;
  GaussianSumUtilities1D qpUtils(qpState) ;
  float errorTrackMomentum = trackMomentum*trackMomentum*sqrt(qpUtils.mode().variance()) ;


  //=======================================================================================
  // combination
  //=======================================================================================

  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 ; }
//           }
//      }

    bool eleIsNotInCombination = false ;
    if ( (eOverP  > 1 + 2.5*errorEOverP) || (eOverP  < 1 - 2.5*errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3) )
     { eleIsNotInCombination = true ; }
    if (eleIsNotInCombination)
     {
      if (eOverP > 1)
       { finalMomentum = scEnergy ; finalMomentumError = errorEnergy ; }
      else
       {
        if (elClass == reco::GsfElectron::GOLDEN)
         { finalMomentum = scEnergy; finalMomentumError = errorEnergy; }
        if (elClass == reco::GsfElectron::BIGBREM)
         {
          if (scEnergy<36)
           { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum ; }
          else
           { finalMomentum = scEnergy ; finalMomentumError = errorEnergy ; }
         }
        if (elClass == reco::GsfElectron::BADTRACK)
         { finalMomentum = scEnergy; finalMomentumError = errorEnergy ; }
        if (elClass == reco::GsfElectron::SHOWERING)
         {
          if (scEnergy<30)
           { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum; }
          else
           { finalMomentum = scEnergy; finalMomentumError = errorEnergy;}
         }
        if (elClass == reco::GsfElectron::GAP)
         {
          if (scEnergy<60)
           { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum ; }
          else
           { finalMomentum = scEnergy; finalMomentumError = errorEnergy ; }
         }
       }
     }

    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);
     }

  }


  //=======================================================================================
  // final set
  //=======================================================================================

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

  electron.correctMomentum(newMomentum,errorTrackMomentum,finalMomentumError);

 }

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