Commit 7ac93034 authored by Christos Anastopoulos's avatar Christos Anastopoulos Committed by Graeme Stewart
Browse files

Coverity fixes (egammaMVACalib-01-00-16)

parent a269c5ca
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include "egammaMVACalib/BDT.h"
#include "egammaMVACalib/Node.h"
#include "TMVA/MethodBDT.h"
#include "TMVA/DecisionTree.h"
#include "TTree.h"
using namespace egammaMVACalibNmsp;
/** c-tor from TTree **/
BDT::BDT(TTree *tree)
{
std::vector<int> *vars = 0;
std::vector<float> *values = 0;
tree->SetBranchAddress("offset", &m_offset);
tree->SetBranchAddress("vars", &vars);
tree->SetBranchAddress("values", &values);
for (int i = 0; i < tree->GetEntries(); ++i)
{
tree->GetEntry(i);
assert (vars);
assert (values);
int j = 0;
m_forest.push_back( Node::newNode(*vars, *values, j) );
}
}
/** c-tor from TMVA::MethodBDT **/
BDT::BDT(TMVA::MethodBDT* bdt)
{
assert(bdt);
m_offset = bdt->GetBoostWeights().size() ? bdt->GetBoostWeights()[0] : 0.;
std::vector<TMVA::DecisionTree*>::const_iterator it;
for(it = bdt->GetForest().begin(); it != bdt->GetForest().end(); ++it)
m_forest.push_back( new Node( (*it)->GetRoot()) );
}
BDT::~BDT()
{
std::vector<Node*>::iterator it;
for(it = m_forest.begin(); it != m_forest.end(); ++it)
delete *it;
}
/** Return offset + the sum of the response of each tree **/
float BDT::GetResponse(const std::vector<float>& values) const
{
float result = m_offset;
std::vector<Node*>::const_iterator it;
for (it = m_forest.begin(); it != m_forest.end(); ++it)
result += (*it)->GetResponse(values);
return result;
}
/** Return offset + the sum of the response of each tree **/
float BDT::GetResponse(const std::vector<float*>& pointers) const
{
float result = m_offset;
std::vector<Node*>::const_iterator it;
for (it = m_forest.begin(); it != m_forest.end(); ++it)
result += (*it)->GetResponse(pointers);
return result;
}
/** Return the values corresponding to m_pointers (or an empty vector) **/
std::vector<float> BDT::GetValues() const
{
std::vector<float> result;
std::vector<float*>::const_iterator it;
for (it = m_pointers.begin(); it != m_pointers.end(); ++it)
{
float *ptr = *it;
assert (ptr);
assert (*ptr);
result.push_back(*ptr);
}
return result;
}
/** Return a TTree representing the BDT:
* each entry is a binary tree, each element of the vectors is a node
**/
TTree* BDT::WriteTree(TString name)
{
std::vector<int> vars;
std::vector<float> values;
TTree *tree = new TTree(name.Data(), "BDT tree");
tree->Branch("offset", &m_offset);
tree->Branch("vars", &vars);
tree->Branch("values", &values);
std::vector<Node*>::iterator it;
for(it = m_forest.begin(); it != m_forest.end(); ++it)
{
vars.clear();
values.clear();
(*it)->AddToVectors(vars, values);
tree->Fill();
}
return tree;
}
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifdef __CINT__
#include "egammaMVACalib/egammaMVACalib.h"
#include "egammaMVACalib/Node.h"
#include "egammaMVACalib/BDT.h"
#include "egammaMVACalib/egammaMVATool.h"
#include <vector>
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclass;
#pragma link C++ nestedtypedef;
#pragma link C++ class egammaMVACalib+;
#pragma link C++ class egammaMVACalib::ReaderID+;
#pragma link C++ class egammaMVACalibNmsp::Node+;
#pragma link C++ class egammaMVACalibNmsp::BDT+;
#pragma link C++ class std::vector<egammaMVACalibNmsp::Node*>+;
#pragma link C++ class egammaMVATool;
#pragma link C++ enum egammaMVACalib::ParticleType+;
#pragma link C++ enum egammaMVACalib::CalibrationType+;
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include "egammaMVACalib/Node.h"
#include <cassert>
using namespace std; // has to come before TMVA include (errors on header)
#include "TMVA/DecisionTreeNode.h"
using namespace egammaMVACalibNmsp;
/** c-tor from TMVA::DecisionTree node.
* Creates the full tree structure when called from the top node
**/
Node::Node(TMVA::DecisionTreeNode *node) : m_left(0), m_right(0)
{
if (!node->GetLeft())
{
m_var = -1;
m_val = node->GetResponse();
}
else
{
m_var = node->GetSelector();
m_val = node->GetCutValue();
m_left = new Node(node->GetCutType() ? node->GetRight() : node->GetLeft() );
m_right = new Node(node->GetCutType() ? node->GetLeft() : node->GetRight() );
}
}
/** Recursively deletes the node and its children **/
Node::~Node()
{
delete m_left;
delete m_right;
m_left = 0;
m_right = 0;
}
/** Recursive function to create nodes and their children.
Creates the full tree structure when called with i=0 **/
Node* Node::newNode(std::vector<int>& vars, std::vector<float>& values, int& i)
{
assert ( (int) vars.size() >= i);
Node *node = new Node(vars.at(i), values.at(i));
if (node->GetVar() >= 0)
{
node->m_left = newNode(vars, values, ++i);
node->m_right = newNode(vars, values, ++i);
}
return node;
}
/** Recursive function to evaluate the response of the decision tree given a vector with
* values for each variable
**/
float Node::GetResponse(const std::vector<float>& values) const
{
if (m_var < 0) return m_val;
assert ( (int) values.size() >= m_var);
Node *node = GetNext(values[m_var]);
assert (node);
return node->GetResponse(values);
}
/** Recursive function to evaluate the response of the decision tree given a vector with
* pointers to each variable
**/
float Node::GetResponse(const std::vector<float*>& pointers) const
{
if (m_var < 0) return m_val;
assert ( (int) pointers.size() >= m_var && pointers[m_var]);
Node *node = GetNext( *(pointers[m_var]) );
assert (node);
return node->GetResponse(pointers);
}
/** Recursive function to add the node info to vectors **/
void Node::AddToVectors(std::vector<int>& vars, std::vector<float>& values) const
{
vars.push_back(m_var);
values.push_back(m_val);
if (m_left) m_left->AddToVectors(vars, values);
if (m_right) m_right->AddToVectors(vars, values);
}
This diff is collapsed.
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
//#include "GaudiKernel/MsgStream.h"
// C++
// Stuff be here...
#include "egammaMVACalib/egammaMVATool.h"
#include "egammaMVACalib/egammaMVACalib.h"
// xAOD
#include "xAODEgamma/Egamma.h"
#include "xAODEgamma/EgammaDefs.h"
#include "xAODEgamma/Electron.h"
#include "xAODEgamma/Photon.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "PATInterfaces/CorrectionCode.h"
#include "xAODEgamma/EgammaxAODHelpers.h"
// tracking
#include "xAODTracking/Vertex.h"
#include "xAODTracking/TrackParticle.h"
#include "xAODTracking/TrackingPrimitives.h"
egammaMVATool::egammaMVATool( const std::string &name )
: asg::AsgTool(name)
{
// Configurable properties...
declareProperty("folder", m_folder="");//, "folder with weight files");
}
StatusCode egammaMVATool::initialize(){
ATH_MSG_DEBUG("In initialize of " << name() << "..." );
m_mvaElectron = new egammaMVACalib(egammaMVACalib::egELECTRON, // particle type
true, // use new BDT (not TMVA)
m_folder, // folder with weight files
"BDTG", // method
1, // Full Calib
false, // not debugging
"", // that
"",
"",
"", // file pattern of xml files
true // ignore spectators
);
m_mvaElectron->InitTree(0);
m_mvaPhoton = new egammaMVACalib(egammaMVACalib::egPHOTON, // particle type
true, // use new BDT (not TMVA)
m_folder, // folder with weight files
"BDTG", // method
1 , // Full Calib
false, // not debugging
"", // that
"",
"",
"", // file pattern of xml files
true // ignore spectators
);
m_mvaPhoton->InitTree(0);
return StatusCode::SUCCESS;
}
StatusCode egammaMVATool::finalize(){
ATH_MSG_DEBUG( "in finalize" );
return StatusCode::SUCCESS;
}
StatusCode egammaMVATool::execute(xAOD::CaloCluster* cluster,const xAOD::Egamma* eg){
ATH_MSG_DEBUG( "in finalize" );
if(!eg || !cluster){
ATH_MSG_ERROR("Invalid Pointer to egamma or cluster object");
return StatusCode::FAILURE;
}
double mvaE = getEnergy(cluster, eg);
ATH_MSG_DEBUG( "Calculated MVA calibrated energy = " << mvaE );
if(mvaE>0){
cluster->setCalE(mvaE);
}
else{
cluster->setCalE(cluster->e());
}
return StatusCode::SUCCESS;
}
StatusCode egammaMVATool::hltexecute(xAOD::CaloCluster* cluster,std::string egType) {
ATH_MSG_DEBUG( "in finalize" );
if(!cluster){
ATH_MSG_ERROR("Invalid Pointer to cluster object");
return StatusCode::FAILURE;
}
double mvaE = getEnergy(cluster, egType);
ATH_MSG_DEBUG( "Calculated MVA calibrated energy = " << mvaE );
if(mvaE>0){
cluster->setCalE(mvaE);
}
else{
cluster->setCalE(cluster->e());
}
return StatusCode::SUCCESS;
}
float egammaMVATool::getEnergy(xAOD::CaloCluster* cluster, std::string egType){
ATH_MSG_DEBUG("In execute...");
// Check for errors...
if ( !cluster ){
ATH_MSG_WARNING("no xAOD::CaloCluster object provided");
return 0;
}
getClusterVariables(cluster);
if(egType == "Electron"){
ATH_MSG_DEBUG("Processing for electron...");
m_mvaCalibratedEnergy = m_mvaElectron->getMVAEnergyElectron(m_rawcl_Es0,
m_rawcl_Es1,
m_rawcl_Es2,
m_rawcl_Es3,
m_cl_eta,
m_cl_E,
m_cl_etaCalo,
m_cl_phiCalo);
return m_mvaCalibratedEnergy;
}
else if(egType == "Photon"){
ATH_MSG_DEBUG("Processing for photon...");
m_ptconv = 0;
m_pt1conv = 0;
m_pt2conv = 0;
m_convtrk1nPixHits = 0;
m_convtrk1nSCTHits = 0;
m_convtrk2nPixHits = 0;
m_convtrk2nSCTHits = 0;
m_Rconv = 0;
m_mvaCalibratedEnergy = m_mvaPhoton->getMVAEnergyPhoton(m_rawcl_Es0,
m_rawcl_Es1,
m_rawcl_Es2,
m_rawcl_Es3,
m_cl_eta,
m_cl_E,
m_cl_etaCalo,
m_cl_phiCalo,
m_ptconv,
m_pt1conv,
m_pt2conv,
m_convtrk1nPixHits,
m_convtrk1nSCTHits,
m_convtrk2nPixHits,
m_convtrk2nSCTHits,
m_Rconv
);
return m_mvaCalibratedEnergy;
}
else{
ATH_MSG_INFO("Unknown Type");
}
return 0;
}
float egammaMVATool::getEnergy(xAOD::CaloCluster* cluster, const xAOD::Egamma* eg){
ATH_MSG_DEBUG("In execute...");
// Check for errors...
if ( !eg ){
ATH_MSG_WARNING("no xAOD::Egamma object provided");
return 0;
}
if( eg->type() == xAOD::Type::Electron ){
ATH_MSG_DEBUG("Processing for electron...");
return getEnergy(cluster, static_cast<const xAOD::Electron*>(eg));
}
else if (eg->type() == xAOD::Type::Photon ){
ATH_MSG_DEBUG("Processing for photon...");
return getEnergy(cluster, static_cast<const xAOD::Photon*>(eg));
}
else{
ATH_MSG_INFO("Unknown Type");
}
return 0;
}
float egammaMVATool::getEnergy(xAOD::CaloCluster* cluster, const xAOD::Electron* el){
if(!el){
ATH_MSG_ERROR("No electron passed");
return 0;
}
getClusterVariables(cluster);
m_mvaCalibratedEnergy = m_mvaElectron->getMVAEnergyElectron(m_rawcl_Es0,
m_rawcl_Es1,
m_rawcl_Es2,
m_rawcl_Es3,
m_cl_eta,
m_cl_E,
m_cl_etaCalo,
m_cl_phiCalo);
return m_mvaCalibratedEnergy;
}
float egammaMVATool::getEnergy(xAOD::CaloCluster* cluster, const xAOD::Photon* ph){
if(!ph){
ATH_MSG_ERROR("No photon passed");
return 0;
}
if( !getClusterVariables(cluster)){
ATH_MSG_ERROR("Could not retrieve cluster variables.");
return 0;
}
const xAOD::Vertex* phVertex = ph->vertex();
if(!phVertex){
m_ptconv = 0;
m_pt1conv = 0;
m_pt2conv = 0;
m_convtrk1nPixHits = 0;
m_convtrk1nSCTHits = 0;
m_convtrk2nPixHits = 0;
m_convtrk2nSCTHits = 0;
m_Rconv = 0;
} else{
getConversionVariables(phVertex);
}
m_mvaCalibratedEnergy = m_mvaPhoton->getMVAEnergyPhoton(m_rawcl_Es0,
m_rawcl_Es1,
m_rawcl_Es2,
m_rawcl_Es3,
m_cl_eta,
m_cl_E,
m_cl_etaCalo,
m_cl_phiCalo,
m_ptconv,
m_pt1conv,
m_pt2conv,
m_convtrk1nPixHits,
m_convtrk1nSCTHits,
m_convtrk2nPixHits,
m_convtrk2nSCTHits,
m_Rconv
);
return m_mvaCalibratedEnergy;
}
bool egammaMVATool::getClusterVariables(const xAOD::CaloCluster* cluster){
// The cluster energy in the sampling layers
m_rawcl_Es0 = cluster->energyBE(0);
m_rawcl_Es1 = cluster->energyBE(1);
m_rawcl_Es2 = cluster->energyBE(2);
m_rawcl_Es3 = cluster->energyBE(3);
m_cl_eta = cluster->eta();//el_cl_eta;
m_cl_E = cluster->e(); //el_cl_E;
// Correct this with the right variables
m_cl_etaCalo = 0.;
if(! (cluster->retrieveMoment(xAOD::CaloCluster::ETACALOFRAME,m_cl_etaCalo) ) ){
ATH_MSG_ERROR("etaCalo not found in CaloCluster object");
return false;
}
m_cl_phiCalo = 0.;
if(! (cluster->retrieveMoment(xAOD::CaloCluster::PHICALOFRAME,m_cl_phiCalo) ) ){
ATH_MSG_ERROR("phiCalo not found in CaloCluster object");
return false;
}
return true;
}
bool egammaMVATool::getConversionVariables(const xAOD::Vertex *phVertex){
const Amg::Vector3D pos = phVertex->position();
m_Rconv = static_cast<float>(hypot (pos.x(), pos.y()));
const xAOD::TrackParticle* tp0 = phVertex->trackParticle(0);
const xAOD::TrackParticle* tp1 = phVertex->trackParticle(1);
TLorentzVector sum;
if(tp0){
sum += tp0->p4();
uint8_t hits;
tp0->summaryValue(hits,xAOD::numberOfPixelHits);
m_convtrk1nPixHits = hits;
tp0->summaryValue(hits,xAOD::numberOfSCTHits);
m_convtrk1nSCTHits = hits;
m_pt1conv = static_cast<float>(tp0->pt());
}
if(tp1){
sum += tp1->p4();
uint8_t hits;
tp1->summaryValue(hits,xAOD::numberOfPixelHits);
m_convtrk2nPixHits = hits;
tp1->summaryValue(hits,xAOD::numberOfSCTHits);
m_convtrk2nSCTHits = hits;
m_pt2conv = static_cast<float>(tp1->pt());
}
m_ptconv = sum.Perp();
return true;
}
# this makefile also gets parsed by shell scripts
# therefore it does not support full make syntax and features
# edit with care
PACKAGE = egammaMVACalib
PACKAGE_PRELOAD = Hist TMVA PATInterfaces
PACKAGE_CXXFLAGS =
PACKAGE_LDFLAGS =
PACKAGE_BINFLAGS =
PACKAGE_DEP = AsgTools xAODEgamma xAODCaloEvent xAODEventInfo xAODRootAccess EventLoop PathResolver
PACKAGE_NOOPT = 0
PACKAGE_PEDANTIC = 1
include $(ROOTCOREDIR)/Makefile-common
#################################################
package egammaMVACalib
author Christos Anastopoulos
author Ruggero Turra <ruggero.turra@cern.ch>
author Bruno Lenzi <Bruno.Lenzi@cern.ch>
author Javier Salazar <jsalazar@cern.ch>
use AtlasPolicy AtlasPolicy-*
use AtlasROOT AtlasROOT-* External
use GaudiInterface GaudiInterface-* External
############################################################
# For the interface for the athena tool
use AsgTools AsgTools-* Control/AthToolSupport
use xAODEgamma xAODEgamma-* Event/xAOD
use xAODCaloEvent xAODCaloEvent-* Event/xAOD
use xAODTracking xAODTracking-* Event/xAOD
use PATInterfaces PATInterfaces-* PhysicsAnalysis/AnalysisCommon
# Add cmake compatibility (doesn't do anything on CMT side of things)
apply_pattern cmake_add_command command="find_package(ROOT COMPONENTS Tree TreePlayer TMVA XMLIO)"
private
use PathResolver PathResolver-* Tools
use AtlasReflex AtlasReflex-* External