Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:25 

0001 #include "RecoJets/JetAlgorithms/interface/HEPTopTaggerV2.h"
0002 
0003 // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
0004 namespace external {
0005 
0006   //optimal_R fit
0007   double R_opt_calc_funct(double pt_filt) { return 327. / pt_filt; }
0008 
0009   void HEPTopTaggerV2_fixed_R::print_banner() {
0010     if (!_first_time) {
0011       return;
0012     }
0013     _first_time = false;
0014 
0015     edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
0016     edm::LogInfo("HEPTopTaggerV2") << "#                   HEPTopTaggerV2 - under construction                      \n";
0017     edm::LogInfo("HEPTopTaggerV2") << "#                                                                          \n";
0018     edm::LogInfo("HEPTopTaggerV2") << "# Please cite JHEP 1010 (2010) 078 [arXiv:1006.2833 [hep-ph]]              \n";
0019     edm::LogInfo("HEPTopTaggerV2") << "# and Phys.Rev. D89 (2014) 074047 [arXiv:1312.1504 [hep-ph]]               \n";
0020     edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
0021   }
0022 
0023   //pt wrt a reference vector
0024   double HEPTopTaggerV2_fixed_R::perp(const fastjet::PseudoJet& vec, const fastjet::PseudoJet& ref) {
0025     double ref_ref = ref.px() * ref.px() + ref.py() * ref.py() + ref.pz() * ref.pz();
0026     double vec_ref = vec.px() * ref.px() + vec.py() * ref.py() + vec.pz() * ref.pz();
0027     double per_per = vec.px() * vec.px() + vec.py() * vec.py() + vec.pz() * vec.pz();
0028     if (ref_ref > 0.)
0029       per_per -= vec_ref * vec_ref / ref_ref;
0030     if (per_per < 0.)
0031       per_per = 0.;
0032     return sqrt(per_per);
0033   }
0034 
0035   //modified Jade distance
0036   double HEPTopTaggerV2_fixed_R::djademod(const fastjet::PseudoJet& subjet_i,
0037                                           const fastjet::PseudoJet& subjet_j,
0038                                           const fastjet::PseudoJet& ref) {
0039     double dj = -1.0;
0040     double delta_phi = subjet_i.delta_phi_to(subjet_j);
0041     double delta_eta = subjet_i.eta() - subjet_j.eta();
0042     double delta_R = sqrt(delta_eta * delta_eta + delta_phi * delta_phi);
0043     dj = perp(subjet_i, ref) * perp(subjet_j, ref) * pow(delta_R, 4.);
0044     return dj;
0045   }
0046 
0047   //minimal |(m_ij / m_123) / (m_w/ m_t) - 1|
0048   double HEPTopTaggerV2_fixed_R::f_rec() {
0049     double m12 = (_top_subs[0] + _top_subs[1]).m();
0050     double m13 = (_top_subs[0] + _top_subs[2]).m();
0051     double m23 = (_top_subs[1] + _top_subs[2]).m();
0052     double m123 = (_top_subs[0] + _top_subs[1] + _top_subs[2]).m();
0053 
0054     double fw12 = fabs((m12 / m123) / (_mwmass / _mtmass) - 1);
0055     double fw13 = fabs((m13 / m123) / (_mwmass / _mtmass) - 1);
0056     double fw23 = fabs((m23 / m123) / (_mwmass / _mtmass) - 1);
0057 
0058     return std::min(fw12, std::min(fw13, fw23));
0059   }
0060 
0061   //find hard substructures
0062   void HEPTopTaggerV2_fixed_R::FindHardSubst(const PseudoJet& this_jet, std::vector<fastjet::PseudoJet>& t_parts) {
0063     PseudoJet parent1(0, 0, 0, 0), parent2(0, 0, 0, 0);
0064     if (this_jet.m() < _max_subjet_mass || !this_jet.validated_cs()->has_parents(this_jet, parent1, parent2)) {
0065       t_parts.push_back(this_jet);
0066     } else {
0067       if (parent1.m() < parent2.m())
0068         std::swap(parent1, parent2);
0069       FindHardSubst(parent1, t_parts);
0070       if (parent1.m() < _mass_drop_threshold * this_jet.m())
0071         FindHardSubst(parent2, t_parts);
0072     }
0073   }
0074 
0075   //store subjets as vector<PseudoJet> with [0]->b [1]->W-jet 1 [2]->W-jet 2
0076   void HEPTopTaggerV2_fixed_R::store_topsubjets(const std::vector<PseudoJet>& top_subs) {
0077     _top_subjets.resize(0);
0078     double m12 = (top_subs[0] + top_subs[1]).m();
0079     double m13 = (top_subs[0] + top_subs[2]).m();
0080     double m23 = (top_subs[1] + top_subs[2]).m();
0081     double dm12 = fabs(m12 - _mwmass);
0082     double dm13 = fabs(m13 - _mwmass);
0083     double dm23 = fabs(m23 - _mwmass);
0084 
0085     if (dm23 <= dm12 && dm23 <= dm13) {
0086       _top_subjets.push_back(top_subs[0]);
0087       _top_subjets.push_back(top_subs[1]);
0088       _top_subjets.push_back(top_subs[2]);
0089     } else if (dm13 <= dm12 && dm13 < dm23) {
0090       _top_subjets.push_back(top_subs[1]);
0091       _top_subjets.push_back(top_subs[0]);
0092       _top_subjets.push_back(top_subs[2]);
0093     } else if (dm12 < dm23 && dm12 < dm13) {
0094       _top_subjets.push_back(top_subs[2]);
0095       _top_subjets.push_back(top_subs[0]);
0096       _top_subjets.push_back(top_subs[1]);
0097     }
0098     _W = _top_subjets[1] + _top_subjets[2];
0099     return;
0100   }
0101 
0102   //check mass plane cuts
0103   bool HEPTopTaggerV2_fixed_R::check_mass_criteria(const std::vector<PseudoJet>& top_subs) const {
0104     bool is_passed = false;
0105     double m12 = (top_subs[0] + top_subs[1]).m();
0106     double m13 = (top_subs[0] + top_subs[2]).m();
0107     double m23 = (top_subs[1] + top_subs[2]).m();
0108     double m123 = (top_subs[0] + top_subs[1] + top_subs[2]).m();
0109     double atan1312 = atan(m13 / m12);
0110     double m23_over_m123 = m23 / m123;
0111     double m23_over_m123_square = m23_over_m123 * m23_over_m123;
0112     double rmin_square = _rmin * _rmin;
0113     double rmax_square = _rmax * _rmax;
0114     double m13m12_square_p1 = (1 + (m13 / m12) * (m13 / m12));
0115     double m12m13_square_p1 = (1 + (m12 / m13) * (m12 / m13));
0116     if ((atan1312 > _m13cutmin && _m13cutmax > atan1312 && (m23_over_m123 > _rmin && _rmax > m23_over_m123)) ||
0117         ((m23_over_m123_square < 1 - rmin_square * m13m12_square_p1) &&
0118          (m23_over_m123_square > 1 - rmax_square * m13m12_square_p1) && (m23_over_m123 > _m23cut)) ||
0119         ((m23_over_m123_square < 1 - rmin_square * m12m13_square_p1) &&
0120          (m23_over_m123_square > 1 - rmax_square * m12m13_square_p1) && (m23_over_m123 > _m23cut))) {
0121       is_passed = true;
0122     }
0123     return is_passed;
0124   }
0125 
0126   double HEPTopTaggerV2_fixed_R::nsub(
0127       fastjet::PseudoJet jet, int order, fastjet::contrib::Njettiness::AxesMode axes, double beta, double R0) {
0128     fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
0129     return nsub.result(jet);
0130   }
0131 
0132   HEPTopTaggerV2_fixed_R::HEPTopTaggerV2_fixed_R()
0133       : _do_qjets(false),
0134         _mass_drop_threshold(0.8),
0135         _max_subjet_mass(30.),
0136         _mode(Mode(0)),
0137         _mtmass(172.3),
0138         _mwmass(80.4),
0139         _mtmin(150.),
0140         _mtmax(200.),
0141         _rmin(0.85 * 80.4 / 172.3),
0142         _rmax(1.15 * 80.4 / 172.3),
0143         _m23cut(0.35),
0144         _m13cutmin(0.2),
0145         _m13cutmax(1.3),
0146         _minpt_tag(200.),
0147         _nfilt(5),
0148         _Rfilt(0.3),
0149         _jet_algorithm_filter(fastjet::cambridge_algorithm),
0150         _minpt_subjet(0.),
0151         _jet_algorithm_recluster(fastjet::cambridge_algorithm),
0152         _zcut(0.1),
0153         _rcut_factor(0.5),
0154         _q_zcut(0.1),
0155         _q_dcut_fctr(0.5),
0156         _q_exp_min(0.),
0157         _q_exp_max(0.),
0158         _q_rigidity(0.1),
0159         _q_truncation_fctr(0.0),
0160         _rnEngine(nullptr),
0161         //_qjet_plugin(_q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr),
0162         _debug(false),
0163         _first_time(true) {
0164     _djsum = 0.;
0165     _delta_top = 1000000000000.0;
0166     _pruned_mass = 0.;
0167     _unfiltered_mass = 0.;
0168     _top_candidate.reset(0., 0., 0., 0.);
0169     _parts_size = 0;
0170     _is_maybe_top = _is_masscut_passed = _is_ptmincut_passed = false;
0171     _top_subs.clear();
0172     _top_subjets.clear();
0173     _top_hadrons.clear();
0174     _top_parts.clear();
0175     _qweight = -1.;
0176   }
0177 
0178   HEPTopTaggerV2_fixed_R::HEPTopTaggerV2_fixed_R(const fastjet::PseudoJet jet)
0179       : _do_qjets(false),
0180         _jet(jet),
0181         _initial_jet(jet),
0182         _mass_drop_threshold(0.8),
0183         _max_subjet_mass(30.),
0184         _mode(Mode(0)),
0185         _mtmass(172.3),
0186         _mwmass(80.4),
0187         _mtmin(150.),
0188         _mtmax(200.),
0189         _rmin(0.85 * 80.4 / 172.3),
0190         _rmax(1.15 * 80.4 / 172.3),
0191         _m23cut(0.35),
0192         _m13cutmin(0.2),
0193         _m13cutmax(1.3),
0194         _minpt_tag(200.),
0195         _nfilt(5),
0196         _Rfilt(0.3),
0197         _jet_algorithm_filter(fastjet::cambridge_algorithm),
0198         _minpt_subjet(0.),
0199         _jet_algorithm_recluster(fastjet::cambridge_algorithm),
0200         _zcut(0.1),
0201         _rcut_factor(0.5),
0202         _q_zcut(0.1),
0203         _q_dcut_fctr(0.5),
0204         _q_exp_min(0.),
0205         _q_exp_max(0.),
0206         _q_rigidity(0.1),
0207         _q_truncation_fctr(0.0),
0208         _fat(jet),
0209         _rnEngine(nullptr),
0210         _debug(false),
0211         _first_time(true) {}
0212 
0213   HEPTopTaggerV2_fixed_R::HEPTopTaggerV2_fixed_R(const fastjet::PseudoJet jet, double mtmass, double mwmass)
0214       : _do_qjets(false),
0215         _jet(jet),
0216         _initial_jet(jet),
0217         _mass_drop_threshold(0.8),
0218         _max_subjet_mass(30.),
0219         _mode(Mode(0)),
0220         _mtmass(mtmass),
0221         _mwmass(mwmass),
0222         _rmin(0.85 * 80.4 / 172.3),
0223         _rmax(1.15 * 80.4 / 172.3),
0224         _m23cut(0.35),
0225         _m13cutmin(0.2),
0226         _m13cutmax(1.3),
0227         _minpt_tag(200.),
0228         _nfilt(5),
0229         _Rfilt(0.3),
0230         _jet_algorithm_filter(fastjet::cambridge_algorithm),
0231         _minpt_subjet(0.),
0232         _jet_algorithm_recluster(fastjet::cambridge_algorithm),
0233         _zcut(0.1),
0234         _rcut_factor(0.5),
0235         _q_zcut(0.1),
0236         _q_dcut_fctr(0.5),
0237         _q_exp_min(0.),
0238         _q_exp_max(0.),
0239         _q_rigidity(0.1),
0240         _q_truncation_fctr(0.0),
0241         _fat(jet),
0242         _rnEngine(nullptr),
0243         _debug(false),
0244         _first_time(true) {}
0245 
0246   void HEPTopTaggerV2_fixed_R::run() {
0247     print_banner();
0248 
0249     if ((_mode != EARLY_MASSRATIO_SORT_MASS) && (_mode != LATE_MASSRATIO_SORT_MASS) &&
0250         (_mode != EARLY_MASSRATIO_SORT_MODDJADE) && (_mode != LATE_MASSRATIO_SORT_MODDJADE) &&
0251         (_mode != TWO_STEP_FILTER)) {
0252       edm::LogError("HEPTopTaggerV2") << "ERROR: UNKNOWN MODE" << std::endl;
0253       return;
0254     }
0255 
0256     //Qjets
0257     QjetsPlugin _qjet_plugin(_q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr);
0258     _qjet_def = fastjet::JetDefinition(&_qjet_plugin);
0259     _qweight = -1;
0260     vector<fastjet::PseudoJet> _q_constits;
0261     ClusterSequence* _qjet_seq;
0262     PseudoJet _qjet;
0263     if (_do_qjets) {
0264       _q_constits = _initial_jet.associated_cluster_sequence()->constituents(_initial_jet);
0265       _qjet_seq = new ClusterSequence(_q_constits, _qjet_def);
0266       _qjet = sorted_by_pt(_qjet_seq->inclusive_jets())[0];
0267       _qjet_seq->delete_self_when_unused();
0268       const QjetsBaseExtras* ext = dynamic_cast<const QjetsBaseExtras*>(_qjet_seq->extras());
0269       _qweight = ext->weight();
0270       _jet = _qjet;
0271       _fat = _qjet;
0272       _qjet_plugin.SetRNEngine(_rnEngine);
0273     }
0274 
0275     //initialization
0276     _djsum = 0.;
0277     _delta_top = 1000000000000.0;
0278     _pruned_mass = 0.;
0279     _unfiltered_mass = 0.;
0280     _top_candidate.reset(0., 0., 0., 0.);
0281     _parts_size = 0;
0282     _is_maybe_top = _is_masscut_passed = _is_ptmincut_passed = false;
0283     _top_subs.clear();
0284     _top_subjets.clear();
0285     _top_hadrons.clear();
0286     _top_parts.clear();
0287 
0288     //find hard substructures
0289     FindHardSubst(_jet, _top_parts);
0290 
0291     if (_top_parts.size() < 3) {
0292       if (_debug) {
0293         edm::LogInfo("HEPTopTaggerV2") << "< 3 hard substructures " << std::endl;
0294       }
0295       return;  //such events are not interesting
0296     }
0297 
0298     // Sort subjets-after-unclustering by pT.
0299     // Necessary so that two-step-filtering can use the leading-three.
0300     _top_parts = sorted_by_pt(_top_parts);
0301 
0302     // loop over triples
0303     _top_parts = sorted_by_pt(_top_parts);
0304     for (unsigned rr = 0; rr < _top_parts.size(); rr++) {
0305       for (unsigned ll = rr + 1; ll < _top_parts.size(); ll++) {
0306         for (unsigned kk = ll + 1; kk < _top_parts.size(); kk++) {
0307           // two-step filtering
0308           // This means that we only look at the triplet formed by the
0309           // three leading-in-pT subjets-after-unclustering.
0310           if ((_mode == TWO_STEP_FILTER) && rr > 0)
0311             continue;
0312           if ((_mode == TWO_STEP_FILTER) && ll > 1)
0313             continue;
0314           if ((_mode == TWO_STEP_FILTER) && kk > 2)
0315             continue;
0316 
0317           //pick triple
0318           PseudoJet triple = join(_top_parts[rr], _top_parts[ll], _top_parts[kk]);
0319 
0320           //filtering
0321           double filt_top_R = std::min(_Rfilt,
0322                                        0.5 * sqrt(std::min(_top_parts[kk].squared_distance(_top_parts[ll]),
0323                                                            std::min(_top_parts[rr].squared_distance(_top_parts[ll]),
0324                                                                     _top_parts[kk].squared_distance(_top_parts[rr])))));
0325           JetDefinition filtering_def(_jet_algorithm_filter, filt_top_R);
0326           fastjet::Filter filter(filtering_def,
0327                                  fastjet::SelectorNHardest(_nfilt) * fastjet::SelectorPtMin(_minpt_subjet));
0328           PseudoJet topcandidate = filter(triple);
0329 
0330           //mass window cut
0331           if (topcandidate.m() < _mtmin || _mtmax < topcandidate.m())
0332             continue;
0333 
0334           // Sanity cut: can't recluster less than 3 objects into three subjets
0335           if (topcandidate.pieces().size() < 3)
0336             continue;
0337 
0338           // Recluster to 3 subjets and apply mass plane cuts
0339           JetDefinition reclustering(_jet_algorithm_recluster, 3.14);
0340           ClusterSequence* cs_top_sub = new ClusterSequence(topcandidate.constituents(), reclustering);
0341           std::vector<PseudoJet> top_subs = sorted_by_pt(cs_top_sub->exclusive_jets(3));
0342           cs_top_sub->delete_self_when_unused();
0343 
0344           // Require the third subjet to be above the pT threshold
0345           if (top_subs[2].perp() < _minpt_subjet)
0346             continue;
0347 
0348           // Modes with early 2d-massplane cuts
0349           if (_mode == EARLY_MASSRATIO_SORT_MASS && !check_mass_criteria(top_subs)) {
0350             continue;
0351           }
0352           if (_mode == EARLY_MASSRATIO_SORT_MODDJADE && !check_mass_criteria(top_subs)) {
0353             continue;
0354           }
0355 
0356           //is this candidate better than the other? -> update
0357           double deltatop = fabs(topcandidate.m() - _mtmass);
0358           double djsum = djademod(top_subs[0], top_subs[1], topcandidate) +
0359                          djademod(top_subs[0], top_subs[2], topcandidate) +
0360                          djademod(top_subs[1], top_subs[2], topcandidate);
0361           bool better = false;
0362 
0363           // Modes 0 and 1 sort by top mass
0364           if ((_mode == EARLY_MASSRATIO_SORT_MASS) || (_mode == LATE_MASSRATIO_SORT_MASS)) {
0365             if (deltatop < _delta_top)
0366               better = true;
0367           }
0368           // Modes 2 and 3 sort by modified jade distance
0369           else if ((_mode == EARLY_MASSRATIO_SORT_MODDJADE) || (_mode == LATE_MASSRATIO_SORT_MODDJADE)) {
0370             if (djsum > _djsum)
0371               better = true;
0372           }
0373           // Mode 4 is the two-step filtering. No sorting necessary as
0374           // we just look at the triplet of highest pT objects after
0375           // unclustering
0376           else if (_mode == TWO_STEP_FILTER) {
0377             better = true;
0378           } else {
0379             edm::LogError("HEPTopTaggerV2") << "ERROR: UNKNOWN MODE (IN DISTANCE MEASURE SELECTION)" << std::endl;
0380             return;
0381           }
0382 
0383           if (better) {
0384             _djsum = djsum;
0385             _delta_top = deltatop;
0386             _is_maybe_top = true;
0387             _top_candidate = topcandidate;
0388             _top_subs = top_subs;
0389             store_topsubjets(top_subs);
0390             _top_hadrons = topcandidate.constituents();
0391             // Pruning
0392             double _Rprun = _initial_jet.validated_cluster_sequence()->jet_def().R();
0393             JetDefinition jet_def_prune(fastjet::cambridge_algorithm, _Rprun);
0394             fastjet::Pruner pruner(jet_def_prune, _zcut, _rcut_factor);
0395             PseudoJet prunedjet = pruner(triple);
0396             _pruned_mass = prunedjet.m();
0397             _unfiltered_mass = triple.m();
0398 
0399             //are all criteria fulfilled?
0400             _is_masscut_passed = false;
0401             if (check_mass_criteria(top_subs)) {
0402               _is_masscut_passed = true;
0403             }
0404             _is_ptmincut_passed = false;
0405             if (_top_candidate.pt() > _minpt_tag) {
0406               _is_ptmincut_passed = true;
0407             }
0408           }  //end better
0409         }    //end kk
0410       }      //end ll
0411     }        //end rr
0412 
0413     return;
0414   }
0415 
0416   void HEPTopTaggerV2_fixed_R::get_info() const {
0417     edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
0418     edm::LogInfo("HEPTopTaggerV2") << "#                          HEPTopTaggerV2 Result" << std::endl;
0419     edm::LogInfo("HEPTopTaggerV2") << "#" << std::endl;
0420     edm::LogInfo("HEPTopTaggerV2") << "# is top candidate: " << _is_maybe_top << std::endl;
0421     edm::LogInfo("HEPTopTaggerV2") << "# mass plane cuts passed: " << _is_masscut_passed << std::endl;
0422     edm::LogInfo("HEPTopTaggerV2") << "# top candidate mass: " << _top_candidate.m() << std::endl;
0423     edm::LogInfo("HEPTopTaggerV2") << "# top candidate (pt, eta, phi): (" << _top_candidate.perp() << ", "
0424                                    << _top_candidate.eta() << ", " << _top_candidate.phi_std() << ")" << std::endl;
0425     edm::LogInfo("HEPTopTaggerV2") << "# top hadrons: " << _top_hadrons.size() << std::endl;
0426     edm::LogInfo("HEPTopTaggerV2") << "# hard substructures: " << _parts_size << std::endl;
0427     edm::LogInfo("HEPTopTaggerV2") << "# |m - mtop| : " << _delta_top << std::endl;
0428     edm::LogInfo("HEPTopTaggerV2") << "# djsum : " << _djsum << std::endl;
0429     edm::LogInfo("HEPTopTaggerV2") << "# is consistency cut passed: " << _is_ptmincut_passed << std::endl;
0430     edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
0431     return;
0432   }
0433 
0434   void HEPTopTaggerV2_fixed_R::get_setting() const {
0435     edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
0436     edm::LogInfo("HEPTopTaggerV2") << "#                         HEPTopTaggerV2 Settings" << std::endl;
0437     edm::LogInfo("HEPTopTaggerV2") << "#" << std::endl;
0438     edm::LogInfo("HEPTopTaggerV2") << "# mode: " << _mode << " (0 = EARLY_MASSRATIO_SORT_MASS) " << std::endl;
0439     edm::LogInfo("HEPTopTaggerV2") << "#        "
0440                                    << " (1 = LATE_MASSRATIO_SORT_MASS)  " << std::endl;
0441     edm::LogInfo("HEPTopTaggerV2") << "#        "
0442                                    << " (2 = EARLY_MASSRATIO_SORT_MODDJADE)  " << std::endl;
0443     edm::LogInfo("HEPTopTaggerV2") << "#        "
0444                                    << " (3 = LATE_MASSRATIO_SORT_MODDJADE)  " << std::endl;
0445     edm::LogInfo("HEPTopTaggerV2") << "#        "
0446                                    << " (4 = TWO_STEP_FILTER)  " << std::endl;
0447     edm::LogInfo("HEPTopTaggerV2") << "# top mass: " << _mtmass << "    ";
0448     edm::LogInfo("HEPTopTaggerV2") << "W mass: " << _mwmass << std::endl;
0449     edm::LogInfo("HEPTopTaggerV2") << "# top mass window: [" << _mtmin << ", " << _mtmax << "]" << std::endl;
0450     edm::LogInfo("HEPTopTaggerV2") << "# W mass ratio: [" << _rmin << ", " << _rmax << "] (["
0451                                    << _rmin * _mtmass / _mwmass << "%, " << _rmax * _mtmass / _mwmass << "%])"
0452                                    << std::endl;
0453     edm::LogInfo("HEPTopTaggerV2") << "# mass plane cuts: (m23cut, m13min, m13max) = (" << _m23cut << ", " << _m13cutmin
0454                                    << ", " << _m13cutmax << ")" << std::endl;
0455     edm::LogInfo("HEPTopTaggerV2") << "# mass_drop_threshold: " << _mass_drop_threshold << "    ";
0456     edm::LogInfo("HEPTopTaggerV2") << "max_subjet_mass: " << _max_subjet_mass << std::endl;
0457     edm::LogInfo("HEPTopTaggerV2") << "# R_filt: " << _Rfilt << "    ";
0458     edm::LogInfo("HEPTopTaggerV2") << "n_filt: " << _nfilt << std::endl;
0459     edm::LogInfo("HEPTopTaggerV2") << "# minimal subjet pt: " << _minpt_subjet << std::endl;
0460     edm::LogInfo("HEPTopTaggerV2") << "# minimal reconstructed pt: " << _minpt_tag << std::endl;
0461     edm::LogInfo("HEPTopTaggerV2") << "# internal jet algorithms (0 = kt, 1 = C/A, 2 = anti-kt): " << std::endl;
0462     edm::LogInfo("HEPTopTaggerV2") << "#   filtering: " << _jet_algorithm_filter << std::endl;
0463     edm::LogInfo("HEPTopTaggerV2") << "#   reclustering: " << _jet_algorithm_recluster << std::endl;
0464     edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
0465 
0466     return;
0467   }
0468 
0469   //uncluster a fat jet to subjets of given cone size
0470   void HEPTopTaggerV2::UnclusterFatjets(const vector<fastjet::PseudoJet>& big_fatjets,
0471                                         vector<fastjet::PseudoJet>& small_fatjets,
0472                                         const ClusterSequence& cseq,
0473                                         const double small_radius) {
0474     for (unsigned i = 0; i < big_fatjets.size(); i++) {
0475       const PseudoJet& this_jet = big_fatjets[i];
0476       PseudoJet parent1(0, 0, 0, 0), parent2(0, 0, 0, 0);
0477       bool test = cseq.has_parents(this_jet, parent1, parent2);
0478       double dR = 100;
0479 
0480       if (test)
0481         dR = sqrt(parent1.squared_distance(parent2));
0482 
0483       if (!test || dR < small_radius) {
0484         small_fatjets.push_back(this_jet);
0485       } else {
0486         vector<fastjet::PseudoJet> parents;
0487         parents.push_back(parent1);
0488         parents.push_back(parent2);
0489         UnclusterFatjets(parents, small_fatjets, cseq, small_radius);
0490       }
0491     }
0492   }
0493 
0494   HEPTopTaggerV2::HEPTopTaggerV2()
0495       : _do_optimalR(false),
0496         _do_qjets(false),
0497         _mass_drop_threshold(0.8),
0498         _max_subjet_mass(30.),
0499         _mode(Mode(0)),
0500         _mtmass(172.3),
0501         _mwmass(80.4),
0502         _mtmin(150.),
0503         _mtmax(200.),
0504         _rmin(0.85 * 80.4 / 172.3),
0505         _rmax(1.15 * 80.4 / 172.3),
0506         _m23cut(0.35),
0507         _m13cutmin(0.2),
0508         _m13cutmax(1.3),
0509         _minpt_tag(200.),
0510         _nfilt(5),
0511         _Rfilt(0.3),
0512         _jet_algorithm_filter(fastjet::cambridge_algorithm),
0513         _minpt_subjet(0.),
0514         _jet_algorithm_recluster(fastjet::cambridge_algorithm),
0515         _zcut(0.1),
0516         _rcut_factor(0.5),
0517         _max_fatjet_R(1.8),
0518         _min_fatjet_R(0.5),
0519         _step_R(0.1),
0520         _optimalR_threshold(0.2),
0521         _R_filt_optimalR_calc(0.2),
0522         _N_filt_optimalR_calc(10),
0523         _r_min_exp_function(&R_opt_calc_funct),
0524         _optimalR_mmin(150.),
0525         _optimalR_mmax(200.),
0526         _optimalR_fw(0.175),
0527         _R_opt_diff(0.3),
0528         _R_opt_reject_min(false),
0529         _R_filt_optimalR_pass(0.2),
0530         _N_filt_optimalR_pass(5),
0531         _R_filt_optimalR_fail(0.3),
0532         _N_filt_optimalR_fail(3),
0533         _q_zcut(0.1),
0534         _q_dcut_fctr(0.5),
0535         _q_exp_min(0.),
0536         _q_exp_max(0.),
0537         _q_rigidity(0.1),
0538         _q_truncation_fctr(0.0),
0539         _rnEngine(nullptr),
0540         _debug(false) {}
0541 
0542   HEPTopTaggerV2::HEPTopTaggerV2(const fastjet::PseudoJet& jet)
0543       : _do_optimalR(false),
0544         _do_qjets(false),
0545         _jet(jet),
0546         _initial_jet(jet),
0547         _mass_drop_threshold(0.8),
0548         _max_subjet_mass(30.),
0549         _mode(Mode(0)),
0550         _mtmass(172.3),
0551         _mwmass(80.4),
0552         _mtmin(150.),
0553         _mtmax(200.),
0554         _rmin(0.85 * 80.4 / 172.3),
0555         _rmax(1.15 * 80.4 / 172.3),
0556         _m23cut(0.35),
0557         _m13cutmin(0.2),
0558         _m13cutmax(1.3),
0559         _minpt_tag(200.),
0560         _nfilt(5),
0561         _Rfilt(0.3),
0562         _jet_algorithm_filter(fastjet::cambridge_algorithm),
0563         _minpt_subjet(0.),
0564         _jet_algorithm_recluster(fastjet::cambridge_algorithm),
0565         _zcut(0.1),
0566         _rcut_factor(0.5),
0567         _max_fatjet_R(jet.validated_cluster_sequence()->jet_def().R()),
0568         _min_fatjet_R(0.5),
0569         _step_R(0.1),
0570         _optimalR_threshold(0.2),
0571         _R_filt_optimalR_calc(0.2),
0572         _N_filt_optimalR_calc(10),
0573         _r_min_exp_function(&R_opt_calc_funct),
0574         _optimalR_mmin(150.),
0575         _optimalR_mmax(200.),
0576         _optimalR_fw(0.175),
0577         _R_opt_diff(0.3),
0578         _R_opt_reject_min(false),
0579         _R_filt_optimalR_pass(0.2),
0580         _N_filt_optimalR_pass(5),
0581         _R_filt_optimalR_fail(0.3),
0582         _N_filt_optimalR_fail(3),
0583         _q_zcut(0.1),
0584         _q_dcut_fctr(0.5),
0585         _q_exp_min(0.),
0586         _q_exp_max(0.),
0587         _q_rigidity(0.1),
0588         _q_truncation_fctr(0.0),
0589         _fat(jet),
0590         _rnEngine(nullptr),
0591         _debug(false) {}
0592 
0593   HEPTopTaggerV2::HEPTopTaggerV2(const fastjet::PseudoJet& jet, double mtmass, double mwmass)
0594       : _do_optimalR(false),
0595         _do_qjets(false),
0596         _jet(jet),
0597         _initial_jet(jet),
0598         _mass_drop_threshold(0.8),
0599         _max_subjet_mass(30.),
0600         _mode(Mode(0)),
0601         _mtmass(mtmass),
0602         _mwmass(mwmass),
0603         _mtmin(150.),
0604         _mtmax(200.),
0605         _rmin(0.85 * 80.4 / 172.3),
0606         _rmax(1.15 * 80.4 / 172.3),
0607         _m23cut(0.35),
0608         _m13cutmin(0.2),
0609         _m13cutmax(1.3),
0610         _minpt_tag(200.),
0611         _nfilt(5),
0612         _Rfilt(0.3),
0613         _jet_algorithm_filter(fastjet::cambridge_algorithm),
0614         _minpt_subjet(0.),
0615         _jet_algorithm_recluster(fastjet::cambridge_algorithm),
0616         _zcut(0.1),
0617         _rcut_factor(0.5),
0618         _max_fatjet_R(jet.validated_cluster_sequence()->jet_def().R()),
0619         _min_fatjet_R(0.5),
0620         _step_R(0.1),
0621         _optimalR_threshold(0.2),
0622         _R_filt_optimalR_calc(0.2),
0623         _N_filt_optimalR_calc(10),
0624         _r_min_exp_function(&R_opt_calc_funct),
0625         _optimalR_mmin(150.),
0626         _optimalR_mmax(200.),
0627         _optimalR_fw(0.175),
0628         _R_opt_diff(0.3),
0629         _R_opt_reject_min(false),
0630         _R_filt_optimalR_pass(0.2),
0631         _N_filt_optimalR_pass(5),
0632         _R_filt_optimalR_fail(0.3),
0633         _N_filt_optimalR_fail(3),
0634         _q_zcut(0.1),
0635         _q_dcut_fctr(0.5),
0636         _q_exp_min(0.),
0637         _q_exp_max(0.),
0638         _q_rigidity(0.1),
0639         _q_truncation_fctr(0.0),
0640         _fat(jet),
0641         _rnEngine(nullptr),
0642         _debug(false) {}
0643 
0644   void HEPTopTaggerV2::run() {
0645     QjetsPlugin _qjet_plugin(_q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr);
0646     int maxR = int(_max_fatjet_R * 10);
0647     int minR = int(_min_fatjet_R * 10);
0648     int stepR = int(_step_R * 10);
0649     _qweight = -1;
0650     _qjet_plugin.SetRNEngine(_rnEngine);
0651 
0652     if (!_do_optimalR) {
0653       HEPTopTaggerV2_fixed_R htt(_jet);
0654       htt.set_mass_drop_threshold(_mass_drop_threshold);
0655       htt.set_max_subjet_mass(_max_subjet_mass);
0656       htt.set_filtering_n(_nfilt);
0657       htt.set_filtering_R(_Rfilt);
0658       htt.set_filtering_minpt_subjet(_minpt_subjet);
0659       htt.set_filtering_jetalgorithm(_jet_algorithm_filter);
0660       htt.set_reclustering_jetalgorithm(_jet_algorithm_recluster);
0661       htt.set_mode(_mode);
0662       htt.set_mt(_mtmass);
0663       htt.set_mw(_mwmass);
0664       htt.set_top_mass_range(_mtmin, _mtmax);
0665       htt.set_mass_ratio_range(_rmin, _rmax);
0666       htt.set_mass_ratio_cut(_m23cut, _m13cutmin, _m13cutmax);
0667       htt.set_top_minpt(_minpt_tag);
0668       htt.set_pruning_zcut(_zcut);
0669       htt.set_pruning_rcut_factor(_rcut_factor);
0670       htt.set_debug(_debug);
0671       htt.set_qjets(_q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr);
0672       htt.run();
0673 
0674       _HEPTopTaggerV2[maxR] = htt;
0675       _Ropt = maxR;
0676       _qweight = htt.q_weight();
0677       _HEPTopTaggerV2_opt = _HEPTopTaggerV2[_Ropt];
0678     } else {
0679       _qjet_def = fastjet::JetDefinition(&_qjet_plugin);
0680       vector<fastjet::PseudoJet> _q_constits;
0681       ClusterSequence* _qjet_seq;
0682       PseudoJet _qjet;
0683       const ClusterSequence* _seq;
0684       _seq = _initial_jet.validated_cluster_sequence();
0685       if (_do_qjets) {
0686         _q_constits = _initial_jet.associated_cluster_sequence()->constituents(_initial_jet);
0687         _qjet_seq = new ClusterSequence(_q_constits, _qjet_def);
0688         _qjet = sorted_by_pt(_qjet_seq->inclusive_jets())[0];
0689         _qjet_seq->delete_self_when_unused();
0690         const QjetsBaseExtras* ext = dynamic_cast<const QjetsBaseExtras*>(_qjet_seq->extras());
0691         _qweight = ext->weight();
0692         _jet = _qjet;
0693         _seq = _qjet_seq;
0694         _fat = _qjet;
0695       }
0696 
0697       // Do MultiR procedure
0698       vector<fastjet::PseudoJet> big_fatjets;
0699       vector<fastjet::PseudoJet> small_fatjets;
0700 
0701       big_fatjets.push_back(_jet);
0702       _Ropt = 0;
0703 
0704       for (int R = maxR; R >= minR; R -= stepR) {
0705         UnclusterFatjets(big_fatjets, small_fatjets, *_seq, R / 10.);
0706 
0707         if (_debug) {
0708           edm::LogInfo("HEPTopTaggerV2") << "R = " << R << " -> n_small_fatjets = " << small_fatjets.size() << endl;
0709         }
0710 
0711         _n_small_fatjets[R] = small_fatjets.size();
0712 
0713         // We are sorting by pt - so start with a negative dummy
0714         double dummy = -99999;
0715 
0716         for (unsigned i = 0; i < small_fatjets.size(); i++) {
0717           HEPTopTaggerV2_fixed_R htt(small_fatjets[i]);
0718           htt.set_mass_drop_threshold(_mass_drop_threshold);
0719           htt.set_max_subjet_mass(_max_subjet_mass);
0720           htt.set_filtering_n(_nfilt);
0721           htt.set_filtering_R(_Rfilt);
0722           htt.set_filtering_minpt_subjet(_minpt_subjet);
0723           htt.set_filtering_jetalgorithm(_jet_algorithm_filter);
0724           htt.set_reclustering_jetalgorithm(_jet_algorithm_recluster);
0725           htt.set_mode(_mode);
0726           htt.set_mt(_mtmass);
0727           htt.set_mw(_mwmass);
0728           htt.set_top_mass_range(_mtmin, _mtmax);
0729           htt.set_mass_ratio_range(_rmin, _rmax);
0730           htt.set_mass_ratio_cut(_m23cut, _m13cutmin, _m13cutmax);
0731           htt.set_top_minpt(_minpt_tag);
0732           htt.set_pruning_zcut(_zcut);
0733           htt.set_pruning_rcut_factor(_rcut_factor);
0734           htt.set_debug(_debug);
0735           htt.set_qjets(_q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr);
0736 
0737           htt.run();
0738 
0739           if (htt.t().perp() > dummy) {
0740             dummy = htt.t().perp();
0741             _HEPTopTaggerV2[R] = htt;
0742           }
0743         }  //End of loop over small_fatjets
0744 
0745         // Only check if we have not found Ropt yet
0746         if (_Ropt == 0 && R < maxR) {
0747           // If the new mass is OUTSIDE the window ..
0748           if (_HEPTopTaggerV2[R].t().m() < (1 - _optimalR_threshold) * _HEPTopTaggerV2[maxR].t().m())
0749             // .. set _Ropt to the previous mass
0750             _Ropt = R + stepR;
0751         }
0752 
0753         big_fatjets = small_fatjets;
0754         small_fatjets.clear();
0755       }  //End of loop over R
0756 
0757       // if we did not find Ropt in the loop:
0758       if (_Ropt == 0 && _HEPTopTaggerV2[maxR].t().m() > 0) {
0759         // either pick the last value (_R_opt_reject_min=false)
0760         // or leace Ropt at zero (_R_opt_reject_min=true)
0761         if (_R_opt_reject_min == false)
0762           _Ropt = minR;
0763       }
0764 
0765       //for the case that there is no tag at all (< 3 hard substructures)
0766       if (_Ropt == 0 && _HEPTopTaggerV2[maxR].t().m() == 0)
0767         _Ropt = maxR;
0768 
0769       _HEPTopTaggerV2_opt = _HEPTopTaggerV2[_Ropt];
0770 
0771       Filter filter_optimalR_calc(_R_filt_optimalR_calc, SelectorNHardest(_N_filt_optimalR_calc));
0772       _R_opt_calc = _r_min_exp_function(filter_optimalR_calc(_fat).pt());
0773       _pt_for_R_opt_calc = filter_optimalR_calc(_fat).pt();
0774 
0775       Filter filter_optimalR_pass(_R_filt_optimalR_pass, SelectorNHardest(_N_filt_optimalR_pass));
0776       Filter filter_optimalR_fail(_R_filt_optimalR_fail, SelectorNHardest(_N_filt_optimalR_fail));
0777       if (optimalR_type() == 1) {
0778         _filt_fat = filter_optimalR_pass(_fat);
0779       } else {
0780         _filt_fat = filter_optimalR_fail(_fat);
0781       }
0782     }
0783   }
0784 
0785   //optimal_R type
0786   int HEPTopTaggerV2::optimalR_type() {
0787     if (_HEPTopTaggerV2_opt.t().m() < _optimalR_mmin || _HEPTopTaggerV2_opt.t().m() > _optimalR_mmax) {
0788       return 0;
0789     }
0790     if (_HEPTopTaggerV2_opt.f_rec() > _optimalR_fw) {
0791       return 0;
0792     }
0793     if (_Ropt / 10. - _R_opt_calc > _R_opt_diff) {
0794       return 0;
0795     }
0796     return 1;
0797   }
0798 
0799   double HEPTopTaggerV2::nsub_unfiltered(int order,
0800                                          fastjet::contrib::Njettiness::AxesMode axes,
0801                                          double beta,
0802                                          double R0) {
0803     fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
0804     return nsub.result(_fat);
0805   }
0806 
0807   double HEPTopTaggerV2::nsub_filtered(int order, fastjet::contrib::Njettiness::AxesMode axes, double beta, double R0) {
0808     fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
0809     return nsub.result(_filt_fat);
0810   }
0811 
0812   HEPTopTaggerV2::~HEPTopTaggerV2() {}
0813   // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
0814 };  // namespace external
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!