JetAlgorithms/src/HEPTopTaggerWrapperV2.cc

0001 // 2011 Christopher Vermilion
0002 //
0003 //----------------------------------------------------------------------
0004 //  This file is free software; you can redistribute it and/or modify
0005 //  it under the terms of the GNU General Public License as published by
0006 //  the Free Software Foundation; either version 3 of the License, or
0007 //  (at your option) any later version.
0008 //
0009 //  This file is distributed in the hope that it will be useful,
0010 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0011 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0012 //  GNU General Public License for more details.
0013 //
0014 //  The GNU General Public License is available at
0015 //  http://www.gnu.org/licenses/gpl.html or you can write to the Free Software
0016 //  Foundation, Inc.:
0017 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
0018 //----------------------------------------------------------------------
0019
0020 #include "RecoJets/JetAlgorithms/interface/HEPTopTaggerWrapperV2.h"
0021
0022 #include "fastjet/Error.hh"
0023 #include "fastjet/JetDefinition.hh"
0024 #include "fastjet/ClusterSequence.hh"
0025 #include "fastjet/PseudoJet.hh"
0026 #include "fastjet/tools/Pruner.hh"
0027 #include "fastjet/tools/Filter.hh"
0028
0029 #include <cmath>
0030 #include <limits>
0031 #include <cassert>
0032 using namespace std;
0033
0034 #include "RecoJets/JetAlgorithms/interface/HEPTopTaggerV2.h"
0035
0036 FASTJET_BEGIN_NAMESPACE
0037
0038 // Expected R_min for tops (as function of filtered initial fatjet pT in GeV (using CA, R=0.2, n=10)
0039 // From ttbar sample, phys14, n20, bx25
0040 // matched to hadronically decaying top with delta R < 1.2 and true top pT > 200
0041 // Cuts are: fW < 0.175 and  m_top = 120..220
0042 // Input objects are packed pfCandidates with CHS
0043 double R_min_expected_function(double x) {
0044   if (x > 700)
0045     x = 700;
0046
0047   double A = -9.42052;
0048   double B = 0.202773;
0049   double C = 4498.45;
0050   double D = -1.05737e+06;
0051   double E = 9.95494e+07;
0052
0053   return A + B * sqrt(x) + C / x + D / (x * x) + E / (x * x * x);
0054 }
0055
0056 //------------------------------------------------------------------------
0057 // returns the tagged PseudoJet if successful, 0 otherwise
0058 //  - jet   the PseudoJet to tag
0059 PseudoJet HEPTopTaggerV2::result(const PseudoJet &jet) const {
0060   // make sure that there is a "regular" cluster sequence associated
0061   // with the jet. Note that we also check it is valid (to avoid a
0062   // more criptic error later on)
0063   if (!jet.has_valid_cluster_sequence()) {
0064     throw Error("HEPTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
0065   }
0066
0067   external::HEPTopTaggerV2 tagger(jet);
0068
0069   external::HEPTopTaggerV2 best_tagger;
0070
0071   // translate the massRatioWidth (which should be the half-width given in %)
0072   // to values useful for the A-shape cuts
0073   double mw_over_mt = 80.4 / 172.3;
0074   double ratio_min = mw_over_mt * (100. - massRatioWidth_) / 100.;
0075   double ratio_max = mw_over_mt * (100. + massRatioWidth_) / 100.;
0076
0077   // Unclustering, Filtering & Subjet Settings
0078   tagger.set_max_subjet_mass(subjetMass_);
0079   tagger.set_mass_drop_threshold(muCut_);
0080   tagger.set_filtering_R(filtR_);
0081   tagger.set_filtering_n(filtN_);
0082   tagger.set_filtering_minpt_subjet(minSubjetPt_);
0083
0084   // Optimal R
0085   tagger.do_optimalR(DoOptimalR_);
0086   tagger.set_optimalR_reject_minimum(optRrejectMin_);
0087
0088   // How to select among candidates
0089   tagger.set_mode((external::Mode)mode_);
0090
0091   // Requirements to accept a candidate
0092   tagger.set_top_minpt(minCandPt_);
0093   tagger.set_top_mass_range(minCandMass_, maxCandMass_);
0094   tagger.set_mass_ratio_cut(minM23Cut_, minM13Cut_, maxM13Cut_);
0095   tagger.set_mass_ratio_range(ratio_min, ratio_max);
0096
0097   // Set function to calculate R_min_expected
0098   tagger.set_optimalR_calc_fun(R_min_expected_function);
0099
0100   double qweight = -1;
0101   double qepsilon = -1;
0102   double qsigmaM = -1;
0103
0104   if (DoQjets_) {
0105     int niter(100);
0106     double q_zcut(0.1);
0107     double q_dcut_fctr(0.5);
0108     double q_exp_min(0.);
0109     double q_exp_max(0.);
0110     double q_rigidity(0.1);
0111     double q_truncation_fctr(0.0);
0112
0113     double weight_q1 = -1.;
0114     double m_sum = 0.;
0115     double m2_sum = 0.;
0116     int qtags = 0;
0117
0118     tagger.set_qjets(q_zcut, q_dcut_fctr, q_exp_min, q_exp_max, q_rigidity, q_truncation_fctr);
0119     tagger.set_qjets_rng(engine_);
0120     tagger.do_qjets(true);
0121     tagger.run();
0122
0123     for (int iq = 0; iq < niter; iq++) {
0124       tagger.run();
0125       if (tagger.is_tagged()) {
0126         qtags++;
0127         m_sum += tagger.t().m();
0128         m2_sum += tagger.t().m() * tagger.t().m();
0129         if (tagger.q_weight() > weight_q1) {
0130           best_tagger = tagger;
0131           weight_q1 = tagger.q_weight();
0132         }
0133       }
0134     }
0135
0136     tagger = best_tagger;
0137     qweight = weight_q1;
0138     qepsilon = float(qtags) / float(niter);
0139
0140     // calculate width of tagged mass distribution if we have at least one candidate
0141     if (qtags > 0) {
0142       double mean_m = m_sum / qtags;
0143       double mean_m2 = m2_sum / qtags;
0144       qsigmaM = sqrt(mean_m2 - mean_m * mean_m);
0145     }
0146   } else {
0147     tagger.run();
0148   }
0149
0150   // Requires:
0151   //   - top mass window
0152   //   - mass ratio cuts
0153   //   - minimal candidate pT
0154   // If this is not intended: use loose top mass and ratio windows
0155   if (!tagger.is_tagged())
0156     return PseudoJet();
0157
0158   // create the result and its structure
0159   const JetDefinition::Recombiner *rec = jet.associated_cluster_sequence()->jet_def().recombiner();
0160
0161   const vector<PseudoJet> &subjets = tagger.top_subjets();
0162   assert(subjets.size() == 3);
0163
0164   PseudoJet non_W = subjets[0];
0165   PseudoJet W1 = subjets[1];
0166   PseudoJet W2 = subjets[2];
0167   PseudoJet W = join(subjets[1], subjets[2], *rec);
0168
0169   PseudoJet result = join<HEPTopTaggerV2Structure>(W1, W2, non_W, *rec);
0170   HEPTopTaggerV2Structure *s = (HEPTopTaggerV2Structure *)result.structure_non_const_ptr();
0171
0172   s->_fj_mass = jet.m();
0173   s->_fj_pt = jet.perp();
0174   s->_fj_eta = jet.eta();
0175   s->_fj_phi = jet.phi();
0176
0177   s->_top_mass = tagger.t().m();
0178   s->_pruned_mass = tagger.pruned_mass();
0179   s->_unfiltered_mass = tagger.unfiltered_mass();
0180   s->_fRec = tagger.f_rec();
0181   s->_mass_ratio_passed = tagger.is_masscut_passed();
0182
0183   if (DoOptimalR_) {
0184     s->_tau1Unfiltered = tagger.nsub_unfiltered(1);
0185     s->_tau2Unfiltered = tagger.nsub_unfiltered(2);
0186     s->_tau3Unfiltered = tagger.nsub_unfiltered(3);
0187     s->_tau1Filtered = tagger.nsub_filtered(1);
0188     s->_tau2Filtered = tagger.nsub_filtered(2);
0189     s->_tau3Filtered = tagger.nsub_filtered(3);
0190   }
0191
0192   s->_qweight = qweight;
0193   s->_qepsilon = qepsilon;
0194   s->_qsigmaM = qsigmaM;
0195
0196   if (DoOptimalR_) {
0197     s->_ropt = tagger.Ropt();
0198     s->_roptCalc = tagger.Ropt_calc();
0199     s->_ptForRoptCalc = tagger.pt_for_Ropt_calc();
0200   } else {
0201     s->_ropt = -1;
0202     s->_roptCalc = -1;
0203     s->_ptForRoptCalc = -1;
0204   }
0205
0206   // Removed selectors as all cuts are applied in HTT
0207
0208   return result;
0209 }
0210
0211 FASTJET_END_NAMESPACE