File indexing completed on 2024-12-20 03:14:04
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include "RecoJets/JetAlgorithms/interface/CMSBoostedTauSeedingAlgorithm.h"
0024
0025 FASTJET_BEGIN_NAMESPACE
0026
0027 namespace contrib {
0028
0029
0030 CMSBoostedTauSeedingAlgorithm::CMSBoostedTauSeedingAlgorithm(double iminPt,
0031 double iminMassDrop,
0032 double imaxMassDrop,
0033 double iminY,
0034 double imaxY,
0035 double iminDeltaR,
0036 double imaxDeltaR,
0037 int maxDepth,
0038 int verbosity)
0039 : ptMin_(iminPt),
0040 muMin_(iminMassDrop),
0041 muMax_(imaxMassDrop),
0042 yMin_(iminY),
0043 yMax_(imaxY),
0044 dRMin_(iminDeltaR),
0045 dRMax_(imaxDeltaR),
0046 maxDepth_(maxDepth),
0047 verbosity_(verbosity) {}
0048
0049
0050
0051 std::string CMSBoostedTauSeedingAlgorithm::description() const {
0052 std::ostringstream oss;
0053 oss << "CMSBoostedTauSeedingAlgorithm algorithm";
0054 return oss.str();
0055 }
0056
0057
0058 void CMSBoostedTauSeedingAlgorithm::dumpSubJetStructure(
0059 const fastjet::PseudoJet& jet, int depth, int maxDepth, const std::string& depth_and_idx_string) const {
0060 if (maxDepth != -1 && depth > maxDepth)
0061 return;
0062 fastjet::PseudoJet subjet1, subjet2;
0063 bool hasSubjets = jet.has_parents(subjet1, subjet2);
0064 if (!hasSubjets)
0065 return;
0066 std::string depth_and_idx_string_subjet1 = depth_and_idx_string;
0067 if (!depth_and_idx_string_subjet1.empty())
0068 depth_and_idx_string_subjet1.append(".");
0069 depth_and_idx_string_subjet1.append("0");
0070 for (int iSpace = 0; iSpace < depth; ++iSpace) {
0071 std::cout << " ";
0072 }
0073 std::cout << " jetConstituent #" << depth_and_idx_string_subjet1 << " (depth = " << depth
0074 << "): Pt = " << subjet1.pt() << ","
0075 << " eta = " << subjet1.eta() << ", phi = " << subjet1.phi() << ", mass = " << subjet1.m()
0076 << " (constituents = " << subjet1.constituents().size() << ")" << std::endl;
0077 dumpSubJetStructure(subjet1, depth + 1, maxDepth, depth_and_idx_string_subjet1);
0078 std::string depth_and_idx_string_subjet2 = depth_and_idx_string;
0079 if (!depth_and_idx_string_subjet2.empty())
0080 depth_and_idx_string_subjet2.append(".");
0081 depth_and_idx_string_subjet2.append("1");
0082 for (int iSpace = 0; iSpace < depth; ++iSpace) {
0083 std::cout << " ";
0084 }
0085 std::cout << " jetConstituent #" << depth_and_idx_string_subjet2 << " (depth = " << depth
0086 << "): Pt = " << subjet2.pt() << ","
0087 << " eta = " << subjet2.eta() << ", phi = " << subjet2.phi() << ", mass = " << subjet2.m()
0088 << " (constituents = " << subjet2.constituents().size() << ")" << std::endl;
0089 dumpSubJetStructure(subjet2, depth + 1, maxDepth, depth_and_idx_string_subjet2);
0090 for (int iSpace = 0; iSpace < depth; ++iSpace) {
0091 std::cout << " ";
0092 }
0093 double dR = subjet1.delta_R(subjet2);
0094 std::cout << " (mass-drop @ " << depth_and_idx_string << " = " << std::max(subjet1.m(), subjet2.m()) / jet.m()
0095 << ", dR = " << dR << ")" << std::endl;
0096 }
0097
0098 std::pair<PseudoJet, PseudoJet> CMSBoostedTauSeedingAlgorithm::findSubjets(
0099 const PseudoJet& jet, int depth, bool& subjetsFound) const {
0100 const float etaMax_ = 3.0;
0101 if (verbosity_ >= 2) {
0102 std::cout << "<CMSBoostedTauSeedingAlgorithm::findSubjets>:" << std::endl;
0103 std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0104 << ", mass = " << jet.m() << std::endl;
0105 }
0106
0107 PseudoJet subjet1, subjet2;
0108 bool hasSubjets = jet.has_parents(subjet1, subjet2);
0109 if (hasSubjets && (maxDepth_ == -1 || depth <= maxDepth_)) {
0110
0111 if (subjet1.m2() < subjet2.m2()) {
0112 std::swap(subjet1, subjet2);
0113 }
0114 double dR = subjet1.delta_R(subjet2);
0115 double kT = subjet1.kt_distance(subjet2);
0116 double mu = (jet.m() > 0.) ? sqrt(std::max(subjet1.m2(), subjet2.m2()) / jet.m2()) : -1.;
0117
0118 if (subjet1.pt() > ptMin_ && fabs(subjet1.eta()) < etaMax_ && subjet2.pt() > ptMin_ &&
0119 fabs(subjet2.eta()) < etaMax_ && dR > dRMin_ && dR < dRMax_ && mu > muMin_ && mu < muMax_ &&
0120 kT < (yMax_ * jet.m2()) && kT > (yMin_ * jet.m2())) {
0121 subjetsFound = true;
0122 return std::make_pair(subjet1, subjet2);
0123 } else if (subjet1.pt() > ptMin_) {
0124 return findSubjets(subjet1, depth + 1, subjetsFound);
0125 } else if (subjet2.pt() > ptMin_) {
0126 return findSubjets(subjet2, depth + 1, subjetsFound);
0127 }
0128 }
0129 subjetsFound = false;
0130 PseudoJet dummy_subjet1, dummy_subjet2;
0131 return std::make_pair(dummy_subjet1, dummy_subjet2);
0132 }
0133
0134 PseudoJet CMSBoostedTauSeedingAlgorithm::result(const PseudoJet& jet) const {
0135 if (verbosity_ >= 1) {
0136 std::cout << "<CMSBoostedTauSeedingAlgorithm::findSubjets>:" << std::endl;
0137 std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0138 << ", mass = " << jet.m() << std::endl;
0139 }
0140
0141 if (verbosity_ >= 2) {
0142 dumpSubJetStructure(jet, 0, maxDepth_, "");
0143 }
0144
0145 bool subjetsFound = false;
0146 std::pair<PseudoJet, PseudoJet> subjets = findSubjets(jet, 0, subjetsFound);
0147 if (subjetsFound) {
0148
0149 PseudoJet subjet1 = subjets.first;
0150 PseudoJet subjet2 = subjets.second;
0151 if (verbosity_ >= 1) {
0152 std::cout << "before recombination:" << std::endl;
0153 std::cout << " subjet #1: Pt = " << subjet1.pt() << ", eta = " << subjet1.eta() << ", phi = " << subjet1.phi()
0154 << ", mass = " << subjet1.m() << std::endl;
0155 std::cout << " subjet #2: Pt = " << subjet2.pt() << ", eta = " << subjet2.eta() << ", phi = " << subjet2.phi()
0156 << ", mass = " << subjet2.m() << std::endl;
0157 }
0158
0159 const JetDefinition::Recombiner* rec = jet.associated_cluster_sequence()->jet_def().recombiner();
0160 PseudoJet result_local = join(subjet1, subjet2, *rec);
0161 if (verbosity_ >= 1) {
0162 std::cout << "after recombination:" << std::endl;
0163 std::vector<fastjet::PseudoJet> subjets = result_local.pieces();
0164 int idx_subjet = 0;
0165 for (std::vector<fastjet::PseudoJet>::const_iterator subjet = subjets.begin(); subjet != subjets.end();
0166 ++subjet) {
0167 std::cout << " subjet #" << idx_subjet << ": Pt = " << subjet->pt() << ", eta = " << subjet->eta()
0168 << ", phi = " << subjet->phi() << ", mass = " << subjet->m()
0169 << " (#constituents = " << subjet->constituents().size() << ")" << std::endl;
0170 std::vector<fastjet::PseudoJet> constituents = subjet->constituents();
0171 int idx_constituent = 0;
0172 for (std::vector<fastjet::PseudoJet>::const_iterator constituent = constituents.begin();
0173 constituent != constituents.end();
0174 ++constituent) {
0175 if (constituent->pt() < 1.e-3)
0176 continue;
0177 std::cout << " constituent #" << idx_constituent << ": Pt = " << constituent->pt()
0178 << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << ","
0179 << " mass = " << constituent->m() << std::endl;
0180 ++idx_constituent;
0181 }
0182 ++idx_subjet;
0183 }
0184 }
0185
0186 CMSBoostedTauSeedingAlgorithmStructure* s = new CMSBoostedTauSeedingAlgorithmStructure(result_local);
0187
0188 s->_mu = (jet.m2() > 0.) ? sqrt(std::max(subjet1.m2(), subjet2.m2()) / jet.m2()) : 0.;
0189 s->_y = (jet.m2() > 0.) ? subjet1.kt_distance(subjet2) / jet.m2() : 0.;
0190 s->_dR = subjet1.delta_R(subjet2);
0191 s->_pt = subjet2.pt();
0192
0193 result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
0194
0195 return result_local;
0196 } else {
0197
0198 if (verbosity_ >= 1) {
0199 std::cout << "No subjets found." << std::endl;
0200 }
0201 return PseudoJet();
0202 }
0203 }
0204
0205 }
0206
0207 FASTJET_END_NAMESPACE