ElectronEPcombinator.cc

Go to the documentation of this file.

#include "EgammaAnalysis/ElectronTools/interface/ElectronEPcombinator.h"
#include <iostream>

//Accessor to the combination results
void ElectronEPcombinator::combine(SimpleElectron& electron) {
  electron_ = electron;
  computeEPcombination();
  electron.setCombinedMomentum(combinedMomentum_);
  electron.setCombinedMomentumError(combinedMomentumError_);
}

//Core code to compute the EP combination
void ElectronEPcombinator::computeEPcombination() {
  if (mode_ == 1) {
    scEnergy_ = electron_.getNewEnergy();
    scEnergyError_ = electron_.getNewEnergyError();
  }
  if (mode_ == 2) {
    scEnergy_ = electron_.getRegEnergy();
    scEnergyError_ = electron_.getRegEnergyError();
  }
  trackerMomentum_ = electron_.getTrackerMomentum();
  trackerMomentumError_ = electron_.getTrackerMomentumError();
  elClass_ = electron_.getElClass();

  combinedMomentum_ = scEnergy_;  // initial
  combinedMomentumError_ = 999.;

  // first check for large errors

  if (trackerMomentumError_ / trackerMomentum_ > 0.5 && scEnergyError_ / scEnergy_ <= 0.5) {
    combinedMomentum_ = scEnergy_;
    combinedMomentumError_ = scEnergyError_;
  } else if (trackerMomentumError_ / trackerMomentum_ <= 0.5 && scEnergyError_ / scEnergy_ > 0.5) {
    combinedMomentum_ = trackerMomentum_;
    combinedMomentumError_ = trackerMomentumError_;
  } else if (trackerMomentumError_ / trackerMomentum_ > 0.5 && scEnergyError_ / scEnergy_ > 0.5) {
    if (trackerMomentumError_ / trackerMomentum_ < scEnergyError_ / scEnergy_) {
      combinedMomentum_ = trackerMomentum_;
      combinedMomentumError_ = trackerMomentumError_;
    } else {
      combinedMomentum_ = scEnergy_;
      combinedMomentumError_ = scEnergyError_;
    }
  }

  // then apply the combination algorithm
  else {
    // calculate E/p and corresponding error
    float eOverP = scEnergy_ / trackerMomentum_;
    float errorEOverP = sqrt((scEnergyError_ / trackerMomentum_) * (scEnergyError_ / trackerMomentum_) +
                             (scEnergy_ * trackerMomentumError_ / trackerMomentum_ / trackerMomentum_) *
                                 (scEnergy_ * trackerMomentumError_ / trackerMomentum_ / trackerMomentum_));

    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) {
        combinedMomentum_ = scEnergy_;
        combinedMomentumError_ = scEnergyError_;
      } else {
        if (elClass_ == 0)  // == reco::GsfElectron::GOLDEN)
        {
          combinedMomentum_ = scEnergy_;
          combinedMomentumError_ = scEnergyError_;
        }
        if (elClass_ == 1)  //reco::GsfElectron::BIGBREM)
        {
          if (scEnergy_ < 36) {
            combinedMomentum_ = trackerMomentum_;
            combinedMomentumError_ = trackerMomentumError_;
          } else {
            combinedMomentum_ = scEnergy_;
            combinedMomentumError_ = scEnergyError_;
          }
        }
        if (elClass_ == 2)  // == reco::GsfElectron::BADTRACK)
        {
          combinedMomentum_ = scEnergy_;
          combinedMomentumError_ = scEnergyError_;
        }
        if (elClass_ == 3)  //reco::GsfElectron::SHOWERING)
        {
          if (scEnergy_ < 30) {
            combinedMomentum_ = trackerMomentum_;
            combinedMomentumError_ = trackerMomentumError_;
          } else {
            combinedMomentum_ = scEnergy_;
            combinedMomentumError_ = scEnergyError_;
          }
        }
        if (elClass_ == 4)  //reco::GsfElectron::GAP)
        {
          if (scEnergy_ < 60) {
            combinedMomentum_ = trackerMomentum_;
            combinedMomentumError_ = trackerMomentumError_;
          } else {
            combinedMomentum_ = scEnergy_;
            combinedMomentumError_ = scEnergyError_;
          }
        }
      }
    } else {
      // combination
      combinedMomentum_ = (scEnergy_ / scEnergyError_ / scEnergyError_ +
                           trackerMomentum_ / trackerMomentumError_ / trackerMomentumError_) /
                          (1 / scEnergyError_ / scEnergyError_ + 1 / trackerMomentumError_ / trackerMomentumError_);
      float combinedMomentum_Variance =
          1 / (1 / scEnergyError_ / scEnergyError_ + 1 / trackerMomentumError_ / trackerMomentumError_);
      combinedMomentumError_ = sqrt(combinedMomentum_Variance);
    }
  }
}