JetMET/src/AlphaT.cc

0001 #include <numeric>
0002 #include <cmath>
0003 #include "HLTrigger/JetMET/interface/AlphaT.h"
0004
0005 double AlphaT::value_(std::vector<bool>* jet_sign) const {
0006   // Clear pseudo-jet container
0007   if (jet_sign) {
0008     jet_sign->clear();
0009     jet_sign->resize(et_.size());
0010   }
0011
0012   // check the size of the input collection
0013   if (et_.empty())
0014     // empty jet collection, return AlphaT = 0
0015     return 0.;
0016
0017   if (et_.size() > (unsigned int)std::numeric_limits<unsigned int>::digits)
0018     // too many jets, return AlphaT = a very large number
0019     return std::numeric_limits<double>::max();
0020
0021   // Momentum sums in transverse plane
0022   const double sum_et = std::accumulate(et_.begin(), et_.end(), 0.);
0023   const double sum_px = std::accumulate(px_.begin(), px_.end(), 0.);
0024   const double sum_py = std::accumulate(py_.begin(), py_.end(), 0.);
0025
0026   // Minimum Delta Et for two pseudo-jets
0027   double min_delta_sum_et = sum_et;
0028
0029   if (setDHtZero_) {
0030     min_delta_sum_et = 0.;
0031   } else {
0032     for (unsigned int i = 0; i < (1U << (et_.size() - 1)); i++) {  //@@ iterate through different combinations
0033       double delta_sum_et = 0.;
0034       for (unsigned int j = 0; j < et_.size(); ++j) {  //@@ iterate through jets
0035         if (i & (1U << j))
0036           delta_sum_et -= et_[j];
0037         else
0038           delta_sum_et += et_[j];
0039       }
0040       delta_sum_et = std::abs(delta_sum_et);
0041       if (delta_sum_et < min_delta_sum_et) {
0042         min_delta_sum_et = delta_sum_et;
0043         if (jet_sign) {
0044           for (unsigned int j = 0; j < et_.size(); ++j)
0045             (*jet_sign)[j] = ((i & (1U << j)) == 0);
0046         }
0047       }
0048     }
0049   }
0050   // Alpha_T
0051   return (0.5 * (sum_et - min_delta_sum_et) / std::sqrt(sum_et * sum_et - (sum_px * sum_px + sum_py * sum_py)));
0052 }