METAlgorithms/src/significanceAlgo.cc

0001 #include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
0002 #include <iostream>
0003 #include <string>
0004 #include "TMath.h"
0005 #include <cmath>
0006 #include "TROOT.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 // -*- C++ -*-
0010 //
0011 // Package:    METAlgorithms
0012 // Class:      SignAlgoResolutions
0013 //
0014 /**\class METSignificance SignAlgoResolutions.cc RecoMET/METAlgorithms/src/SignAlgoResolutions.cc
0015
0016  Description: <one line class summary>
0017
0018  Implementation:
0019      <Notes on implementation>
0020 */
0021 //
0022 // Original Author:  Kyle Story, Freya Blekman (Cornell University)
0023 //         Created:  Fri Apr 18 11:58:33 CEST 2008
0024 //
0025 //
0026
0027 metsig::significanceAlgo::significanceAlgo()
0028     :  //  eventVec_(0),
0029       set_worker_(0),
0030       xmet_(0),
0031       ymet_(0)
0032
0033 {
0034   //  std::cout << "in constructor ! " << std::endl;
0035   signifmatrix_(0, 0) = signifmatrix_(1, 0) = signifmatrix_(0, 1) = signifmatrix_(1, 1) = 0;
0036 }
0037
0038 //******* 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)!
0039 const void metsig::significanceAlgo::addSignifMatrix(const reco::METCovMatrix &input) {
0040   signifmatrix_ += input;
0041   return;
0042 }
0043 ////////////////////////
0044 /// reset the signficance matrix (this is the most likely case), so that the vector sum can be continued
0045
0046 const void metsig::significanceAlgo::setSignifMatrix(const reco::METCovMatrix &input,
0047                                                      const double &met_r,
0048                                                      const double &met_phi,
0049                                                      const double &met_set) {
0050   signifmatrix_ = input;
0051   set_worker_ = met_set;
0052   xmet_ = met_r * cos(met_phi);
0053   ymet_ = met_r * sin(met_phi);
0054
0055   return;
0056 }
0057
0058 // ********destructor ********
0059
0060 metsig::significanceAlgo::~significanceAlgo() {}
0061
0062 // //*** rotate a 2D matrix by angle theta **********************//
0063
0064 // void
0065 // metsig::significanceAlgo::rotateMatrix( Double_t theta, reco::METCovMatrix &v)
0066 // {
0067 //   // I suggest not using this to rotate trivial matrices.
0068 //   reco::METCovMatrix r;
0069 //   reco::METCovMatrix rInv;
0070
0071 //   r(0,0) = cos(theta); r(0,1) = sin(theta); r(1,0) = -sin(theta); r(1,1) = cos(theta);
0072 //   rInv = r;
0073 //   rInv.Invert();
0074 //   //-- Rotate v --//
0075 //   v = rInv * v * r;
0076 // }
0077 //************************************************************//
0078
0079 const void metsig::significanceAlgo::subtractObjects(const std::vector<metsig::SigInputObj> &eventVec) {
0080   reco::METCovMatrix v_tot = signifmatrix_;
0081   //--- Loop over physics objects in the event ---//
0082   //  for(unsigned int objnum=1; objnum < EventVec.size(); objnum++ ) {
0083   for (std::vector<SigInputObj>::const_iterator obj = eventVec.begin(); obj != eventVec.end(); ++obj) {
0084     double et_tmp = obj->get_energy();
0085     double phi_tmp = obj->get_phi();
0086     double sigma_et = obj->get_sigma_e();
0087     double sigma_tan = obj->get_sigma_tan();
0088
0089     double cosphi = cos(phi_tmp);
0090     double sinphi = sin(phi_tmp);
0091     double xval = et_tmp * cosphi;
0092     double yval = et_tmp * sinphi;
0093     xmet_ += xval;
0094     ymet_ += yval;
0095     set_worker_ -= et_tmp;
0096
0097     double sigma0_2 = sigma_et * sigma_et;
0098     double sigma1_2 = sigma_tan * sigma_tan;
0099
0100     v_tot(0, 0) -= sigma0_2 * cosphi * cosphi + sigma1_2 * sinphi * sinphi;
0101     v_tot(0, 1) -= cosphi * sinphi * (sigma0_2 - sigma1_2);
0102     v_tot(1, 0) -= cosphi * sinphi * (sigma0_2 - sigma1_2);
0103     v_tot(1, 1) -= sigma1_2 * cosphi * cosphi + sigma0_2 * sinphi * sinphi;
0104   }
0105   signifmatrix_ = v_tot;
0106 }
0107 //************************************************************//
0108
0109 const void metsig::significanceAlgo::addObjects(const std::vector<metsig::SigInputObj> &eventVec) {
0110   reco::METCovMatrix v_tot = signifmatrix_;
0111   //--- Loop over physics objects in the event ---//
0112   //  for(unsigned int objnum=1; objnum < EventVec.size(); objnum++ ) {
0113   for (std::vector<SigInputObj>::const_iterator obj = eventVec.begin(); obj != eventVec.end(); ++obj) {
0114     double et_tmp = obj->get_energy();
0115     double phi_tmp = obj->get_phi();
0116     double sigma_et = obj->get_sigma_e();
0117     double sigma_tan = obj->get_sigma_tan();
0118
0119     double cosphi = cos(phi_tmp);
0120     double sinphi = sin(phi_tmp);
0121     double xval = et_tmp * cosphi;
0122     double yval = et_tmp * sinphi;
0123     xmet_ -= xval;
0124     ymet_ -= yval;
0125     set_worker_ += et_tmp;
0126
0127     double sigma0_2 = sigma_et * sigma_et;
0128     double sigma1_2 = sigma_tan * sigma_tan;
0129
0130     v_tot(0, 0) += sigma0_2 * cosphi * cosphi + sigma1_2 * sinphi * sinphi;
0131     v_tot(0, 1) += cosphi * sinphi * (sigma0_2 - sigma1_2);
0132     v_tot(1, 0) += cosphi * sinphi * (sigma0_2 - sigma1_2);
0133     v_tot(1, 1) += sigma1_2 * cosphi * cosphi + sigma0_2 * sinphi * sinphi;
0134   }
0135   signifmatrix_ = v_tot;
0136 }
0137
0138 //************************************************************//
0139 const double metsig::significanceAlgo::significance(double &met_r, double &met_phi, double &met_set) {
0140   if (signifmatrix_(0, 0) == 0 && signifmatrix_(1, 1) == 0 && signifmatrix_(1, 0) == 0 && signifmatrix_(0, 1) == 0) {
0141     //edm::LogWarning("SignCaloSpecificAlgo") << "Event Vector is empty!  Return significance -1";
0142     return (-1);
0143   }
0144
0145   //--- Temporary variables ---//
0146
0147   reco::METCovMatrix v_tot;
0148   ROOT::Math::SVector<double, 2> metvec;
0149
0150   //--- Initialize sum of rotated covariance matrices ---//
0151   v_tot = signifmatrix_;
0152   //  std::cout << "INPUT:\n"<< v_tot(0,0) << "," << v_tot(0,1) << "\n" << v_tot(1,0) << "," << v_tot(1,1) << std::endl;
0153
0154   //--- Calculate magnitude and angle of MET, store in returned variables ---//
0155   met_r = sqrt(xmet_ * xmet_ + ymet_ * ymet_);
0156   met_set = set_worker_;
0157   met_phi = TMath::ATan2(ymet_, xmet_);
0158
0159   // one other option: if particles cancel there could be small numbers.
0160   // this check fixes this, added by F.Blekman
0161   double det = 0;
0162   v_tot.Det(det);
0163   if (fabs(det) < 0.000001)
0164     return -1;
0165
0166   // save matrix into object:
0167   //  std::cout << "SAVED:\n"<< v_tot(0,0) << "," << v_tot(0,1) << "\n" << v_tot(1,0) << "," << v_tot(1,1) << std::endl;
0168   v_tot.Invert();
0169   //  std::cout << "INVERTED:\n"<< v_tot(0,0) << "," << v_tot(0,1) << "\n" << v_tot(1,0) << "," << v_tot(1,1) << std::endl;
0170
0171   metvec(0) = xmet_;
0172   metvec(1) = ymet_;
0173   double lnSignificance = ROOT::Math::Dot(metvec, (v_tot * metvec));
0174
0175   //  v_tot.Invert();
0176   //  std::cout << "INVERTED AGAIN:\n"<< v_tot(0,0) << "," << v_tot(0,1) << "\n" << v_tot(1,0) << "," << v_tot(1,1) << std::endl;
0177   return lnSignificance;
0178 }
0179 //*** End of Significance calculation ********************************//