EgammaElectronAlgos/src/ElectronMomentumCorrector.cc

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronMomentumCorrector.h"
0002
0003 #include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
0004 #include "TrackingTools/GsfTools/interface/GaussianSumUtilities1D.h"
0005 #include "TrackingTools/GsfTools/interface/MultiGaussianStateTransform.h"
0006
0007 /****************************************************************************
0008  *
0009  * Class based E-p combination for the final electron momentum. It relies on
0010  * the electron classification and on the dtermination of track momentum and ecal
0011  * supercluster energy errors. The track momentum error is taken from the gsf fit.
0012  * The ecal supercluster energy error is taken from a class dependant parametrisation
0013  * of the energy resolution.
0014  *
0015  *
0016  * \author Federico Ferri - INFN Milano, Bicocca university
0017  * \author Ivica Puljak - FESB, Split
0018  * \author Stephanie Baffioni - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3
0019  *
0020  *
0021  ****************************************************************************/
0022
0023 egamma::ElectronMomentum egamma::correctElectronMomentum(reco::GsfElectron const& electron,
0024                                                          TrajectoryStateOnSurface const& vtxTsos) {
0025   int elClass = electron.classification();
0026
0027   //=======================================================================================
0028   // cluster energy
0029   //=======================================================================================
0030
0031   float scEnergy = electron.correctedEcalEnergy();
0032   float errorEnergy = electron.correctedEcalEnergyError();
0033
0034   //=======================================================================================
0035   // track  momentum
0036   //=======================================================================================
0037
0038   // basic values
0039   float trackMomentum = electron.trackMomentumAtVtx().R();
0040   //float errorTrackMomentum = 999. ;
0041
0042   // tracker momentum scale corrections (Mykhailo Dalchenko)
0043   double scale = 1.;
0044   if (electron.isEB()) {
0045     if (elClass == 0) {
0046       scale = 1. / (0.00104 * sqrt(trackMomentum) + 1);
0047     }
0048     if (elClass == 1) {
0049       scale = 1. / (0.0017 * sqrt(trackMomentum) + 0.9986);
0050     }
0051     if (elClass == 3) {
0052       scale = 1. / (1.004 - 0.00021 * trackMomentum);
0053     }
0054     if (elClass == 4) {
0055       scale = 0.995;
0056     }
0057   } else if (electron.isEE()) {
0058     if (elClass == 3) {
0059       scale = 1. / (1.01432 - 0.00201872 * trackMomentum + 0.0000142621 * trackMomentum * trackMomentum);
0060     }
0061     if (elClass == 4) {
0062       scale = 1. / (0.996859 - 0.000345347 * trackMomentum);
0063     }
0064   }
0065   if (scale < 0.)
0066     scale = 1.;  // CC added protection
0067   trackMomentum = trackMomentum * scale;
0068
0069   // error (must be done after trackMomentum rescaling)
0070   MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(vtxTsos, 0));
0071   GaussianSumUtilities1D qpUtils(qpState);
0072   float errorTrackMomentum = trackMomentum * trackMomentum * sqrt(qpUtils.mode().variance());
0073
0074   //=======================================================================================
0075   // combination
0076   //=======================================================================================
0077
0078   float finalMomentum = electron.p4().t();  // initial
0079   float finalMomentumError = 999.;
0080
0081   // first check for large errors
0082
0083   if (errorTrackMomentum / trackMomentum > 0.5 && errorEnergy / scEnergy <= 0.5) {
0084     finalMomentum = scEnergy;
0085     finalMomentumError = errorEnergy;
0086   } else if (errorTrackMomentum / trackMomentum <= 0.5 && errorEnergy / scEnergy > 0.5) {
0087     finalMomentum = trackMomentum;
0088     finalMomentumError = errorTrackMomentum;
0089   } else if (errorTrackMomentum / trackMomentum > 0.5 && errorEnergy / scEnergy > 0.5) {
0090     if (errorTrackMomentum / trackMomentum < errorEnergy / scEnergy) {
0091       finalMomentum = trackMomentum;
0092       finalMomentumError = errorTrackMomentum;
0093     } else {
0094       finalMomentum = scEnergy;
0095       finalMomentumError = errorEnergy;
0096     }
0097   }
0098
0099   // then apply the combination algorithm
0100   else {
0101     // calculate E/p and corresponding error
0102     float eOverP = scEnergy / trackMomentum;
0103     float errorEOverP = sqrt((errorEnergy / trackMomentum) * (errorEnergy / trackMomentum) +
0104                              (scEnergy * errorTrackMomentum / trackMomentum / trackMomentum) *
0105                                  (scEnergy * errorTrackMomentum / trackMomentum / trackMomentum));
0106
0107     bool eleIsNotInCombination = false;
0108     if ((eOverP > 1 + 2.5 * errorEOverP) || (eOverP < 1 - 2.5 * errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3)) {
0109       eleIsNotInCombination = true;
0110     }
0111     if (eleIsNotInCombination) {
0112       if (eOverP > 1) {
0113         finalMomentum = scEnergy;
0114         finalMomentumError = errorEnergy;
0115       } else {
0116         if (elClass == reco::GsfElectron::GOLDEN) {
0117           finalMomentum = scEnergy;
0118           finalMomentumError = errorEnergy;
0119         }
0120         if (elClass == reco::GsfElectron::BIGBREM) {
0121           if (scEnergy < 36) {
0122             finalMomentum = trackMomentum;
0123             finalMomentumError = errorTrackMomentum;
0124           } else {
0125             finalMomentum = scEnergy;
0126             finalMomentumError = errorEnergy;
0127           }
0128         }
0129         if (elClass == reco::GsfElectron::BADTRACK) {
0130           finalMomentum = scEnergy;
0131           finalMomentumError = errorEnergy;
0132         }
0133         if (elClass == reco::GsfElectron::SHOWERING) {
0134           if (scEnergy < 30) {
0135             finalMomentum = trackMomentum;
0136             finalMomentumError = errorTrackMomentum;
0137           } else {
0138             finalMomentum = scEnergy;
0139             finalMomentumError = errorEnergy;
0140           }
0141         }
0142         if (elClass == reco::GsfElectron::GAP) {
0143           if (scEnergy < 60) {
0144             finalMomentum = trackMomentum;
0145             finalMomentumError = errorTrackMomentum;
0146           } else {
0147             finalMomentum = scEnergy;
0148             finalMomentumError = errorEnergy;
0149           }
0150         }
0151       }
0152     }
0153
0154     else {
0155       // combination
0156       finalMomentum = (scEnergy / errorEnergy / errorEnergy + trackMomentum / errorTrackMomentum / errorTrackMomentum) /
0157                       (1 / errorEnergy / errorEnergy + 1 / errorTrackMomentum / errorTrackMomentum);
0158       float finalMomentumVariance = 1 / (1 / errorEnergy / errorEnergy + 1 / errorTrackMomentum / errorTrackMomentum);
0159       finalMomentumError = sqrt(finalMomentumVariance);
0160     }
0161   }
0162
0163   //=======================================================================================
0164   // final set
0165   //=======================================================================================
0166
0167   auto const& oldMomentum = electron.p4();
0168
0169   return {{oldMomentum.x() * finalMomentum / oldMomentum.t(),
0170            oldMomentum.y() * finalMomentum / oldMomentum.t(),
0171            oldMomentum.z() * finalMomentum / oldMomentum.t(),
0172            finalMomentum},
0173           errorTrackMomentum,
0174           finalMomentumError};
0175 }