Reference Manual

Go to the documentation of this file.
00001 // Original author: Brock Tweedie (JHU)
00002 // Ported to CMSSW by: Sal Rappoccio (JHU)
00003 // $Id: CATopJetAlgorithm.cc,v 1.7 2010/03/23 21:02:09 srappocc Exp $
00004 
00005 #include "RecoJets/JetAlgorithms/interface/CATopJetAlgorithm.h"
00006 #include "FWCore/Framework/interface/ESHandle.h"
00007 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00009 #include "FWCore/Utilities/interface/Exception.h"
00010 #include "DataFormats/Math/interface/deltaPhi.h"
00011 #include "DataFormats/Math/interface/angle.h"
00012 
00013 
00014 using namespace std;
00015 using namespace reco;
00016 using namespace edm;
00017 
00018 
00019 //  Run the algorithm
00020 //  ------------------
00021 void CATopJetAlgorithm::run( const vector<fastjet::PseudoJet> & cell_particles, 
00022                                                         vector<CompoundPseudoJet> & hardjetsOutput,
00023                                                         edm::EventSetup const & c )  
00024 {
00025         if ( verbose_ ) cout << "Welcome to CATopSubJetAlgorithm::run" << endl;
00026         
00027         // Sum Et of the event
00028         double sumEt = 0.;
00029         
00030         //make a list of input objects ordered by ET and calculate sum et
00031         // list of fastjet pseudojet constituents
00032         for (unsigned i = 0; i < cell_particles.size(); ++i) {
00033                 sumEt += cell_particles[i].perp();
00034         }
00035         
00036         // Determine which bin we are in for et clustering
00037         int sumEtBinId = -1;
00038         for ( unsigned int i = 0; i < sumEtBins_.size(); ++i ) {
00039                 if ( sumEt > sumEtBins_[i] ) sumEtBinId = i;
00040         }
00041         if ( verbose_ ) cout << "Using sumEt = " << sumEt << ", bin = " << sumEtBinId << endl;
00042         
00043         // If the sum et is too low, exit
00044         if ( sumEtBinId < 0 ) {    
00045                 return;
00046         }
00047         
00048         // empty 4-vector
00049         fastjet::PseudoJet blankJetA(0,0,0,0);
00050         blankJetA.set_user_index(-1);
00051         const fastjet::PseudoJet blankJet = blankJetA;
00052         
00053         // Define adjacency variables which depend on which sumEtBin we are in
00054         double deltarcut = deltarBins_[sumEtBinId];
00055         double nCellMin = nCellBins_[sumEtBinId];
00056         
00057         if ( verbose_ )cout<<"useAdjacency_ = "<<useAdjacency_<<endl;
00058         if ( verbose_ && useAdjacency_==0)cout<<"No Adjacency"<<endl;
00059         if ( verbose_ && useAdjacency_==1)cout<<"using deltar adjacency"<<endl;
00060         if ( verbose_ && useAdjacency_==2)cout<<"using modified adjacency"<<endl;
00061         if ( verbose_ && useAdjacency_==3)cout<<"using calorimeter nearest neighbor based adjacency"<<endl;
00062         if ( verbose_ && useAdjacency_==1)cout << "Using deltarcut = " << deltarcut << endl;
00063         if ( verbose_ && useAdjacency_==3)cout << "Using nCellMin = " << nCellMin << endl;
00064         
00065         
00066         // Define strategy, recombination scheme, and jet definition
00067         fastjet::Strategy strategy = fastjet::Best;
00068         fastjet::RecombinationScheme recombScheme = fastjet::E_scheme;
00069         
00070         // pick algorithm
00071         fastjet::JetAlgorithm algorithm = static_cast<fastjet::JetAlgorithm>( algorithm_ );
00072         
00073         fastjet::JetDefinition jetDef( algorithm, 
00074                                                                   rBins_[sumEtBinId], recombScheme, strategy);
00075         
00076         if ( verbose_ ) cout << "About to do jet clustering in CA" << endl;
00077         // run the jet clustering
00078         fastjet::ClusterSequence clusterSeq(cell_particles, jetDef);
00079         
00080         if ( verbose_ ) cout << "Getting inclusive jets" << endl;
00081         // Get the transient inclusive jets
00082         vector<fastjet::PseudoJet> inclusiveJets = clusterSeq.inclusive_jets(ptMin_);
00083         
00084         if ( verbose_ ) cout << "Getting central jets" << endl;
00085         // Find the transient central jets
00086         vector<fastjet::PseudoJet> centralJets;
00087         for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
00088                 
00089                 if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].rapidity()) < centralEtaCut_) {
00090                         centralJets.push_back(inclusiveJets[i]);
00091                 }
00092         }
00093         // Sort the transient central jets in Et
00094         GreaterByEtPseudoJet compEt;
00095         sort( centralJets.begin(), centralJets.end(), compEt );
00096         
00097         // These will store the 4-vectors of each hard jet
00098         vector<math::XYZTLorentzVector> p4_hardJets;
00099         
00100         // These will store the indices of each subjet that 
00101         // are present in each jet
00102         vector<vector<int> > indices( centralJets.size() );
00103         
00104         // Loop over central jets, attempt to find substructure
00105         vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(),
00106     centralJetsBegin = centralJets.begin(),
00107     centralJetsEnd = centralJets.end();
00108         if ( verbose_ )cout<<"Loop over jets"<<endl;
00109         int i=0;
00110         for ( ; jetIt != centralJetsEnd; ++jetIt ) {
00111                 if ( verbose_ )cout<<"\nJet "<<i<<endl;
00112                 i++;
00113                 fastjet::PseudoJet localJet = *jetIt;
00114                 
00115                 // Get the 4-vector for this jet
00116                 p4_hardJets.push_back( math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e() ));
00117                 
00118                 // jet decomposition.  try to find 3 or 4 hard, well-localized subjets, characteristic of a boosted top.
00119                 if ( verbose_ )cout<<"local jet pt = "<<localJet.perp()<<endl;
00120                 if ( verbose_ )cout<<"deltap = "<<ptFracBins_[sumEtBinId]<<endl;
00121                 
00122                 double ptHard = ptFracBins_[sumEtBinId]*localJet.perp();
00123                 vector<fastjet::PseudoJet> leftoversAll;
00124                 
00125                 // stage 1:  primary decomposition.  look for when the jet declusters into two hard subjets
00126                 if ( verbose_ ) cout << "Doing decomposition 1" << endl;
00127                 fastjet::PseudoJet ja, jb;
00128                 vector<fastjet::PseudoJet> leftovers1;
00129                 bool hardBreak1 = decomposeJet(localJet,clusterSeq,cell_particles,ptHard,nCellMin,deltarcut,ja,jb,leftovers1);
00130                 leftoversAll.insert(leftoversAll.end(),leftovers1.begin(),leftovers1.end());
00131                 
00132                 // stage 2:  secondary decomposition.  look for when the hard subjets found above further decluster into two hard sub-subjets
00133                 //
00134                 // ja -> jaa+jab ?
00135                 if ( verbose_ ) cout << "Doing decomposition 2. ja->jaa+jab?" << endl;
00136                 fastjet::PseudoJet jaa, jab;
00137                 vector<fastjet::PseudoJet> leftovers2a;
00138                 bool hardBreak2a = false;
00139                 if (hardBreak1)  hardBreak2a = decomposeJet(ja,clusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jaa,jab,leftovers2a);
00140                 leftoversAll.insert(leftoversAll.end(),leftovers2a.begin(),leftovers2a.end());
00141                 // jb -> jba+jbb ?
00142                 if ( verbose_ ) cout << "Doing decomposition 2. ja->jba+jbb?" << endl;
00143                 fastjet::PseudoJet jba, jbb;
00144                 vector<fastjet::PseudoJet> leftovers2b;
00145                 bool hardBreak2b = false;
00146                 if (hardBreak1)  hardBreak2b = decomposeJet(jb,clusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jba,jbb,leftovers2b);
00147                 leftoversAll.insert(leftoversAll.end(),leftovers2b.begin(),leftovers2b.end());
00148                 
00149                 // NOTE:  it might be good to consider some checks for whether these subjets can be further decomposed.  e.g., the above procedure leaves
00150                 //        open the possibility of "subjets" that actually consist of two or more distinct hard clusters.  however, this kind of thing
00151                 //        is a rarity for the simulations so far considered.
00152                 
00153                 // proceed if one or both of the above hard subjets successfully decomposed
00154                 if ( verbose_ ) cout << "Done with decomposition" << endl;
00155                 
00156                 int nBreak2 = 0;
00157                 fastjet::PseudoJet hardA = blankJet, hardB = blankJet, hardC = blankJet, hardD = blankJet;
00158                 if (!hardBreak2a && !hardBreak2b) { nBreak2 = 0; hardA = ja;  hardB = jb;  hardC = blankJet; hardD = blankJet; }
00159                 if ( hardBreak2a && !hardBreak2b) { nBreak2 = 1; hardA = jaa; hardB = jab; hardC = jb;       hardD = blankJet; }
00160                 if (!hardBreak2a &&  hardBreak2b) { nBreak2 = 1; hardA = jba; hardB = jbb; hardC = ja;       hardD = blankJet;}
00161                 if ( hardBreak2a &&  hardBreak2b) { nBreak2 = 2; hardA = jaa; hardB = jab; hardC = jba;      hardD = jbb; }
00162                 
00163                 // check if we are left with >= 3 hard subjets
00164                 fastjet::PseudoJet subjet1 = blankJet;
00165                 fastjet::PseudoJet subjet2 = blankJet;
00166                 fastjet::PseudoJet subjet3 = blankJet;
00167                 fastjet::PseudoJet subjet4 = blankJet;
00168                 subjet1 = hardA; subjet2 = hardB; subjet3 = hardC; subjet4 = hardD;
00169                 
00170                 // record the hard subjets
00171                 vector<fastjet::PseudoJet> hardSubjets;
00172                 
00173                 if ( subjet1.user_index() >= 0 )
00174                         hardSubjets.push_back(subjet1);
00175                 if ( subjet2.user_index() >= 0 )
00176                         hardSubjets.push_back(subjet2);
00177                 if ( subjet3.user_index() >= 0 )
00178                         hardSubjets.push_back(subjet3);
00179                 if ( subjet4.user_index() >= 0 )
00180                         hardSubjets.push_back(subjet4);
00181                 sort(hardSubjets.begin(), hardSubjets.end(), compEt );
00182                 
00183                 // create the subjets objects to put into the "output" objects
00184                 vector<CompoundPseudoSubJet>  subjetsOutput;
00185                 std::vector<fastjet::PseudoJet>::const_iterator itSubJetBegin = hardSubjets.begin(),
00186                 itSubJet = itSubJetBegin, itSubJetEnd = hardSubjets.end();
00187                 for (; itSubJet != itSubJetEnd; ++itSubJet ){
00188                         //       if ( verbose_ ) cout << "Adding input collection element " << (*itSubJet).user_index() << endl;
00189                         //       if ( (*itSubJet).user_index() >= 0 && (*itSubJet).user_index() < cell_particles.size() )
00190                         
00191                         // Get the transient subjet constituents from fastjet
00192                         vector<fastjet::PseudoJet> subjetFastjetConstituents = clusterSeq.constituents( *itSubJet );
00193                         
00194                         // Get the indices of the constituents
00195                         vector<int> constituents;
00196                         
00197                         // Loop over the constituents and get the indices
00198                         vector<fastjet::PseudoJet>::const_iterator fastSubIt = subjetFastjetConstituents.begin(),
00199                         transConstBegin = subjetFastjetConstituents.begin(),
00200                         transConstEnd = subjetFastjetConstituents.end();
00201                         for ( ; fastSubIt != transConstEnd; ++fastSubIt ) {
00202                                 if ( fastSubIt->user_index() >= 0 && static_cast<unsigned int>(fastSubIt->user_index()) < cell_particles.size() ) {
00203                                         constituents.push_back( fastSubIt->user_index() );
00204                                 }
00205                         }
00206                         
00207                         // Make a CompoundPseudoSubJet object to hold this subjet and the indices of its constituents
00208                         subjetsOutput.push_back( CompoundPseudoSubJet( *itSubJet, constituents ) );
00209                 }
00210                 
00211                 
00212                 // Make a CompoundPseudoJet object to hold this hard jet, and the subjets that make it up
00213                 hardjetsOutput.push_back( CompoundPseudoJet( *jetIt, subjetsOutput ) );
00214                 
00215         }
00216 }
00217 
00218 
00219 
00220 
00221 //-----------------------------------------------------------------------
00222 // determine whether two clusters (made of calorimeter towers) are living on "adjacent" cells.  if they are, then
00223 // we probably shouldn't consider them to be independent objects!
00224 //
00225 // From Sal: Ignoring genjet case
00226 //
00227 bool CATopJetAlgorithm::adjacentCells(const fastjet::PseudoJet & jet1, const fastjet::PseudoJet & jet2, 
00228                                                                           const vector<fastjet::PseudoJet> & cell_particles,
00229                                                                           const fastjet::ClusterSequence & theClusterSequence,
00230                                                                           double nCellMin ) const {
00231         
00232         
00233         double eta1 = jet1.rapidity();
00234         double phi1 = jet1.phi();
00235         double eta2 = jet2.rapidity();
00236         double phi2 = jet2.phi();
00237         
00238         double deta = abs(eta2 - eta1) / 0.087;
00239         double dphi = fabs( reco::deltaPhi(phi2,phi1) ) / 0.087;
00240         
00241         return ( ( deta + dphi ) <= nCellMin );
00242 }
00243 
00244 
00245 
00246 //-------------------------------------------------------------------------
00247 // attempt to decompose a jet into "hard" subjets, where hardness is set by ptHard
00248 //
00249 bool CATopJetAlgorithm::decomposeJet(const fastjet::PseudoJet & theJet, 
00250                                                                          const fastjet::ClusterSequence & theClusterSequence, 
00251                                                                          const vector<fastjet::PseudoJet> & cell_particles,
00252                                                                          double ptHard, double nCellMin, double deltarcut,
00253                                                                          fastjet::PseudoJet & ja, fastjet::PseudoJet & jb, 
00254                                                                          vector<fastjet::PseudoJet> & leftovers) const {
00255         
00256         bool goodBreak;
00257         fastjet::PseudoJet j = theJet;
00258         double InputObjectPt = j.perp();
00259         if ( verbose_ )cout<<"Input Object Pt = "<<InputObjectPt<<endl;
00260         if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
00261         leftovers.clear();
00262         if ( verbose_ )cout<<"start while loop"<<endl;
00263         
00264         while (1) {                                                      // watch out for infinite loop!
00265                 goodBreak = theClusterSequence.has_parents(j,ja,jb);
00266                 if (!goodBreak){
00267                         if ( verbose_ )cout<<"bad break. this is one cell. can't decluster anymore."<<endl;
00268                         break;         // this is one cell, can't decluster anymore
00269                 }
00270                 
00271                 if ( verbose_ )cout<<"good break. ja Pt = "<<ja.perp()<<" jb Pt = "<<jb.perp()<<endl;
00272                 
00274                 
00275                 // check if clusters are adjacent using a constant deltar adjacency.
00276                 double clusters_deltar=fabs(ja.eta()-jb.eta())+fabs(deltaPhi(ja.phi(),jb.phi()));
00277                 
00278                 if ( verbose_  && useAdjacency_ ==1)cout<<"clusters_deltar = "<<clusters_deltar<<endl;
00279                 if ( verbose_  && useAdjacency_ ==1)cout<<"deltar cut = "<<deltarcut<<endl;
00280                 
00281                 if ( useAdjacency_==1 && clusters_deltar < deltarcut){
00282                         if ( verbose_ )cout<<"clusters too close. consant adj. break."<<endl;
00283                         break;
00284                 } 
00285                 
00286                 // Check if clusters are adjacent using a DeltaR adjacency which is a function of pT.
00287                 double clusters_deltaR=deltaR( ja.eta(), ja.phi(), jb.eta(), jb.phi() );
00288                 
00289                 if ( verbose_  && useAdjacency_ ==2)cout<<"clusters_deltaR = "<<clusters_deltaR<<endl;
00290                 if ( verbose_  && useAdjacency_ ==2)cout<<"0.4-0.0004*InputObjectPt = "<<0.4-0.0004*InputObjectPt<<endl;
00291                 
00292                 if ( useAdjacency_==2 && clusters_deltaR < 0.4-0.0004*InputObjectPt)
00293                 {
00294                         if ( verbose_ )cout<<"clusters too close. modified adj. break."<<endl;
00295                         break;
00296                 } 
00297 
00298                 // Check if clusters are adjacent in the calorimeter. 
00299                 if ( useAdjacency_==3 &&  adjacentCells(ja,jb,cell_particles,theClusterSequence,nCellMin) ){                  
00300                         if ( verbose_ )cout<<"clusters too close in the calorimeter. calorimeter adj. break."<<endl;
00301                         break;         // the clusters are "adjacent" in the calorimeter => shouldn't have decomposed
00302                 }
00303                                 
00304                 if ( verbose_ )cout<<"clusters pass distance cut"<<endl;
00305                 
00307                 
00308                 if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
00309                 
00310                 if (ja.perp() < ptHard && jb.perp() < ptHard){
00311                         if ( verbose_ )cout<<"two soft clusters. dead end"<<endl;
00312                         break;         // broke into two soft clusters, dead end
00313                 }
00314                 
00315                 if (ja.perp() > ptHard && jb.perp() > ptHard){
00316                         if ( verbose_ )cout<<"two hard clusters. done"<<endl;
00317                         return true;   // broke into two hard clusters, we're done!
00318                 }
00319                 
00320                 else if (ja.perp() > jb.perp()) {                              // broke into one hard and one soft, ditch the soft one and try again
00321                         if ( verbose_ )cout<<"ja hard jb soft. try to split hard. j = ja"<<endl; 
00322                         j = ja;
00323                         vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(jb);
00324                         leftovers.insert(leftovers.end(),particles.begin(),particles.end());
00325                 }
00326                 else {
00327                         if ( verbose_ )cout<<"ja hard jb soft. try to split hard. j = jb"<<endl; 
00328                         j = jb;
00329                         vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(ja);
00330                         leftovers.insert(leftovers.end(),particles.begin(),particles.end());
00331                 }
00332         }
00333         
00334         if ( verbose_ )cout<<"did not decluster."<<endl;  // did not decluster into hard subjets
00335         
00336         ja.reset(0,0,0,0);
00337         jb.reset(0,0,0,0);
00338         leftovers.clear();
00339         return false;
00340 }
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoJets/JetAlgorithms/src/CATopJetAlgorithm.cc
