Heppy/src/Apc.cc

0001 #include "PhysicsTools/Heppy/interface/Apc.h"
0002 #include "DataFormats/Math/interface/deltaPhi.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 #include <cmath>
0006 #include <numeric>
0007 #include <vector>
0008 #include <algorithm>
0009 #include <cstdlib>
0010
0011 namespace heppy {
0012
0013   double Apc::getApcJetMetMin(const std::vector<double>& et,
0014                               const std::vector<double>& px,
0015                               const std::vector<double>& py,
0016                               const double metx,
0017                               const double mety) {
0018     if (et.empty())
0019       return -1.;
0020
0021     // Momentum sums in transverse plane
0022     const double ht = accumulate(et.begin(), et.end(), 0.);
0023
0024     // jets are pt-sorted
0025     double apcjetmetmin(0.);
0026
0027     std::vector<double> apcjetvector;
0028     std::vector<double> apcjetmetvector;
0029     for (size_t j = 0; j < et.size(); j++) {
0030       apcjetvector.push_back(0.);
0031       apcjetmetvector.push_back(0.);
0032       double jet_phi_j = atan2(py[j], px[j]);
0033       for (size_t i = 0; i < et.size(); i++) {
0034         double jet_phi_i = atan2(py[i], px[i]);
0035         double dphi_jet = fabs(deltaPhi(jet_phi_i, jet_phi_j));
0036         double met_phi = atan2(mety, metx);
0037         double dphimet = fabs(deltaPhi(jet_phi_i, met_phi));
0038
0039         apcjetvector.back() += et[i] * cos(dphi_jet / 2.0);
0040         apcjetmetvector.back() += et[i] * cos(dphi_jet / 2.0) * sin(dphimet / 2.0);
0041       }
0042     }
0043     if (!apcjetvector.empty() && !apcjetmetvector.empty()) {
0044       apcjetmetmin = *min_element(apcjetmetvector.begin(), apcjetmetvector.end());
0045     }
0046
0047     if (ht != 0)
0048       return apcjetmetmin / ht;
0049     else
0050       return -1.;
0051   }
0052 }  // namespace heppy