Replacement of TMatrix and TVector by SMatrix and SVector in MET code

DataFormats/METReco/interface/MET.h

@@ -31,11 +31,14 @@
 #include <cmath>
 #include <vector>
 #include <cstring>
-#include "TMatrixD.h"
+#include <Math/SMatrix.h>
+#include <Math/SVector.h>
 
 //____________________________________________________________________________||
 namespace reco
 {
+  typedef ROOT::Math::SMatrix<double,2> METCovMatrix;
+
   class MET : public RecoCandidate
   {
   public:
@@ -67,8 +70,8 @@ namespace reco
     std::vector<CorrMETData> mEtCorr() const { return corr; }
 
     //________________________________________________________________________||
-    void setSignificanceMatrix(const TMatrixD& matrix);
-    TMatrixD getSignificanceMatrix(void) const;
+    void setSignificanceMatrix(const reco::METCovMatrix& matrix);
+    reco::METCovMatrix getSignificanceMatrix(void) const;
 
   private:
     virtual bool overlap( const Candidate & ) const;

DataFormats/METReco/interface/SigInputObj.h

@@ -27,6 +27,7 @@
 
 //=== Class SigInputObj ==============================//
 namespace metsig{
+
   class SigInputObj{
 
   public:

DataFormats/METReco/src/MET.cc

@@ -78,14 +78,16 @@ MET * MET::clone() const {
 double MET::significance() const {
   if(signif_dxx==0 && signif_dyy==0 && signif_dxy==0 && signif_dyx==0)
     return -1;
-  TMatrixD metmat = getSignificanceMatrix();
-  TVectorD metvec(2);
+  METCovMatrix metmat = getSignificanceMatrix();
+  ROOT::Math::SVector<double,2> metvec;
   metvec(0)=this->px();
   metvec(1)=this->py();
   double signif = -1;
-  if(std::fabs(metmat.Determinant())>0.000001){
+  double det=0;
+  metmat.Det(det);
+  if(std::fabs(det)>0.000001){
     metmat.Invert();
-    signif = metvec * (metmat * metvec);
+    signif = ROOT::Math::Dot(metvec, (metmat * metvec) );
   }
   return signif;
 }
@@ -134,9 +136,9 @@ std::vector<double> MET::dsumEt() const
 
 // returns the significance matrix
 //____________________________________________________________________________||
-TMatrixD MET::getSignificanceMatrix(void) const
+METCovMatrix MET::getSignificanceMatrix(void) const
 {
-  TMatrixD result(2,2);
+  METCovMatrix result;
   result(0,0)=signif_dxx;
   result(0,1)=signif_dxy;
   result(1,0)=signif_dyx;
@@ -152,7 +154,7 @@ bool MET::overlap( const Candidate & ) const
 }
 
 //____________________________________________________________________________||
-void MET::setSignificanceMatrix(const TMatrixD &inmatrix)
+void MET::setSignificanceMatrix(const METCovMatrix &inmatrix)
 {
   signif_dxx=inmatrix(0,0);
   signif_dxy=inmatrix(0,1);

RecoMET/METAlgorithms/interface/SignCaloSpecificAlgo.h

@@ -19,12 +19,12 @@
 #define METProducers_SignCaloMETAlgo_h
 
 //____________________________________________________________________________||
+#include "DataFormats/METReco/interface/MET.h"
 #include "DataFormats/METReco/interface/CommonMETData.h"
 #include "DataFormats/Common/interface/Handle.h"
 #include "DataFormats/Common/interface/View.h"
 #include "DataFormats/Candidate/interface/Candidate.h"
 #include "DataFormats/METReco/interface/SigInputObj.h"
-#include "TMatrixD.h"
 
 namespace metsig {
   class SignAlgoResolutions;
@@ -39,9 +39,9 @@ class SignCaloSpecificAlgo {
   ~SignCaloSpecificAlgo();
 
   void usePreviousSignif(const std::vector<double> &values);
-  void usePreviousSignif(const TMatrixD &matrix) { matrix_ = matrix; }
+  void usePreviousSignif(const reco::METCovMatrix &matrix) { matrix_ = matrix; }
   double getSignificance(){return significance_;}
-  TMatrixD getSignificanceMatrix()const {return matrix_;}
+  reco::METCovMatrix getSignificanceMatrix()const {return matrix_;}
 
   void calculateBaseCaloMET(edm::Handle<edm::View<reco::Candidate> > towers,  const CommonMETData& met, const metsig::SignAlgoResolutions & resolutions, bool noHF, double globalthreshold);
 
@@ -50,7 +50,7 @@ class SignCaloSpecificAlgo {
   std::vector<metsig::SigInputObj> makeVectorOutOfCaloTowers(edm::Handle<edm::View<reco::Candidate> > towers, const metsig::SignAlgoResolutions& resolutions, bool noHF, double globalthreshold);
 
   double significance_;
-  TMatrixD matrix_;
+  reco::METCovMatrix matrix_;
 };
 
 

RecoMET/METAlgorithms/interface/SignPFSpecificAlgo.h

@@ -21,8 +21,8 @@
 #include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
 #include "RecoMET/METAlgorithms/interface/SignAlgoResolutions.h"
 #include "DataFormats/JetReco/interface/PFJet.h"
+#include "DataFormats/METReco/interface/MET.h"
 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
-#include "TMatrixD.h"
 
 //____________________________________________________________________________||
 namespace metsig
@@ -39,8 +39,8 @@ namespace metsig
     void addPFJets(const edm::View<reco::PFJet>* PFJets);
     void addPFCandidate(reco::PFCandidatePtr pf);
     void useOriginalPtrs(const edm::ProductID& productID);
-    TMatrixD getSignifMatrix() const {return algo_.getSignifMatrix();}
-    TMatrixD mkSignifMatrix(edm::Handle<edm::View<reco::Candidate> > &PFCandidates);
+    reco::METCovMatrix getSignifMatrix() const {return algo_.getSignifMatrix();}
+    reco::METCovMatrix mkSignifMatrix(edm::Handle<edm::View<reco::Candidate> > &PFCandidates);
 
   private:
     metsig::SignAlgoResolutions *resolutions_;

RecoMET/METAlgorithms/interface/significanceAlgo.h

@@ -77,6 +77,7 @@
 #include <sstream>
 
 #include "DataFormats/METReco/interface/SigInputObj.h"
+#include "DataFormats/METReco/interface/MET.h"
 #include "TMatrixTBase.h"
 #include "TMatrixD.h"
 #include "TVectorD.h"
@@ -87,18 +88,18 @@ namespace metsig{
     significanceAlgo();
     ~significanceAlgo();
 
-    const void addSignifMatrix(const TMatrixD &input);
-    const void setSignifMatrix(const TMatrixD &input,const double &met_r, const double &met_phi, const double &met_set);
+    const void addSignifMatrix(const reco::METCovMatrix &input);
+    const void setSignifMatrix(const reco::METCovMatrix &input,const double &met_r, const double &met_phi, const double &met_set);
     const double significance(double& met_r, double& met_phi, double& met_set);
     const void addObjects(const std::vector<metsig::SigInputObj>& EventVec);
     const void subtractObjects(const std::vector<metsig::SigInputObj>& EventVec);
-    TMatrixD getSignifMatrix() const {return signifmatrix_;}
+    reco::METCovMatrix getSignifMatrix() const {return signifmatrix_;}
     //    const std::vector<metsig::SigInputObj> eventVec(){return eventVec_;}
   private:
-    void rotateMatrix( Double_t theta, TMatrixD &v);  
+    void rotateMatrix( Double_t theta, reco::METCovMatrix  &v);  
 
     //    std::vector<metsig::SigInputObj> eventVec_;
-    TMatrixD signifmatrix_;
+    reco::METCovMatrix signifmatrix_;
     // workers:
     double set_worker_;
     double xmet_;

RecoMET/METAlgorithms/src/SignCaloSpecificAlgo.cc

@@ -23,8 +23,7 @@ using namespace std;
 
 //____________________________________________________________________________||
 SignCaloSpecificAlgo::SignCaloSpecificAlgo():
-  significance_(0.),
-  matrix_(2,2)
+  significance_(0.)
 {
   matrix_(0,0)=matrix_(1,0)=matrix_(0,1)=matrix_(1,1)=0.;
 }

RecoMET/METAlgorithms/src/SignPFSpecificAlgo.cc

@@ -78,7 +78,7 @@ void metsig::SignPFSpecificAlgo::addPFCandidate(reco::PFCandidatePtr pf)
 }
 
 //____________________________________________________________________________||
-TMatrixD metsig::SignPFSpecificAlgo::mkSignifMatrix(edm::Handle<edm::View<reco::Candidate> > &PFCandidates)
+reco::METCovMatrix metsig::SignPFSpecificAlgo::mkSignifMatrix(edm::Handle<edm::View<reco::Candidate> > &PFCandidates)
 {
   useOriginalPtrs(PFCandidates.id());
   for(edm::View<reco::Candidate>::const_iterator iParticle = (PFCandidates.product())->begin(); iParticle != (PFCandidates.product())->end(); ++iParticle )

RecoMET/METAlgorithms/src/significanceAlgo.cc

@@ -26,7 +26,6 @@
 
 metsig::significanceAlgo::significanceAlgo():
   //  eventVec_(0),
-  signifmatrix_(2,2),
   set_worker_(0),
   xmet_(0),
   ymet_(0)
@@ -39,25 +38,20 @@ metsig::significanceAlgo::significanceAlgo():
 
 //******* Add an existing significance matrix to the algo, so that the vector sum can be continued. Only makes sense if matrix is empty or you want to purposefully increase uncertainties (for example in systematic studies)!
 const void
-metsig::significanceAlgo::addSignifMatrix(const TMatrixD &input){
-  // check that the matrix is size 2:
-  if(input.GetNrows()==2 && input.GetNcols()==2) {
-    signifmatrix_+=input;
-  }
+metsig::significanceAlgo::addSignifMatrix(const reco::METCovMatrix& input){
+  signifmatrix_+=input;
   return;
 }
 ////////////////////////
 /// reset the signficance matrix (this is the most likely case), so that the vector sum can be continued
 
 const void
-metsig::significanceAlgo::setSignifMatrix(const TMatrixD &input,const double &met_r, const double &met_phi, const double &met_set){
-  // check that the matrix is size 2:
-  if(input.GetNrows()==2 && input.GetNcols()==2) {
-    signifmatrix_=input;
-    set_worker_=met_set;
-    xmet_=met_r*cos(met_phi);
-    ymet_=met_r*sin(met_phi);
-  }
+metsig::significanceAlgo::setSignifMatrix(const reco::METCovMatrix &input,const double &met_r, const double &met_phi, const double &met_set){
+  signifmatrix_=input;
+  set_worker_=met_set;
+  xmet_=met_r*cos(met_phi);
+  ymet_=met_r*sin(met_phi);
+
   return;
 }
 
@@ -69,11 +63,11 @@ metsig::significanceAlgo::~significanceAlgo(){
 // //*** rotate a 2D matrix by angle theta **********************//
 
 void
-metsig::significanceAlgo::rotateMatrix( Double_t theta, TMatrixD &v)
+metsig::significanceAlgo::rotateMatrix( Double_t theta, reco::METCovMatrix &v)
 {
   // I suggest not using this to rotate trivial matrices.
-  TMatrixD r(2,2);
-  TMatrixD rInv(2,2);
+  reco::METCovMatrix r;
+  reco::METCovMatrix rInv;
 
   r(0,0) = cos(theta); r(0,1) = sin(theta); r(1,0) = -sin(theta); r(1,1) = cos(theta);
   rInv = r;
@@ -87,7 +81,7 @@ metsig::significanceAlgo::rotateMatrix( Double_t theta, TMatrixD &v)
 const void 
 metsig::significanceAlgo::subtractObjects(const std::vector<metsig::SigInputObj>& eventVec)
 { 
-  TMatrixD v_tot = signifmatrix_;
+  reco::METCovMatrix v_tot = signifmatrix_;
   //--- Loop over physics objects in the event ---//
   //  for(unsigned int objnum=1; objnum < EventVec.size(); objnum++ ) {
   for(std::vector<SigInputObj>::const_iterator obj = eventVec.begin(); obj!= eventVec.end(); ++obj){
@@ -121,7 +115,7 @@ metsig::significanceAlgo::subtractObjects(const std::vector<metsig::SigInputObj>
 const void 
 metsig::significanceAlgo::addObjects(const std::vector<metsig::SigInputObj>& eventVec)
 { 
-  TMatrixD v_tot = signifmatrix_;
+  reco::METCovMatrix v_tot = signifmatrix_;
   //--- Loop over physics objects in the event ---//
   //  for(unsigned int objnum=1; objnum < EventVec.size(); objnum++ ) {
   for(std::vector<SigInputObj>::const_iterator obj = eventVec.begin(); obj!= eventVec.end(); ++obj){
@@ -162,8 +156,8 @@ metsig::significanceAlgo::significance(double &met_r, double &met_phi, double &m
 
   //--- Temporary variables ---//
 
-  TMatrixD v_tot(2,2);
-  TVectorD metvec(2);
+  reco::METCovMatrix v_tot;
+  ROOT::Math::SVector<double, 2> metvec;
 
  //--- Initialize sum of rotated covariance matrices ---//  
   v_tot=signifmatrix_;
@@ -178,7 +172,9 @@ metsig::significanceAlgo::significance(double &met_r, double &met_phi, double &m
 
   // one other option: if particles cancel there could be small numbers.
   // this check fixes this, added by F.Blekman
-  if(fabs(v_tot.Determinant())<0.000001)
+  double det=0;
+  v_tot.Det(det);
+  if(fabs(det)<0.000001)
     return -1;
 
 
@@ -190,7 +186,7 @@ metsig::significanceAlgo::significance(double &met_r, double &met_phi, double &m
 
 
   metvec(0) = xmet_; metvec(1) = ymet_;
-  double lnSignificance = metvec * (v_tot * metvec);
+  double lnSignificance = ROOT::Math::Dot(metvec, (v_tot * metvec) );
 
   //  v_tot.Invert();
   //  std::cout << "INVERTED AGAIN:\n"<< v_tot(0,0) << "," << v_tot(0,1) << "\n" << v_tot(1,0) << "," << v_tot(1,1) << std::endl;