CMS 3D CMS Logo

HEPTopTaggerWrapperV2.cc
Go to the documentation of this file.
1 // 2011 Christopher Vermilion
2 //
3 //----------------------------------------------------------------------
4 // This file is free software; you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation; either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This file is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // The GNU General Public License is available at
15 // http://www.gnu.org/licenses/gpl.html or you can write to the Free Software
16 // Foundation, Inc.:
17 // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 //----------------------------------------------------------------------
19 
21 
22 #include "fastjet/Error.hh"
23 #include "fastjet/JetDefinition.hh"
24 #include "fastjet/ClusterSequence.hh"
25 #include "fastjet/PseudoJet.hh"
26 #include "fastjet/tools/Pruner.hh"
27 #include "fastjet/tools/Filter.hh"
28 
29 #include <cmath>
30 #include <limits>
31 #include <cassert>
32 using namespace std;
33 
35 
36 
37 FASTJET_BEGIN_NAMESPACE
38 
39 
40 // Expected R_min for tops (as function of filtered initial fatjet pT in GeV (using CA, R=0.2, n=10)
41 // From ttbar sample, phys14, n20, bx25
42 // matched to hadronically decaying top with delta R < 1.2 and true top pT > 200
43 // Cuts are: fW < 0.175 and m_top = 120..220
44 // Input objects are packed pfCandidates with CHS
45 double R_min_expected_function(double x){
46 
47 
48  if (x>700)
49  x=700;
50 
51  double A = -9.42052;
52  double B = 0.202773;
53  double C = 4498.45;
54  double D = -1.05737e+06;
55  double E = 9.95494e+07;
56 
57  return A+B*sqrt(x)+C/x+D/(x*x)+E/(x*x*x);
58 }
59 
60 
61 
62 //------------------------------------------------------------------------
63 // returns the tagged PseudoJet if successful, 0 otherwise
64 // - jet the PseudoJet to tag
65 PseudoJet HEPTopTaggerV2::result(const PseudoJet & jet) const{
66 
67  // make sure that there is a "regular" cluster sequence associated
68  // with the jet. Note that we also check it is valid (to avoid a
69  // more criptic error later on)
70  if (!jet.has_valid_cluster_sequence()){
71  throw Error("HEPTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
72  }
73 
74  external::HEPTopTaggerV2 tagger(jet);
75 
76  external::HEPTopTaggerV2 best_tagger;
77 
78  // translate the massRatioWidth (which should be the half-width given in %)
79  // to values useful for the A-shape cuts
80  double mw_over_mt = 80.4/172.3;
81  double ratio_min = mw_over_mt * (100.-massRatioWidth_)/100.;
82  double ratio_max = mw_over_mt * (100.+massRatioWidth_)/100.;
83 
84  // Unclustering, Filtering & Subjet Settings
85  tagger.set_max_subjet_mass(subjetMass_);
86  tagger.set_mass_drop_threshold(muCut_);
87  tagger.set_filtering_R(filtR_);
88  tagger.set_filtering_n(filtN_);
89  tagger.set_filtering_minpt_subjet(minSubjetPt_);
90 
91  // Optimal R
92  tagger.do_optimalR(DoOptimalR_);
93  tagger.set_optimalR_reject_minimum(optRrejectMin_);
94 
95  // How to select among candidates
96  tagger.set_mode((external::Mode)mode_);
97 
98  // Requirements to accept a candidate
99  tagger.set_top_minpt(minCandPt_);
100  tagger.set_top_mass_range(minCandMass_, maxCandMass_);
101  tagger.set_mass_ratio_cut(minM23Cut_, minM13Cut_, maxM13Cut_);
102  tagger.set_mass_ratio_range(ratio_min, ratio_max);
103 
104  // Set function to calculate R_min_expected
106 
107 
108  double qweight = -1;
109  double qepsilon = -1;
110  double qsigmaM = -1;
111 
112  if (DoQjets_){
113 
114  int niter(100);
115  double q_zcut(0.1);
116  double q_dcut_fctr(0.5);
117  double q_exp_min(0.);
118  double q_exp_max(0.);
119  double q_rigidity(0.1);
120  double q_truncation_fctr(0.0);
121 
122  double weight_q1 = -1.;
123  double m_sum = 0.;
124  double m2_sum = 0.;
125  int qtags = 0;
126 
127  tagger.set_qjets(q_zcut,
128  q_dcut_fctr,
129  q_exp_min,
130  q_exp_max,
131  q_rigidity,
132  q_truncation_fctr);
133  tagger.set_qjets_rng(engine_);
134  tagger.do_qjets(true);
135  tagger.run();
136 
137  for (int iq = 0; iq < niter; iq++) {
138  tagger.run();
139  if (tagger.is_tagged()) {
140  qtags++;
141  m_sum += tagger.t().m();
142  m2_sum += tagger.t().m() * tagger.t().m();
143  if (tagger.q_weight() > weight_q1) {
144  best_tagger = tagger;
145  weight_q1=tagger.q_weight();
146  }
147  }
148  }
149 
150  tagger = best_tagger;
151  qweight = weight_q1;
152  qepsilon = float(qtags)/float(niter);
153 
154  // calculate width of tagged mass distribution if we have at least one candidate
155  if (qtags > 0){
156  double mean_m = m_sum / qtags;
157  double mean_m2 = m2_sum / qtags;
158  qsigmaM = sqrt(mean_m2 - mean_m*mean_m);
159  }
160  }
161  else{
162  tagger.run();
163  }
164 
165  // Requires:
166  // - top mass window
167  // - mass ratio cuts
168  // - minimal candidate pT
169  // If this is not intended: use loose top mass and ratio windows
170  if (!tagger.is_tagged())
171  return PseudoJet();
172 
173  // create the result and its structure
174  const JetDefinition::Recombiner *rec
175  = jet.associated_cluster_sequence()->jet_def().recombiner();
176 
177  const vector<PseudoJet>& subjets = tagger.top_subjets();
178  assert(subjets.size() == 3);
179 
180  PseudoJet non_W = subjets[0];
181  PseudoJet W1 = subjets[1];
182  PseudoJet W2 = subjets[2];
183  PseudoJet W = join(subjets[1], subjets[2], *rec);
184 
185  PseudoJet result = join<HEPTopTaggerV2Structure>( W1, W2, non_W, *rec);
186  HEPTopTaggerV2Structure *s = (HEPTopTaggerV2Structure*) result.structure_non_const_ptr();
187 
188  s->_fj_mass = jet.m();
189  s->_fj_pt = jet.perp();
190  s->_fj_eta = jet.eta();
191  s->_fj_phi = jet.phi();
192 
193  s->_top_mass = tagger.t().m();
194  s->_pruned_mass = tagger.pruned_mass();
195  s->_unfiltered_mass = tagger.unfiltered_mass();
196  s->_fRec = tagger.f_rec();
197  s->_mass_ratio_passed = tagger.is_masscut_passed();
198 
199  if (DoOptimalR_){
200  s->_tau1Unfiltered = tagger.nsub_unfiltered(1);
201  s->_tau2Unfiltered = tagger.nsub_unfiltered(2);
202  s->_tau3Unfiltered = tagger.nsub_unfiltered(3);
203  s->_tau1Filtered = tagger.nsub_filtered(1);
204  s->_tau2Filtered = tagger.nsub_filtered(2);
205  s->_tau3Filtered = tagger.nsub_filtered(3);
206  }
207 
208  s->_qweight = qweight;
209  s->_qepsilon = qepsilon;
210  s->_qsigmaM = qsigmaM;
211 
212 
213  if (DoOptimalR_){
214  s->_ropt = tagger.Ropt();
215  s->_roptCalc = tagger.Ropt_calc();
216  s->_ptForRoptCalc = tagger.pt_for_Ropt_calc();
217  }
218  else {
219  s->_ropt = -1;
220  s->_roptCalc = -1;
221  s->_ptForRoptCalc = -1;
222  }
223 
224  // Removed selectors as all cuts are applied in HTT
225 
226  return result;
227 }
228 
229 FASTJET_END_NAMESPACE
edm::ErrorSummaryEntry Error
void set_top_mass_range(double xmin, double xmax)
void do_qjets(bool qjets)
void set_optimalR_reject_minimum(bool x)
const PseudoJet & t() const
void set_qjets_rng(CLHEP::HepRandomEngine *engine)
void set_filtering_minpt_subjet(double x)
void set_top_minpt(double x)
void set_mass_ratio_range(double rmin, double rmax)
void do_optimalR(bool optimalR)
double nsub_filtered(int order, fastjet::contrib::Njettiness::AxesMode axes=fastjet::contrib::Njettiness::kt_axes, double beta=1., double R0=1.)
double pt_for_Ropt_calc() const
T sqrt(T t)
Definition: SSEVec.h:18
void set_qjets(double q_zcut, double q_dcut_fctr, double q_exp_min, double q_exp_max, double q_rigidity, double q_truncation_fctr)
double unfiltered_mass() const
void set_filtering_R(double Rfilt)
double nsub_unfiltered(int order, fastjet::contrib::Njettiness::AxesMode axes=fastjet::contrib::Njettiness::kt_axes, double beta=1., double R0=1.)
void set_filtering_n(unsigned nfilt)
double Ropt_calc() const
void set_mode(enum Mode mode)
static const std::string B
const std::vector< PseudoJet > & top_subjets() const
PseudoJet result(const PseudoJet &jet) const override
void set_mass_drop_threshold(double x)
void set_max_subjet_mass(double x)
static std::string join(char **cmd)
Definition: RemoteFile.cc:18
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:152
bool is_masscut_passed() const
double pruned_mass() const
void set_mass_ratio_cut(double m23cut, double m13cutmin, double m13cutmax)
FASTJET_BEGIN_NAMESPACE double R_min_expected_function(double x)
void set_optimalR_calc_fun(double(*f)(double))