From 1868dbbf9555111a465b1b50bfe2909a4b6e5bd0 Mon Sep 17 00:00:00 2001
From: Christos Anastopoulos <Christos.Anastopoulos@cern.ch>
Date: Tue, 20 May 2014 16:51:21 +0200
Subject: [PATCH] Make dup removal const correct (egammaUtils-00-06-09)
---
.../egammaUtils/Root/EMConversionUtils.cxx | 119 ++++++++++++++++++
.../Root/ParameterDownWeighter.cxx | 31 +++++
.../egammaUtils/Root/WeightedMeanCalc.cxx | 70 +++++++++++
.../Root/egammaDuplicateRemoval.cxx | 45 +++++++
.../egamma/egammaUtils/cmt/requirements | 27 ++++
.../egamma/egammaUtils/doc/mainpage.h | 43 +++++++
.../egammaUtils/EMConversionUtils.h | 26 ++++
.../egammaUtils/EigenVectorAndMatrix.h | 33 +++++
.../egammaUtils/ParameterDownWeighter.h | 14 +++
.../egammaUtils/WeightedMeanCalc.h | 14 +++
.../egammaUtils/egammaDuplicateRemoval.h | 37 ++++++
11 files changed, 459 insertions(+)
create mode 100644 Reconstruction/egamma/egammaUtils/Root/EMConversionUtils.cxx
create mode 100644 Reconstruction/egamma/egammaUtils/Root/ParameterDownWeighter.cxx
create mode 100644 Reconstruction/egamma/egammaUtils/Root/WeightedMeanCalc.cxx
create mode 100644 Reconstruction/egamma/egammaUtils/Root/egammaDuplicateRemoval.cxx
create mode 100644 Reconstruction/egamma/egammaUtils/cmt/requirements
create mode 100644 Reconstruction/egamma/egammaUtils/doc/mainpage.h
create mode 100644 Reconstruction/egamma/egammaUtils/egammaUtils/EMConversionUtils.h
create mode 100644 Reconstruction/egamma/egammaUtils/egammaUtils/EigenVectorAndMatrix.h
create mode 100644 Reconstruction/egamma/egammaUtils/egammaUtils/ParameterDownWeighter.h
create mode 100644 Reconstruction/egamma/egammaUtils/egammaUtils/WeightedMeanCalc.h
create mode 100644 Reconstruction/egamma/egammaUtils/egammaUtils/egammaDuplicateRemoval.h
diff --git a/Reconstruction/egamma/egammaUtils/Root/EMConversionUtils.cxx b/Reconstruction/egamma/egammaUtils/Root/EMConversionUtils.cxx
new file mode 100644
index 00000000000..f8b5a44c49f
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/Root/EMConversionUtils.cxx
@@ -0,0 +1,119 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "egammaUtils/EMConversionUtils.h"
+#include "VxVertex/VxTrackAtVertex.h"
+#include "TrkParameters/TrackParameters.h"
+#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "TrkSurfaces/PerigeeSurface.h"
+
+namespace EMConversionUtils {
+
+ /** @brief struct of Particle Masses */
+ const Trk::ParticleMasses particleMasses;
+
+ /** @brief Get the position and momentum of the conversion candidate */
+ void getConversionPositionAndPt(const Trk::VxCandidate* myCandidate, Amg::Vector3D& pos, double& pt)
+ {
+
+ const std::vector<Trk::VxTrackAtVertex*> trkAtVx(*(myCandidate->vxTrackAtVertex()));
+ const Trk::RecVertex & vertex = myCandidate->recVertex();
+
+ pos = vertex.position();
+
+ //invariant mass, photon pT
+ if(trkAtVx.size() == 2) {
+ const Trk::Perigee* perigee1 = dynamic_cast<const Trk::Perigee*>(trkAtVx[0]->perigeeAtVertex());
+ const Trk::Perigee* perigee2 = dynamic_cast<const Trk::Perigee*>(trkAtVx[1]->perigeeAtVertex());
+ if (perigee1 && perigee2) {
+ const Amg::Vector3D sum_mom = perigee1->momentum() + perigee2->momentum();
+ pt = sum_mom.perp();
+ }
+ } else {
+ const Trk::Perigee* perigee1 = dynamic_cast<const Trk::Perigee*>(trkAtVx[0]->perigeeAtVertex());
+ if (perigee1) {
+ pt = perigee1->momentum().perp();
+ }
+ }
+ }
+
+ /** @brief Get the momentum of the conversion candidate */
+ double getConversionP(const Trk::VxCandidate* myCandidate)
+ {
+
+ const std::vector<Trk::VxTrackAtVertex*> trkAtVx(*(myCandidate->vxTrackAtVertex()));
+
+ //invariant mass, photon pT
+ if(trkAtVx.size() == 2) {
+ const Trk::Perigee* perigee1 = dynamic_cast<const Trk::Perigee*>(trkAtVx[0]->perigeeAtVertex());
+ const Trk::Perigee* perigee2 = dynamic_cast<const Trk::Perigee*>(trkAtVx[1]->perigeeAtVertex());
+ if (perigee1 && perigee2) {
+ const Amg::Vector3D sum_mom = perigee1->momentum() + perigee2->momentum();
+ return sum_mom.mag();
+ } else {
+ return -999.;
+ }
+ } else {
+ const Trk::Perigee* perigee1 = dynamic_cast<const Trk::Perigee*>(trkAtVx[0]->perigeeAtVertex());
+ if (perigee1) {
+ return perigee1->momentum().mag();
+ } else {
+ return -999.;
+ }
+ }
+ }
+
+ /** @brief Get the position of the conversion candidate */
+ Amg::Vector3D getConversionPosition(const Trk::VxCandidate* myCandidate)
+ {
+
+ const std::vector<Trk::VxTrackAtVertex*> trkAtVx(*(myCandidate->vxTrackAtVertex()));
+ const Trk::RecVertex & vertex = myCandidate->recVertex();
+
+ return vertex.position();
+ }
+
+ /** @brief Create neutral track parameters for the reconstructed converted photon. */
+ const Trk::NeutralPerigee* createPhotonParameters(const Trk::VxCandidate* myCandidate)
+ {
+ const std::vector<Trk::VxTrackAtVertex*> trkAtVx(*(myCandidate->vxTrackAtVertex()));
+
+ if(trkAtVx.size() == 2){
+ const Trk::Perigee* perigee1 = dynamic_cast<const Trk::Perigee*>(trkAtVx[0]->perigeeAtVertex());
+ const Trk::Perigee* perigee2 = dynamic_cast<const Trk::Perigee*>(trkAtVx[1]->perigeeAtVertex());
+ if (!perigee1 || !perigee2)
+ return 0;
+ const Amg::Vector3D sum_mom = perigee1->momentum() + perigee2->momentum();
+
+ //Get the photon inverse momentum
+ const double p = sum_mom.mag();
+ const double inv_p = (p != 0) ? 1.0/p : 0;
+
+ //const Trk::RecVertex & vertex = myCandidate->recVertex();
+
+ const AmgSymMatrix(5) cov1 = *(perigee1->covariance());
+ const AmgSymMatrix(5) cov2 = *(perigee2->covariance());
+ AmgSymMatrix(5)* eM = new AmgSymMatrix(5);
+ (*eM) = cov1+cov2;
+
+ //
+ const Trk::RecVertex & vertex = myCandidate->recVertex();
+ Amg::Vector3D pos(vertex.position().x(),vertex.position().y(),vertex.position().z());
+ const Trk::PerigeeSurface surface(pos);
+ //Construct the photon neutral parameters as neutral perigee
+ const Trk::NeutralPerigee* neutPar = const_cast<Trk::NeutralPerigee*>
+ (surface.createParameters<5,Trk::Neutral>(
+ 0.0 ,
+ 0.0,
+ sum_mom.phi(),
+ sum_mom.theta(),
+ inv_p,
+ eM));
+ return neutPar;
+ } else {
+ return 0;
+ }
+
+ }
+}
diff --git a/Reconstruction/egamma/egammaUtils/Root/ParameterDownWeighter.cxx b/Reconstruction/egamma/egammaUtils/Root/ParameterDownWeighter.cxx
new file mode 100644
index 00000000000..4b6590a6566
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/Root/ParameterDownWeighter.cxx
@@ -0,0 +1,31 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "egammaUtils/ParameterDownWeighter.h"
+#include <cfloat>
+
+namespace ParameterDownWeighter {
+
+ void getDownWeightedVM(EigenVectorAndMatrix& vm1, const EigenVectorAndMatrix& vm2, int indx) {
+
+ const double factor=10e10;
+
+ Amg::VectorX v1 = vm1.getVector();
+ Amg::MatrixX m1 = vm1.getMatrix();
+
+ Amg::VectorX v2 = vm2.getVector();
+
+ if (indx != 999) {
+ for (int i =0; i<5; i++) {
+ if ( i==indx ) {
+ v1[i] = v2[indx];
+ m1(i,indx)*=factor;
+ } else {
+ m1(i,indx)*=std::sqrt(factor);
+ m1(indx,i)*=std::sqrt(factor);
+ }
+ }
+ }
+ }
+}
diff --git a/Reconstruction/egamma/egammaUtils/Root/WeightedMeanCalc.cxx b/Reconstruction/egamma/egammaUtils/Root/WeightedMeanCalc.cxx
new file mode 100644
index 00000000000..338961a38f6
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/Root/WeightedMeanCalc.cxx
@@ -0,0 +1,70 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "egammaUtils/WeightedMeanCalc.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FourMomUtils/P4Helpers.h"
+#include <Eigen/LU>
+
+/* Note that the vector is of the form:
+
+ v[0] = d0
+ v[1] = z0
+ v[2] = phi
+ v[3] = eta
+ v[4] = E
+
+ The matrix is the associated covariance matrix.
+
+ Note that one has to be careful with phi wraparound when making the calculations
+ June 07, 2010 change to being a function in a namespace instead of a static
+*/
+
+namespace WeightedMeanCalc {
+ EigenVectorAndMatrix getWeightedMean( const EigenVectorAndMatrix& vm1, const EigenVectorAndMatrix& vm2)
+ {
+
+ // let's take care of the input for phi wraparound
+
+ Amg::VectorX v1 = vm1.getVector();
+ Amg::VectorX v2 = vm2.getVector();
+
+ const double diff = v1[2] - v2[2];
+
+ if (diff != P4Helpers::deltaPhi(v1[2], v2[2])) {
+ // temporarily make the smaller to be in [pi,3pi[
+ if (diff > 0) {
+ v2[2] += 2*M_PI;
+ } else {
+ v1[2] += 2*M_PI;
+ }
+ }
+ // end of input phi wraparound
+
+ Amg::MatrixX sumMatrix = vm1.getMatrix()+vm2.getMatrix();
+ Amg::MatrixX invSumMatrix = sumMatrix.inverse() ;
+ Amg::MatrixX weight = vm1.getMatrix()*invSumMatrix;
+
+ // get weighted mean vector and weighted mean matrix
+ Amg::VectorX combinedVec = vm1.getVector() + (weight*(vm2.getVector()-vm1.getVector()));
+ Amg::MatrixX combinedMat = weight*vm2.getMatrix();
+
+ // fix phi wraparound
+ if (combinedVec[2] >= M_PI) {
+ combinedVec[2] = fmod(combinedVec[2] + M_PI, 2*M_PI) - M_PI;
+ } else if (combinedVec[2] < - M_PI) {
+ combinedVec[2] = fmod(combinedVec[2] - M_PI, 2*M_PI) + M_PI;
+ // make M_PI be -M_PI
+ if (combinedVec[2] == M_PI) {
+ combinedVec[2] = - M_PI;
+ }
+ }
+
+ //Make sure combined matrix is symmetric
+ Amg::MatrixX combinedSymMatrix(combinedMat);
+
+ return EigenVectorAndMatrix(combinedVec, combinedSymMatrix);
+
+ }
+}
diff --git a/Reconstruction/egamma/egammaUtils/Root/egammaDuplicateRemoval.cxx b/Reconstruction/egamma/egammaUtils/Root/egammaDuplicateRemoval.cxx
new file mode 100644
index 00000000000..77dab3a6a50
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/Root/egammaDuplicateRemoval.cxx
@@ -0,0 +1,45 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "egammaUtils/egammaDuplicateRemoval.h"
+#include "FourMomUtils/P4Helpers.h"
+#include "AnalysisUtils/AnalysisMisc.h"
+
+namespace egammaDuplicateRemoval {
+
+ std::vector<const xAOD::CaloCluster *> getClusterDuplicateRemoval(std::vector<const xAOD::CaloCluster* > input_phList, double DETA /*=0.05 */, double DPHI /*=0.1*/) {
+
+ //sort w.r.t. cluster Et
+ AnalysisUtils::Sort::pT (&input_phList);
+
+ std::vector<bool> keepCand(input_phList.size(), true);
+ std::vector<const xAOD::CaloCluster*> output_phList;
+
+ // we need more than 2 clusters to check overlap
+ if ( input_phList.size() > 1 ) {
+ for (unsigned int iph1=0; iph1 < input_phList.size(); iph1++) {
+ if (keepCand[iph1]) {
+ for (unsigned int iph2=iph1+1; iph2 < input_phList.size(); iph2++) {
+ if (keepCand[iph2]) {
+ double dPhi = P4Helpers::deltaPhi( input_phList[iph1]->phi(), input_phList[iph2]->phi() );
+ double dEta = input_phList[iph1]->eta()- input_phList[iph2]->eta();
+ // find out if there are overlaps
+ if ( (fabs(dPhi) <= DPHI) && (fabs(dEta) <= DETA) ) {
+ keepCand[iph2]=false;
+ }
+ }
+ } // second loop on clusters
+ }
+ } // first loop on clusters
+ } // check input size
+
+ // fill the new vector with non-overlapped candidates
+ for (unsigned int i=0; i < input_phList.size(); i++) {
+ if ( keepCand[i]) output_phList.push_back(input_phList[i]);
+ }
+
+ return output_phList;
+}
+
+} // end of namespace
diff --git a/Reconstruction/egamma/egammaUtils/cmt/requirements b/Reconstruction/egamma/egammaUtils/cmt/requirements
new file mode 100644
index 00000000000..25813579526
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/cmt/requirements
@@ -0,0 +1,27 @@
+package egammaUtils
+
+use AtlasPolicy AtlasPolicy-*
+
+use AtlasROOT AtlasROOT-* External
+use xAODCaloEvent xAODCaloEvent-* Event/xAOD
+use VxVertex VxVertex-* Tracking/TrkEvent
+use TrkNeutralParameters TrkNeutralParameters-* Tracking/TrkEvent
+use EventPrimitives EventPrimitives-* Event
+
+private
+
+use AthenaBaseComps AthenaBaseComps-* Control
+use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
+use TrkParameters TrkParameters-* Tracking/TrkEvent
+use AtlasEigen AtlasEigen-* External
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+use FourMomUtils FourMomUtils-* Event
+use AnalysisUtils AnalysisUtils-* PhysicsAnalysis/AnalysisCommon
+
+apply_tag ROOTMathLibs
+
+end_private
+
+library egammaUtils ../Root/*.cxx
+apply_pattern installed_library
+
diff --git a/Reconstruction/egamma/egammaUtils/doc/mainpage.h b/Reconstruction/egamma/egammaUtils/doc/mainpage.h
new file mode 100644
index 00000000000..491a99496b5
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/doc/mainpage.h
@@ -0,0 +1,43 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+
+@mainpage egammaUtils Package
+
+This package provides various genaral tools for egamma use,
+but also called from other packages, such as
+
+BoostedDecisionTree
+Tool to calculate BDT - called e.g from egammaPIDTools package
+
+FisherWrapper
+TFermaAbsEngine
+TFermaFisher
+TFermaSample
+TFermaTemplate
+Tools to calculate Fisher discriminant - called e.g from egammaPIDTools package
+
+
+egammaMCTruthClassifier
+Starting from electron MC truth match or a track-match, the source of electron is defined.
+
+egammaBitMaskGenerator
+Tool to generate IsEM bitmask according to the user's requirements
+
+egammaTrackIsolationTool
+Tool to call to get track isolation: starting from CaloCluster
+Loops over TrackParticleContainer container calculating pT of the tracks which match certain
+selection criteria
+
+
+@author RD Schaffer <R.D.Schaffer@cern.ch>
+@author Srini Rajagopalan <srini@sun2.bnl.gov>
+@author Hong Ma <hma@bnl.gov>
+@author Peter Loch <loch@physics.arizona.edu>
+@author Frederic Derue <derue@lpnhep.in2p3.fr>
+@author Alexander Khodinov <khodinov@bnl.gov>
+
+*/
+
diff --git a/Reconstruction/egamma/egammaUtils/egammaUtils/EMConversionUtils.h b/Reconstruction/egamma/egammaUtils/egammaUtils/EMConversionUtils.h
new file mode 100644
index 00000000000..a7da84409e6
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/egammaUtils/EMConversionUtils.h
@@ -0,0 +1,26 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EGAMMAUTILS_EMCONVERSIONUTILS_H
+#define EGAMMAUTILS_EMCONVERSIONUTILS_H
+
+#include "VxVertex/ExtendedVxCandidate.h"
+#include "TrkNeutralParameters/NeutralParameters.h"
+
+namespace EMConversionUtils {
+
+ /** @brief Get the position and momentum of the conversion candidate */
+ void getConversionPositionAndPt(const Trk::VxCandidate*, Amg::Vector3D& p, double& pt);
+
+ /** @brief Get the momentum of the conversion candidate */
+ double getConversionP(const Trk::VxCandidate*);
+
+ /** @brief Get the position of the conversion candidate */
+ Amg::Vector3D getConversionPosition(const Trk::VxCandidate*);
+
+ /** @brief Create neutral track parameters for the reconstructed converted photon. */
+ const Trk::NeutralPerigee* createPhotonParameters(const Trk::VxCandidate * myCandidate);
+}
+
+#endif
diff --git a/Reconstruction/egamma/egammaUtils/egammaUtils/EigenVectorAndMatrix.h b/Reconstruction/egamma/egammaUtils/egammaUtils/EigenVectorAndMatrix.h
new file mode 100644
index 00000000000..e18b898c8b3
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/egammaUtils/EigenVectorAndMatrix.h
@@ -0,0 +1,33 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EGAMMAUTILS_HEPVECTORANDMATRIX_H
+#define EGAMMAUTILS_HEPVECTORANDMATRIX_H
+
+#include "EventPrimitives/EventPrimitives.h"
+
+class EigenVectorAndMatrix {
+
+public:
+ EigenVectorAndMatrix(Amg::VectorX &v, Amg::MatrixX &m) : vector(v), matrix(m) {
+ }
+ Amg::VectorX& getVector() {
+ return vector;
+ }
+ Amg::MatrixX& getMatrix() {
+ return matrix;
+ }
+ const Amg::VectorX& getVector() const {
+ return vector;
+ }
+ const Amg::MatrixX& getMatrix() const {
+ return matrix;
+ }
+
+private:
+ Amg::VectorX vector;
+ Amg::MatrixX matrix;
+};
+
+#endif
diff --git a/Reconstruction/egamma/egammaUtils/egammaUtils/ParameterDownWeighter.h b/Reconstruction/egamma/egammaUtils/egammaUtils/ParameterDownWeighter.h
new file mode 100644
index 00000000000..bbf15363a43
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/egammaUtils/ParameterDownWeighter.h
@@ -0,0 +1,14 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EGAMMAUTILS_PARAMETERDOWNWEIGHTER_H
+#define EGAMMAUTILS_PARAMETERDOWNWEIGHTER_H
+
+#include "egammaUtils/EigenVectorAndMatrix.h"
+
+namespace ParameterDownWeighter {
+ void getDownWeightedVM(EigenVectorAndMatrix& vm1, const EigenVectorAndMatrix& vm2, int indx);
+}
+
+#endif
diff --git a/Reconstruction/egamma/egammaUtils/egammaUtils/WeightedMeanCalc.h b/Reconstruction/egamma/egammaUtils/egammaUtils/WeightedMeanCalc.h
new file mode 100644
index 00000000000..eb86b36afbb
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/egammaUtils/WeightedMeanCalc.h
@@ -0,0 +1,14 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EGAMMAUTILS_WEIGHTEDMEANCALC_H
+#define EGAMMAUTILS_WEIGHTEDMEANCALC_H
+
+#include "egammaUtils/EigenVectorAndMatrix.h"
+
+namespace WeightedMeanCalc {
+ EigenVectorAndMatrix getWeightedMean(const EigenVectorAndMatrix& vm1, const EigenVectorAndMatrix& vm2);
+}
+
+#endif
diff --git a/Reconstruction/egamma/egammaUtils/egammaUtils/egammaDuplicateRemoval.h b/Reconstruction/egamma/egammaUtils/egammaUtils/egammaDuplicateRemoval.h
new file mode 100644
index 00000000000..0be4a4b8804
--- /dev/null
+++ b/Reconstruction/egamma/egammaUtils/egammaUtils/egammaDuplicateRemoval.h
@@ -0,0 +1,37 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef EGAMMAUTILS_EGAMMADUPLICATEREMOVAL_H
+#define EGAMMAUTILS_EGAMMADUPLICATEREMOVAL_H
+
+
+/// @brief egammaDuplicateRemoval
+/// @brief Remove duplicate using cluster
+
+/// @author T. Cuhadar Donszelmann
+/// @author Christos Anastopoulos
+///
+/// $Revision:$
+/// $Date: 2014-02-11 17:40:48 +0100 (Tue, 11 Feb 2014) $
+///
+
+/********************************************************************
+NAME:
+PACKAGE: offline/Reconstruction/egammaUtils
+AUTHORS:
+CREATED: Nov 2011
+PURPOSE: removal
+UPDATED:
+********************************************************************/
+
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "xAODCaloEvent/CaloCluster.h"
+#include <vector>
+
+
+namespace egammaDuplicateRemoval {
+ std::vector<const xAOD::CaloCluster *> getClusterDuplicateRemoval(std::vector<const xAOD::CaloCluster*> input_Ptr, double DETA=0.05, double DPHI=0.1);
+}
+
+#endif
--
GitLab