ElectronTools/src/ElectronEPcombinator.cc

0001 #include "EgammaAnalysis/ElectronTools/interface/ElectronEPcombinator.h"
0002 #include <iostream>
0003
0004 //Accessor to the combination results
0005 void ElectronEPcombinator::combine(SimpleElectron& electron) {
0006   electron_ = electron;
0007   computeEPcombination();
0008   electron.setCombinedMomentum(combinedMomentum_);
0009   electron.setCombinedMomentumError(combinedMomentumError_);
0010 }
0011
0012 //Core code to compute the EP combination
0013 void ElectronEPcombinator::computeEPcombination() {
0014   if (mode_ == 1) {
0015     scEnergy_ = electron_.getNewEnergy();
0016     scEnergyError_ = electron_.getNewEnergyError();
0017   }
0018   if (mode_ == 2) {
0019     scEnergy_ = electron_.getRegEnergy();
0020     scEnergyError_ = electron_.getRegEnergyError();
0021   }
0022   trackerMomentum_ = electron_.getTrackerMomentum();
0023   trackerMomentumError_ = electron_.getTrackerMomentumError();
0024   elClass_ = electron_.getElClass();
0025
0026   combinedMomentum_ = scEnergy_;  // initial
0027   combinedMomentumError_ = 999.;
0028
0029   // first check for large errors
0030
0031   if (trackerMomentumError_ / trackerMomentum_ > 0.5 && scEnergyError_ / scEnergy_ <= 0.5) {
0032     combinedMomentum_ = scEnergy_;
0033     combinedMomentumError_ = scEnergyError_;
0034   } else if (trackerMomentumError_ / trackerMomentum_ <= 0.5 && scEnergyError_ / scEnergy_ > 0.5) {
0035     combinedMomentum_ = trackerMomentum_;
0036     combinedMomentumError_ = trackerMomentumError_;
0037   } else if (trackerMomentumError_ / trackerMomentum_ > 0.5 && scEnergyError_ / scEnergy_ > 0.5) {
0038     if (trackerMomentumError_ / trackerMomentum_ < scEnergyError_ / scEnergy_) {
0039       combinedMomentum_ = trackerMomentum_;
0040       combinedMomentumError_ = trackerMomentumError_;
0041     } else {
0042       combinedMomentum_ = scEnergy_;
0043       combinedMomentumError_ = scEnergyError_;
0044     }
0045   }
0046
0047   // then apply the combination algorithm
0048   else {
0049     // calculate E/p and corresponding error
0050     float eOverP = scEnergy_ / trackerMomentum_;
0051     float errorEOverP = sqrt((scEnergyError_ / trackerMomentum_) * (scEnergyError_ / trackerMomentum_) +
0052                              (scEnergy_ * trackerMomentumError_ / trackerMomentum_ / trackerMomentum_) *
0053                                  (scEnergy_ * trackerMomentumError_ / trackerMomentum_ / trackerMomentum_));
0054
0055     bool eleIsNotInCombination = false;
0056     if ((eOverP > 1 + 2.5 * errorEOverP) || (eOverP < 1 - 2.5 * errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3)) {
0057       eleIsNotInCombination = true;
0058     }
0059
0060     if (eleIsNotInCombination) {
0061       if (eOverP > 1) {
0062         combinedMomentum_ = scEnergy_;
0063         combinedMomentumError_ = scEnergyError_;
0064       } else {
0065         if (elClass_ == 0)  // == reco::GsfElectron::GOLDEN)
0066         {
0067           combinedMomentum_ = scEnergy_;
0068           combinedMomentumError_ = scEnergyError_;
0069         }
0070         if (elClass_ == 1)  //reco::GsfElectron::BIGBREM)
0071         {
0072           if (scEnergy_ < 36) {
0073             combinedMomentum_ = trackerMomentum_;
0074             combinedMomentumError_ = trackerMomentumError_;
0075           } else {
0076             combinedMomentum_ = scEnergy_;
0077             combinedMomentumError_ = scEnergyError_;
0078           }
0079         }
0080         if (elClass_ == 2)  // == reco::GsfElectron::BADTRACK)
0081         {
0082           combinedMomentum_ = scEnergy_;
0083           combinedMomentumError_ = scEnergyError_;
0084         }
0085         if (elClass_ == 3)  //reco::GsfElectron::SHOWERING)
0086         {
0087           if (scEnergy_ < 30) {
0088             combinedMomentum_ = trackerMomentum_;
0089             combinedMomentumError_ = trackerMomentumError_;
0090           } else {
0091             combinedMomentum_ = scEnergy_;
0092             combinedMomentumError_ = scEnergyError_;
0093           }
0094         }
0095         if (elClass_ == 4)  //reco::GsfElectron::GAP)
0096         {
0097           if (scEnergy_ < 60) {
0098             combinedMomentum_ = trackerMomentum_;
0099             combinedMomentumError_ = trackerMomentumError_;
0100           } else {
0101             combinedMomentum_ = scEnergy_;
0102             combinedMomentumError_ = scEnergyError_;
0103           }
0104         }
0105       }
0106     } else {
0107       // combination
0108       combinedMomentum_ = (scEnergy_ / scEnergyError_ / scEnergyError_ +
0109                            trackerMomentum_ / trackerMomentumError_ / trackerMomentumError_) /
0110                           (1 / scEnergyError_ / scEnergyError_ + 1 / trackerMomentumError_ / trackerMomentumError_);
0111       float combinedMomentum_Variance =
0112           1 / (1 / scEnergyError_ / scEnergyError_ + 1 / trackerMomentumError_ / trackerMomentumError_);
0113       combinedMomentumError_ = sqrt(combinedMomentum_Variance);
0114     }
0115   }
0116 }