| 1 |
// Original author: Brock Tweedie (JHU)
|
| 2 |
// Ported to CMSSW by: Sal Rappoccio (JHU)
|
| 3 |
// $Id: CATopJetAlgorithm.cc,v 1.1 2011/06/08 17:25:32 tboccali Exp $
|
| 4 |
|
| 5 |
#include "VHbbAnalysis/HbbAnalyzer/interface/CATopJetAlgorithm.h"
|
| 6 |
#include "FWCore/Framework/interface/ESHandle.h"
|
| 7 |
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
|
| 8 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
| 9 |
#include "FWCore/Utilities/interface/Exception.h"
|
| 10 |
#include "DataFormats/Math/interface/deltaPhi.h"
|
| 11 |
#include "DataFormats/Math/interface/angle.h"
|
| 12 |
|
| 13 |
|
| 14 |
using namespace std;
|
| 15 |
using namespace reco;
|
| 16 |
using namespace edm;
|
| 17 |
|
| 18 |
|
| 19 |
// Run the algorithm
|
| 20 |
// ------------------
|
| 21 |
void CATopJetAlgorithm::run( const vector<fastjet::PseudoJet> & cell_particles,
|
| 22 |
vector<CompoundPseudoJet> & hardjetsOutput,
|
| 23 |
edm::EventSetup const & c )
|
| 24 |
{
|
| 25 |
if ( verbose_ ) cout << "Welcome to CATopSubJetAlgorithm::run" << endl;
|
| 26 |
|
| 27 |
// Sum Et of the event
|
| 28 |
double sumEt = 0.;
|
| 29 |
|
| 30 |
//make a list of input objects ordered by ET and calculate sum et
|
| 31 |
// list of fastjet pseudojet constituents
|
| 32 |
for (unsigned i = 0; i < cell_particles.size(); ++i) {
|
| 33 |
sumEt += cell_particles[i].perp();
|
| 34 |
}
|
| 35 |
|
| 36 |
// Determine which bin we are in for et clustering
|
| 37 |
int sumEtBinId = -1;
|
| 38 |
for ( unsigned int i = 0; i < sumEtBins_.size(); ++i ) {
|
| 39 |
if ( sumEt > sumEtBins_[i] ) sumEtBinId = i;
|
| 40 |
}
|
| 41 |
if ( verbose_ ) cout << "Using sumEt = " << sumEt << ", bin = " << sumEtBinId << endl;
|
| 42 |
|
| 43 |
// If the sum et is too low, exit
|
| 44 |
if ( sumEtBinId < 0 ) {
|
| 45 |
return;
|
| 46 |
}
|
| 47 |
|
| 48 |
// empty 4-vector
|
| 49 |
fastjet::PseudoJet blankJetA(0,0,0,0);
|
| 50 |
blankJetA.set_user_index(-1);
|
| 51 |
const fastjet::PseudoJet blankJet = blankJetA;
|
| 52 |
|
| 53 |
// Define adjacency variables which depend on which sumEtBin we are in
|
| 54 |
double deltarcut = deltarBins_[sumEtBinId];
|
| 55 |
double nCellMin = nCellBins_[sumEtBinId];
|
| 56 |
|
| 57 |
if ( verbose_ )cout<<"useAdjacency_ = "<<useAdjacency_<<endl;
|
| 58 |
if ( verbose_ && useAdjacency_==0)cout<<"No Adjacency"<<endl;
|
| 59 |
if ( verbose_ && useAdjacency_==1)cout<<"using deltar adjacency"<<endl;
|
| 60 |
if ( verbose_ && useAdjacency_==2)cout<<"using modified adjacency"<<endl;
|
| 61 |
if ( verbose_ && useAdjacency_==3)cout<<"using calorimeter nearest neighbor based adjacency"<<endl;
|
| 62 |
if ( verbose_ && useAdjacency_==1)cout << "Using deltarcut = " << deltarcut << endl;
|
| 63 |
if ( verbose_ && useAdjacency_==3)cout << "Using nCellMin = " << nCellMin << endl;
|
| 64 |
|
| 65 |
|
| 66 |
// Define strategy, recombination scheme, and jet definition
|
| 67 |
fastjet::Strategy strategy = fastjet::Best;
|
| 68 |
fastjet::RecombinationScheme recombScheme = fastjet::E_scheme;
|
| 69 |
|
| 70 |
// pick algorithm
|
| 71 |
fastjet::JetAlgorithm algorithm = static_cast<fastjet::JetAlgorithm>( algorithm_ );
|
| 72 |
|
| 73 |
fastjet::JetDefinition jetDef( algorithm,
|
| 74 |
rBins_[sumEtBinId], recombScheme, strategy);
|
| 75 |
|
| 76 |
if ( verbose_ ) cout << "About to do jet clustering in CA" << endl;
|
| 77 |
// run the jet clustering
|
| 78 |
fastjet::ClusterSequence clusterSeq(cell_particles, jetDef);
|
| 79 |
|
| 80 |
if ( verbose_ ) cout << "Getting inclusive jets" << endl;
|
| 81 |
// Get the transient inclusive jets
|
| 82 |
// vector<fastjet::PseudoJet> inclusiveJets = clusterSeq.inclusive_jets(ptMin_);
|
| 83 |
vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(clusterSeq.inclusive_jets(ptMin_));
|
| 84 |
|
| 85 |
if ( verbose_ ) cout << "Getting central jets" << endl;
|
| 86 |
// Find the transient central jets
|
| 87 |
vector<fastjet::PseudoJet> centralJets;
|
| 88 |
for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
|
| 89 |
|
| 90 |
// if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].rapidity()) < centralEtaCut_) {
|
| 91 |
if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].eta()) < centralEtaCut_) {
|
| 92 |
centralJets.push_back(inclusiveJets[i]);
|
| 93 |
}
|
| 94 |
}
|
| 95 |
// Sort the transient central jets in Et
|
| 96 |
//// GreaterByEtPseudoJet compEt;
|
| 97 |
//// sort( centralJets.begin(), centralJets.end(), compEt );
|
| 98 |
|
| 99 |
// These will store the 4-vectors of each hard jet
|
| 100 |
vector<math::XYZTLorentzVector> p4_hardJets;
|
| 101 |
|
| 102 |
// These will store the indices of each subjet that
|
| 103 |
// are present in each jet
|
| 104 |
vector<vector<int> > indices( centralJets.size() );
|
| 105 |
|
| 106 |
// Loop over central jets, attempt to find substructure
|
| 107 |
vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(),
|
| 108 |
centralJetsBegin = centralJets.begin(),
|
| 109 |
centralJetsEnd = centralJets.end();
|
| 110 |
if ( verbose_ )cout<<"Loop over jets"<<endl;
|
| 111 |
int i=0;
|
| 112 |
for ( ; jetIt != centralJetsEnd; ++jetIt ) {
|
| 113 |
if ( verbose_ )cout<<"\nJet "<<i<<endl;
|
| 114 |
i++;
|
| 115 |
fastjet::PseudoJet localJet = *jetIt;
|
| 116 |
|
| 117 |
// Get the 4-vector for this jet
|
| 118 |
p4_hardJets.push_back( math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e() ));
|
| 119 |
|
| 120 |
// jet decomposition. try to find 3 or 4 hard, well-localized subjets, characteristic of a boosted top.
|
| 121 |
if ( verbose_ )cout<<"local jet pt = "<<localJet.perp()<<endl;
|
| 122 |
if ( verbose_ )cout<<"deltap = "<<ptFracBins_[sumEtBinId]<<endl;
|
| 123 |
|
| 124 |
double ptHard = ptFracBins_[sumEtBinId]*localJet.perp();
|
| 125 |
vector<fastjet::PseudoJet> leftoversAll;
|
| 126 |
|
| 127 |
// stage 1: primary decomposition. look for when the jet declusters into two hard subjets
|
| 128 |
if ( verbose_ ) cout << "Doing decomposition 1" << endl;
|
| 129 |
fastjet::PseudoJet ja, jb;
|
| 130 |
vector<fastjet::PseudoJet> leftovers1;
|
| 131 |
bool hardBreak1 = decomposeJet(localJet,clusterSeq,cell_particles,ptHard,nCellMin,deltarcut,ja,jb,leftovers1);
|
| 132 |
leftoversAll.insert(leftoversAll.end(),leftovers1.begin(),leftovers1.end());
|
| 133 |
|
| 134 |
/* // stage 2: secondary decomposition. look for when the hard subjets found above further decluster into two hard sub-subjets
|
| 135 |
//
|
| 136 |
// ja -> jaa+jab ?
|
| 137 |
if ( verbose_ ) cout << "Doing decomposition 2. ja->jaa+jab?" << endl;
|
| 138 |
fastjet::PseudoJet jaa, jab;
|
| 139 |
vector<fastjet::PseudoJet> leftovers2a;
|
| 140 |
bool hardBreak2a = false;
|
| 141 |
if (hardBreak1) hardBreak2a = decomposeJet(ja,clusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jaa,jab,leftovers2a);
|
| 142 |
leftoversAll.insert(leftoversAll.end(),leftovers2a.begin(),leftovers2a.end());
|
| 143 |
// jb -> jba+jbb ?
|
| 144 |
if ( verbose_ ) cout << "Doing decomposition 2. ja->jba+jbb?" << endl;
|
| 145 |
fastjet::PseudoJet jba, jbb;
|
| 146 |
vector<fastjet::PseudoJet> leftovers2b;
|
| 147 |
bool hardBreak2b = false;
|
| 148 |
if (hardBreak1) hardBreak2b = decomposeJet(jb,clusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jba,jbb,leftovers2b);
|
| 149 |
leftoversAll.insert(leftoversAll.end(),leftovers2b.begin(),leftovers2b.end());
|
| 150 |
*/
|
| 151 |
//// double DeltaR12=TMath::Sqrt((ja.eta()-jb.eta())*(ja.eta()-jb.eta())+(ja.phi()-jb.phi())*(ja.phi()-jb.phi()));
|
| 152 |
//// double y=TMath::Min(ja.perp2(),jb.perp2())*DeltaR12*DeltaR12/(localJet.m()*localJet.m());
|
| 153 |
|
| 154 |
double DeltaR12=TMath::Sqrt(ja.squared_distance(jb));
|
| 155 |
double y=ja.kt_distance(jb)/localJet.m2();
|
| 156 |
|
| 157 |
if(y>0.09 && hardBreak1){
|
| 158 |
|
| 159 |
fastjet::PseudoJet hardSubA = ja;
|
| 160 |
fastjet::PseudoJet hardSubB = jb;
|
| 161 |
|
| 162 |
double Rfilt=TMath::Min(0.3,DeltaR12/2.);
|
| 163 |
fastjet::JetDefinition jetDef2( algorithm,
|
| 164 |
Rfilt, recombScheme, strategy);
|
| 165 |
|
| 166 |
if ( verbose_ ) cout << "About to do sub jet clustering in CA with R=Rfilt" << endl;
|
| 167 |
// vector<fastjet::PseudoJet> hcand;
|
| 168 |
// hcand.push_back(localJet);
|
| 169 |
vector<fastjet::PseudoJet> hcandConstituents = clusterSeq.constituents( localJet );
|
| 170 |
// run the jet clustering
|
| 171 |
fastjet::ClusterSequence clusterSubSeq(hcandConstituents, jetDef2);
|
| 172 |
|
| 173 |
if ( verbose_ ) cout << "Getting inclusive sub jets" << endl;
|
| 174 |
// Get the transient inclusive jets
|
| 175 |
// vector<fastjet::PseudoJet> inclusiveSubJets = clusterSubSeq.inclusive_jets(ptMin_);
|
| 176 |
vector<fastjet::PseudoJet> inclusiveSubJets = fastjet::sorted_by_pt(clusterSubSeq.inclusive_jets(ptMin_));
|
| 177 |
//// vector<fastjet::PseudoJet> inclusiveSubJets = fastjet::sorted_by_pt(clusterSubSeq.inclusive_jets());
|
| 178 |
|
| 179 |
if ( verbose_ ) cout << "Getting central sub jets" << endl;
|
| 180 |
// Find the transient central jets
|
| 181 |
vector<fastjet::PseudoJet> centralSubJets;
|
| 182 |
for (unsigned int i = 0; i < inclusiveSubJets.size(); i++) {
|
| 183 |
// if (inclusiveSubJets[i].perp() > ptMin_ && fabs(inclusiveSubJets[i].rapidity()) < centralEtaCut_) {
|
| 184 |
if (inclusiveSubJets[i].perp() > ptMin_ && fabs(inclusiveSubJets[i].eta()) < centralEtaCut_) {
|
| 185 |
//v2 if (inclusiveSubJets[i].perp() > ptMin_) {
|
| 186 |
centralSubJets.push_back(inclusiveSubJets[i]);
|
| 187 |
}
|
| 188 |
}
|
| 189 |
|
| 190 |
////// sort(centralSubJets.begin(), centralSubJets.end(), compEt );
|
| 191 |
|
| 192 |
if ( hardSubA.user_index() >= 0 ) centralSubJets.push_back(hardSubA);
|
| 193 |
if ( hardSubB.user_index() >= 0 ) centralSubJets.push_back(hardSubB);
|
| 194 |
|
| 195 |
|
| 196 |
// create the subjets objects to put into the "output" objects
|
| 197 |
vector<CompoundPseudoSubJet> subjetsOutput;
|
| 198 |
std::vector<fastjet::PseudoJet>::const_iterator itSubJetBegin = centralSubJets.begin(),
|
| 199 |
itSubJet = itSubJetBegin, itSubJetEnd = centralSubJets.end();
|
| 200 |
for (; itSubJet != itSubJetEnd; ++itSubJet ){
|
| 201 |
// if ( verbose_ ) cout << "Adding input collection element " << (*itSubJet).user_index() << endl;
|
| 202 |
// if ( (*itSubJet).user_index() >= 0 && (*itSubJet).user_index() < cell_particles.size() )
|
| 203 |
|
| 204 |
// Get the transient subjet constituents from fastjet
|
| 205 |
vector<fastjet::PseudoJet> subjetFastjetConstituents = clusterSeq.constituents( *itSubJet );
|
| 206 |
|
| 207 |
// Get the indices of the constituents
|
| 208 |
vector<int> constituents;
|
| 209 |
|
| 210 |
// Loop over the constituents and get the indices
|
| 211 |
vector<fastjet::PseudoJet>::const_iterator fastSubIt = subjetFastjetConstituents.begin(),
|
| 212 |
transConstBegin = subjetFastjetConstituents.begin(),
|
| 213 |
transConstEnd = subjetFastjetConstituents.end();
|
| 214 |
for ( ; fastSubIt != transConstEnd; ++fastSubIt ) {
|
| 215 |
if ( fastSubIt->user_index() >= 0 && static_cast<unsigned int>(fastSubIt->user_index()) < cell_particles.size() ) {
|
| 216 |
constituents.push_back( fastSubIt->user_index() );
|
| 217 |
}
|
| 218 |
}
|
| 219 |
|
| 220 |
// Make a CompoundPseudoSubJet object to hold this subjet and the indices of its constituents
|
| 221 |
subjetsOutput.push_back( CompoundPseudoSubJet( *itSubJet, constituents ) );
|
| 222 |
}
|
| 223 |
|
| 224 |
|
| 225 |
// Make a CompoundPseudoJet object to hold this hard jet, and the subjets that make it up
|
| 226 |
hardjetsOutput.push_back( CompoundPseudoJet( *jetIt, subjetsOutput ) );
|
| 227 |
|
| 228 |
} // y cut
|
| 229 |
} // loop central jets
|
| 230 |
}
|
| 231 |
|
| 232 |
|
| 233 |
|
| 234 |
|
| 235 |
//-----------------------------------------------------------------------
|
| 236 |
// determine whether two clusters (made of calorimeter towers) are living on "adjacent" cells. if they are, then
|
| 237 |
// we probably shouldn't consider them to be independent objects!
|
| 238 |
//
|
| 239 |
// From Sal: Ignoring genjet case
|
| 240 |
//
|
| 241 |
bool CATopJetAlgorithm::adjacentCells(const fastjet::PseudoJet & jet1, const fastjet::PseudoJet & jet2,
|
| 242 |
const vector<fastjet::PseudoJet> & cell_particles,
|
| 243 |
const fastjet::ClusterSequence & theClusterSequence,
|
| 244 |
double nCellMin ) const {
|
| 245 |
|
| 246 |
|
| 247 |
double eta1 = jet1.rapidity();
|
| 248 |
double phi1 = jet1.phi();
|
| 249 |
double eta2 = jet2.rapidity();
|
| 250 |
double phi2 = jet2.phi();
|
| 251 |
|
| 252 |
double deta = abs(eta2 - eta1) / 0.087;
|
| 253 |
double dphi = fabs( reco::deltaPhi(phi2,phi1) ) / 0.087;
|
| 254 |
|
| 255 |
return ( ( deta + dphi ) <= nCellMin );
|
| 256 |
}
|
| 257 |
|
| 258 |
|
| 259 |
|
| 260 |
//-------------------------------------------------------------------------
|
| 261 |
// attempt to decompose a jet into "hard" subjets, where hardness is set by ptHard
|
| 262 |
//
|
| 263 |
bool CATopJetAlgorithm::decomposeJet(const fastjet::PseudoJet & theJet,
|
| 264 |
const fastjet::ClusterSequence & theClusterSequence,
|
| 265 |
const vector<fastjet::PseudoJet> & cell_particles,
|
| 266 |
double ptHard, double nCellMin, double deltarcut,
|
| 267 |
fastjet::PseudoJet & ja, fastjet::PseudoJet & jb,
|
| 268 |
vector<fastjet::PseudoJet> & leftovers) const {
|
| 269 |
|
| 270 |
bool goodBreak;
|
| 271 |
fastjet::PseudoJet j = theJet;
|
| 272 |
double InputObjectPt = j.perp();
|
| 273 |
if ( verbose_ )cout<<"Input Object Pt = "<<InputObjectPt<<endl;
|
| 274 |
if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
|
| 275 |
leftovers.clear();
|
| 276 |
if ( verbose_ )cout<<"start while loop"<<endl;
|
| 277 |
|
| 278 |
while (1) { // watch out for infinite loop!
|
| 279 |
goodBreak = theClusterSequence.has_parents(j,ja,jb);
|
| 280 |
if (!goodBreak){
|
| 281 |
if ( verbose_ )cout<<"bad break. this is one cell. can't decluster anymore."<<endl;
|
| 282 |
break; // this is one cell, can't decluster anymore
|
| 283 |
}
|
| 284 |
|
| 285 |
if ( verbose_ )cout<<"good break. ja Pt = "<<ja.perp()<<" jb Pt = "<<jb.perp()<<endl;
|
| 286 |
|
| 287 |
/// Adjacency Requirement ///
|
| 288 |
|
| 289 |
// check if clusters are adjacent using a constant deltar adjacency.
|
| 290 |
double clusters_deltar=fabs(ja.eta()-jb.eta())+fabs(deltaPhi(ja.phi(),jb.phi()));
|
| 291 |
|
| 292 |
if ( verbose_ && useAdjacency_ ==1)cout<<"clusters_deltar = "<<clusters_deltar<<endl;
|
| 293 |
if ( verbose_ && useAdjacency_ ==1)cout<<"deltar cut = "<<deltarcut<<endl;
|
| 294 |
|
| 295 |
if ( useAdjacency_==1 && clusters_deltar < deltarcut){
|
| 296 |
if ( verbose_ )cout<<"clusters too close. consant adj. break."<<endl;
|
| 297 |
break;
|
| 298 |
}
|
| 299 |
|
| 300 |
// Check if clusters are adjacent using a DeltaR adjacency which is a function of pT.
|
| 301 |
double clusters_deltaR=deltaR( ja.eta(), ja.phi(), jb.eta(), jb.phi() );
|
| 302 |
|
| 303 |
if ( verbose_ && useAdjacency_ ==2)cout<<"clusters_deltaR = "<<clusters_deltaR<<endl;
|
| 304 |
if ( verbose_ && useAdjacency_ ==2)cout<<"0.4-0.0004*InputObjectPt = "<<0.4-0.0004*InputObjectPt<<endl;
|
| 305 |
|
| 306 |
if ( useAdjacency_==2 && clusters_deltaR < 0.4-0.0004*InputObjectPt)
|
| 307 |
{
|
| 308 |
if ( verbose_ )cout<<"clusters too close. modified adj. break."<<endl;
|
| 309 |
break;
|
| 310 |
}
|
| 311 |
|
| 312 |
// Check if clusters are adjacent in the calorimeter.
|
| 313 |
if ( useAdjacency_==3 && adjacentCells(ja,jb,cell_particles,theClusterSequence,nCellMin) ){
|
| 314 |
if ( verbose_ )cout<<"clusters too close in the calorimeter. calorimeter adj. break."<<endl;
|
| 315 |
break; // the clusters are "adjacent" in the calorimeter => shouldn't have decomposed
|
| 316 |
}
|
| 317 |
|
| 318 |
if ( verbose_ )cout<<"clusters pass distance cut"<<endl;
|
| 319 |
|
| 320 |
/// Pt Fraction Requirement ///
|
| 321 |
|
| 322 |
if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
|
| 323 |
|
| 324 |
if(ja.m() > jb.m() && ja.m()< 0.67*j.m() ) return true;
|
| 325 |
else if(jb.m() > ja.m() && jb.m()< 0.67*j.m() ) return true;
|
| 326 |
|
| 327 |
else if (ja.m() > jb.m()) { // broke into one hard and one soft, ditch the soft one and try again
|
| 328 |
j = ja;
|
| 329 |
vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(jb);
|
| 330 |
leftovers.insert(leftovers.end(),particles.begin(),particles.end());
|
| 331 |
}
|
| 332 |
else {
|
| 333 |
j = jb;
|
| 334 |
vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(ja);
|
| 335 |
leftovers.insert(leftovers.end(),particles.begin(),particles.end());
|
| 336 |
}
|
| 337 |
|
| 338 |
|
| 339 |
}
|
| 340 |
|
| 341 |
if ( verbose_ )cout<<"did not decluster."<<endl; // did not decluster into hard subjets
|
| 342 |
|
| 343 |
ja.reset(0,0,0,0);
|
| 344 |
jb.reset(0,0,0,0);
|
| 345 |
leftovers.clear();
|
| 346 |
return false;
|
| 347 |
}
|
