ElectronTools/src/EpCombinationTool.cc

0001 #include "EgammaAnalysis/ElectronTools/interface/EpCombinationTool.h"
0002 #include "CondFormats/GBRForest/interface/GBRForest.h"
0003 #include <TFile.h>
0004 #include <TSystem.h>
0005 #include <cmath>
0006 #include <vector>
0007 #include <iostream>
0008
0009 using namespace std;
0010
0011 /*****************************************************************/
0012 EpCombinationTool::EpCombinationTool()
0013     : m_forest(nullptr),
0014       m_ownForest(false)
0015 /*****************************************************************/
0016 {}
0017
0018 /*****************************************************************/
0019 EpCombinationTool::~EpCombinationTool()
0020 /*****************************************************************/
0021 {
0022   if (m_ownForest)
0023     delete m_forest;
0024 }
0025
0026 /*****************************************************************/
0027 bool EpCombinationTool::init(const std::string& regressionFileName, const std::string& bdtName)
0028 /*****************************************************************/
0029 {
0030   TFile* regressionFile = TFile::Open(regressionFileName.c_str());
0031   if (!regressionFile) {
0032     cout << "ERROR: Cannot open regression file " << regressionFileName << "\n";
0033     return false;
0034   }
0035   if (m_ownForest)
0036     delete m_forest;
0037   m_forest = (GBRForest*)regressionFile->Get(bdtName.c_str());
0038   m_ownForest = true;
0039   //regressionFile->GetObject(bdtName.c_str(), m_forest);
0040   if (!m_forest) {
0041     cout << "ERROR: Cannot find forest " << bdtName << " in " << regressionFileName << "\n";
0042     regressionFile->Close();
0043     return false;
0044   }
0045   regressionFile->Close();
0046   return true;
0047 }
0048
0049 bool EpCombinationTool::init(const GBRForest* forest) {
0050   if (m_ownForest)
0051     delete m_forest;
0052   m_forest = forest;
0053   m_ownForest = false;
0054   return true;
0055 }
0056
0057 /*****************************************************************/
0058 void EpCombinationTool::combine(SimpleElectron& mySimpleElectron) const
0059 /*****************************************************************/
0060 {
0061   if (!m_forest) {
0062     cout << "ERROR: The combination tool is not initialized\n";
0063     return;
0064   }
0065
0066   float energy = mySimpleElectron.getNewEnergy();
0067   float energyError = mySimpleElectron.getNewEnergyError();
0068   float momentum = mySimpleElectron.getTrackerMomentum();
0069   float momentumError = mySimpleElectron.getTrackerMomentumError();
0070   int electronClass = mySimpleElectron.getElClass();
0071   bool isEcalDriven = mySimpleElectron.isEcalDriven();
0072   bool isTrackerDriven = mySimpleElectron.isTrackerDriven();
0073   bool isEB = mySimpleElectron.isEB();
0074
0075   // compute relative errors and ratio of errors
0076   float energyRelError = energyError / energy;
0077   float momentumRelError = momentumError / momentum;
0078   float errorRatio = energyRelError / momentumRelError;
0079
0080   // calculate E/p and corresponding error
0081   float eOverP = energy / momentum;
0082   float eOverPerror =
0083       sqrt((energyError / momentum) * (energyError / momentum) +
0084            (energy * momentumError / momentum / momentum) * (energy * momentumError / momentum / momentum));
0085
0086   // fill input variables
0087   float regressionInputs[11];
0088   regressionInputs[0] = energy;
0089   regressionInputs[1] = energyRelError;
0090   regressionInputs[2] = momentum;
0091   regressionInputs[3] = momentumRelError;
0092   regressionInputs[4] = errorRatio;
0093   regressionInputs[5] = eOverP;
0094   regressionInputs[6] = eOverPerror;
0095   regressionInputs[7] = static_cast<float>(isEcalDriven);
0096   regressionInputs[8] = static_cast<float>(isTrackerDriven);
0097   regressionInputs[9] = static_cast<float>(electronClass);
0098   regressionInputs[10] = static_cast<float>(isEB);
0099
0100   // retrieve combination weight
0101   float weight = 0.;
0102   if (eOverP > 0.025 &&
0103       fabs(momentum - energy) < 15. * sqrt(momentumError * momentumError + energyError * energyError) &&
0104       ((momentumError < 10. * momentum) || (energy < 200.)))  // protection against crazy track measurement
0105   {
0106     weight = m_forest->GetResponse(regressionInputs);
0107     if (weight > 1.)
0108       weight = 1.;
0109     else if (weight < 0.)
0110       weight = 0.;
0111   }
0112
0113   float combinedMomentum = weight * momentum + (1. - weight) * energy;
0114   float combinedMomentumError =
0115       sqrt(weight * weight * momentumError * momentumError + (1. - weight) * (1. - weight) * energyError * energyError);
0116
0117   // FIXME : pure tracker electrons have track momentum error of 999.
0118   // If the combination try to combine such electrons then the original combined momentum is kept
0119   if (momentumError != 999. || weight == 0.) {
0120     mySimpleElectron.setCombinedMomentum(combinedMomentum);
0121     mySimpleElectron.setCombinedMomentumError(combinedMomentumError);
0122   }
0123 }