Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:00:57

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 ********************************//