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.13 2012/02/14 19:43:30 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 using namespace std;
00014 using namespace reco;
00015 using namespace edm;
00016 
00017 
00018 //  Run the algorithm
00019 //  ------------------
00020 void CATopJetAlgorithm::run( const vector<fastjet::PseudoJet> & cell_particles, 
00021                              vector<fastjet::PseudoJet> & hardjetsOutput,
00022                              boost::shared_ptr<fastjet::ClusterSequence> & fjClusterSeq
00023                              )  
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         if ( verbose_ ) cout << "About to do jet clustering in CA" << endl;
00066         // run the jet clustering
00067 
00068         //cluster the jets with the jet definition jetDef:
00069         // run algorithm
00070         // boost::shared_ptr<fastjet::ClusterSequence> fjClusterSeq;
00071         // if ( !doAreaFastjet_ ) {
00072         //   fjClusterSeq = boost::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequence( cell_particles, jetDef ) );
00073         // } else if (voronoiRfact_ <= 0) {
00074         //   fjClusterSeq = boost::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceArea( cell_particles, jetDef , *fjActiveArea_ ) );
00075         // } else {
00076         //   fjClusterSeq = boost::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceVoronoiArea( cell_particles, jetDef , fastjet::VoronoiAreaSpec(voronoiRfact_) ) );
00077         // }
00078         
00079         if ( verbose_ ) cout << "Getting inclusive jets" << endl;
00080         // Get the transient inclusive jets
00081         vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq->inclusive_jets(ptMin_);
00082         
00083         if ( verbose_ ) cout << "Getting central jets" << endl;
00084         // Find the transient central jets
00085         vector<fastjet::PseudoJet> centralJets;
00086         for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
00087                 
00088                 if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].rapidity()) < centralEtaCut_) {
00089                         centralJets.push_back(inclusiveJets[i]);
00090                 }
00091         }
00092         // Sort the transient central jets in Et
00093         GreaterByEtPseudoJet compEt;
00094         sort( centralJets.begin(), centralJets.end(), compEt );
00095         
00096         // These will store the 4-vectors of each hard jet
00097         vector<math::XYZTLorentzVector> p4_hardJets;
00098         
00099         // These will store the indices of each subjet that 
00100         // are present in each jet
00101         vector<vector<int> > indices( centralJets.size() );
00102         
00103         // Loop over central jets, attempt to find substructure
00104         vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(),
00105           centralJetsEnd = centralJets.end();
00106         if ( verbose_ )cout<<"Loop over jets"<<endl;
00107         int i=0;
00108         for ( ; jetIt != centralJetsEnd; ++jetIt ) {
00109                 if ( verbose_ )cout<<"\nJet "<<i<<endl;
00110                 i++;
00111                 fastjet::PseudoJet localJet = *jetIt;
00112                 
00113                 // Get the 4-vector for this jet
00114                 p4_hardJets.push_back( math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e() ));
00115                 
00116                 // jet decomposition.  try to find 3 or 4 hard, well-localized subjets, characteristic of a boosted top.
00117                 if ( verbose_ )cout<<"local jet pt = "<<localJet.perp()<<endl;
00118                 if ( verbose_ )cout<<"deltap = "<<ptFracBins_[sumEtBinId]<<endl;
00119                 
00120                 double ptHard = ptFracBins_[sumEtBinId]*localJet.perp();
00121                 vector<fastjet::PseudoJet> leftoversAll;
00122                 
00123                 // stage 1:  primary decomposition.  look for when the jet declusters into two hard subjets
00124                 if ( verbose_ ) cout << "Doing decomposition 1" << endl;
00125                 fastjet::PseudoJet ja, jb;
00126                 vector<fastjet::PseudoJet> leftovers1;
00127                 bool hardBreak1 = decomposeJet(localJet,*fjClusterSeq,cell_particles,ptHard,nCellMin,deltarcut,ja,jb,leftovers1);
00128                 leftoversAll.insert(leftoversAll.end(),leftovers1.begin(),leftovers1.end());
00129                 
00130                 // stage 2:  secondary decomposition.  look for when the hard subjets found above further decluster into two hard sub-subjets
00131                 //
00132                 // ja -> jaa+jab ?
00133                 if ( verbose_ ) cout << "Doing decomposition 2. ja->jaa+jab?" << endl;
00134                 fastjet::PseudoJet jaa, jab;
00135                 vector<fastjet::PseudoJet> leftovers2a;
00136                 bool hardBreak2a = false;
00137                 if (hardBreak1)  hardBreak2a = decomposeJet(ja,*fjClusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jaa,jab,leftovers2a);
00138                 leftoversAll.insert(leftoversAll.end(),leftovers2a.begin(),leftovers2a.end());
00139                 // jb -> jba+jbb ?
00140                 if ( verbose_ ) cout << "Doing decomposition 2. ja->jba+jbb?" << endl;
00141                 fastjet::PseudoJet jba, jbb;
00142                 vector<fastjet::PseudoJet> leftovers2b;
00143                 bool hardBreak2b = false;
00144                 if (hardBreak1)  hardBreak2b = decomposeJet(jb,*fjClusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jba,jbb,leftovers2b);
00145                 leftoversAll.insert(leftoversAll.end(),leftovers2b.begin(),leftovers2b.end());
00146                 
00147                 // NOTE:  it might be good to consider some checks for whether these subjets can be further decomposed.  e.g., the above procedure leaves
00148                 //        open the possibility of "subjets" that actually consist of two or more distinct hard clusters.  however, this kind of thing
00149                 //        is a rarity for the simulations so far considered.
00150                 
00151                 // proceed if one or both of the above hard subjets successfully decomposed
00152                 if ( verbose_ ) cout << "Done with decomposition" << endl;
00153 
00154                 if ( verbose_ ) cout<<"hardBreak1 = "<<hardBreak1<<endl;
00155                 if ( verbose_ ) cout<<"hardBreak2a = "<<hardBreak2a<<endl;
00156                 if ( verbose_ ) cout<<"hardBreak2b = "<<hardBreak2b<<endl;
00157             
00158                 fastjet::PseudoJet hardA = blankJet, hardB = blankJet, hardC = blankJet, hardD = blankJet;
00159                 if (!hardBreak1) { 
00160                   hardA = localJet;  
00161                   hardB = blankJet;  
00162                   hardC = blankJet; 
00163                   hardD = blankJet; 
00164                   if(verbose_)cout<<"Hardbreak failed. Save subjet1=localJet"<<endl;
00165                 } 
00166                 if (hardBreak1 && !hardBreak2a && !hardBreak2b) { 
00167                   hardA = ja;  
00168                   hardB = jb;  
00169                   hardC = blankJet; 
00170                   hardD = blankJet; 
00171                   if(verbose_)cout<<"First decomposition succeeded, both second decompositions failed. Save subjet1=ja subjet2=jb"<<endl;
00172                 }
00173                 if (hardBreak1 && hardBreak2a && !hardBreak2b) { 
00174                   hardA = jaa; 
00175                   hardB = jab; 
00176                   hardC = jb; 
00177                   hardD = blankJet; 
00178                   if(verbose_)cout<<"First decomposition succeeded, ja split succesfully, jb did not split. Save subjet1=jaa subjet2=jab subjet3=jb"<<endl;
00179                 }
00180                 if (hardBreak1 && !hardBreak2a &&  hardBreak2b) { 
00181                   hardA = jba; 
00182                   hardB = jbb; 
00183                   hardC = ja; 
00184                   hardD = blankJet; 
00185                   if(verbose_)cout<<"First decomposition succeeded, jb split succesfully, ja did not split. Save subjet1=jba subjet2=jbb subjet3=ja"<<endl;
00186                 }
00187                 if (hardBreak1 && hardBreak2a &&  hardBreak2b) { 
00188                   hardA = jaa; 
00189                   hardB = jab; 
00190                   hardC = jba; 
00191                   hardD = jbb; 
00192                   if(verbose_)cout<<"First decomposition and both secondary decompositions succeeded. Save subjet1=jaa subjet2=jab subjet3=jba subjet4=jbb"<<endl;
00193                 }
00194                 
00195                 // check if we are left with >= 3 hard subjets
00196                 fastjet::PseudoJet subjet1 = blankJet;
00197                 fastjet::PseudoJet subjet2 = blankJet;
00198                 fastjet::PseudoJet subjet3 = blankJet;
00199                 fastjet::PseudoJet subjet4 = blankJet;
00200                 subjet1 = hardA; subjet2 = hardB; subjet3 = hardC; subjet4 = hardD;
00201                 
00202                 // record the hard subjets
00203                 vector<fastjet::PseudoJet> hardSubjets;
00204 
00205                 if ( verbose_ ) {
00206                   std::cout << "HardA : user_index = " << hardA.user_index() << ", (Pt,Y,Phi,M) = (" 
00207                             << hardA.pt() << ", " << hardA.rapidity() << ", " 
00208                             << hardA.phi() << ", " << hardA.m() << ")" << std::endl;
00209 
00210                   std::cout << "HardB : user_index = " << hardB.user_index() << ", (Pt,Y,Phi,M) = (" 
00211                             << hardB.pt() << ", " << hardB.rapidity() << ", " 
00212                             << hardB.phi() << ", " << hardB.m() << ")" << std::endl;
00213 
00214                   std::cout << "HardC : user_index = " << hardC.user_index() << ", (Pt,Y,Phi,M) = (" 
00215                             << hardC.pt() << ", " << hardC.rapidity() << ", " 
00216                             << hardC.phi() << ", " << hardC.m() << ")" << std::endl;
00217 
00218                   std::cout << "HardD : user_index = " << hardD.user_index() << ", (Pt,Y,Phi,M) = (" 
00219                             << hardD.pt() << ", " << hardD.rapidity() << ", " 
00220                             << hardD.phi() << ", " << hardD.m() << ")" << std::endl;
00221                 }
00222 
00223                 // Check to see if any subjects are counted amongst the "hard" subjets from previous
00224                 // lines. NOTE: In Fastjet 3.0, the default "user_index" changed from 0 to -1, so
00225                 // this can no longer be used as a designator for the veto of "blankJet" subjets,
00226                 // and now switch to pt > some small value. 
00227                 if ( subjet1.pt() > 0.0001 )
00228                         hardSubjets.push_back(subjet1);
00229                 if ( subjet2.pt() > 0.0001 )
00230                         hardSubjets.push_back(subjet2);
00231                 if ( subjet3.pt() > 0.0001 )
00232                         hardSubjets.push_back(subjet3);
00233                 if ( subjet4.pt() > 0.0001 )
00234                         hardSubjets.push_back(subjet4);
00235                 sort(hardSubjets.begin(), hardSubjets.end(), compEt );
00236 
00237                 // Use new fastjet functionality to create a Pseudojet from constituents
00238                 fastjet::PseudoJet candidate = join(hardSubjets);
00239                 // Reset the jet's 4-vector to the "ungroomed" value
00240                 candidate.reset_momentum( jetIt->px(), jetIt->py(), jetIt->pz(), jetIt->e() );
00241 
00242                 if ( verbose_ ) {
00243                   std::cout << "Final top-jet candidate: (Pt,Y,Phi,M) = (" 
00244                             << candidate.pt() << ", " << candidate.rapidity() << ", " 
00245                             << candidate.phi() << ", " << candidate.m() << ")" << std::endl;
00246                   std::vector<fastjet::PseudoJet> pieces = candidate.pieces();
00247                   std::cout << "Number of pieces = " << pieces.size() << std::endl;
00248                   for ( std::vector<fastjet::PseudoJet>::const_iterator ibegin = pieces.begin(), iend = pieces.end(), i = ibegin;
00249                         i != iend; ++i ) {
00250                     std::cout << "   Piece : " << i - ibegin << ", (Pt,Y,Phi,M) = (" 
00251                               << i->pt() << ", " << i->rapidity() << ", " 
00252                               << i->phi() << ", " << i->m() << ")" << std::endl;
00253                   }
00254                 }
00255                 // Add to the list
00256                 hardjetsOutput.push_back( candidate );          
00257         }
00258 }
00259 
00260 
00261 
00262 
00263 //-----------------------------------------------------------------------
00264 // determine whether two clusters (made of calorimeter towers) are living on "adjacent" cells.  if they are, then
00265 // we probably shouldn't consider them to be independent objects!
00266 //
00267 // From Sal: Ignoring genjet case
00268 //
00269 bool CATopJetAlgorithm::adjacentCells(const fastjet::PseudoJet & jet1, const fastjet::PseudoJet & jet2, 
00270                                                                           const vector<fastjet::PseudoJet> & cell_particles,
00271                                                                           const fastjet::ClusterSequence & theClusterSequence,
00272                                                                           double nCellMin ) const {
00273         
00274         
00275         double eta1 = jet1.rapidity();
00276         double phi1 = jet1.phi();
00277         double eta2 = jet2.rapidity();
00278         double phi2 = jet2.phi();
00279         
00280         double deta = abs(eta2 - eta1) / 0.087;
00281         double dphi = fabs( reco::deltaPhi(phi2,phi1) ) / 0.087;
00282         
00283         return ( ( deta + dphi ) <= nCellMin );
00284 }
00285 
00286 
00287 
00288 //-------------------------------------------------------------------------
00289 // attempt to decompose a jet into "hard" subjets, where hardness is set by ptHard
00290 //
00291 bool CATopJetAlgorithm::decomposeJet(const fastjet::PseudoJet & theJet, 
00292                                                                          const fastjet::ClusterSequence & theClusterSequence, 
00293                                                                          const vector<fastjet::PseudoJet> & cell_particles,
00294                                                                          double ptHard, double nCellMin, double deltarcut,
00295                                                                          fastjet::PseudoJet & ja, fastjet::PseudoJet & jb, 
00296                                                                          vector<fastjet::PseudoJet> & leftovers) const {
00297         
00298         bool goodBreak;
00299         fastjet::PseudoJet j = theJet;
00300         double InputObjectPt = j.perp();
00301         if ( verbose_ )cout<<"Input Object Pt = "<<InputObjectPt<<endl;
00302         if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
00303         leftovers.clear();
00304         if ( verbose_ )cout<<"start while loop"<<endl;
00305         
00306         while (1) {                                                      // watch out for infinite loop!
00307                 goodBreak = theClusterSequence.has_parents(j,ja,jb);
00308                 if (!goodBreak){
00309                         if ( verbose_ )cout<<"bad break. this is one cell. can't decluster anymore."<<endl;
00310                         break;         // this is one cell, can't decluster anymore
00311                 }
00312                 
00313                 if ( verbose_ )cout<<"good break. ja Pt = "<<ja.perp()<<" jb Pt = "<<jb.perp()<<endl;
00314                 
00316                 
00317                 // check if clusters are adjacent using a constant deltar adjacency.
00318                 double clusters_deltar=fabs(ja.eta()-jb.eta())+fabs(deltaPhi(ja.phi(),jb.phi()));
00319                 
00320                 if ( verbose_  && useAdjacency_ ==1)cout<<"clusters_deltar = "<<clusters_deltar<<endl;
00321                 if ( verbose_  && useAdjacency_ ==1)cout<<"deltar cut = "<<deltarcut<<endl;
00322                 
00323                 if ( useAdjacency_==1 && clusters_deltar < deltarcut){
00324                         if ( verbose_ )cout<<"clusters too close. consant adj. break."<<endl;
00325                         break;
00326                 } 
00327                 
00328                 // Check if clusters are adjacent using a DeltaR adjacency which is a function of pT.
00329                 double clusters_deltaR=deltaR( ja.eta(), ja.phi(), jb.eta(), jb.phi() );
00330                 
00331                 if ( verbose_  && useAdjacency_ ==2)cout<<"clusters_deltaR = "<<clusters_deltaR<<endl;
00332                 if ( verbose_  && useAdjacency_ ==2)cout<<"0.4-0.0004*InputObjectPt = "<<0.4-0.0004*InputObjectPt<<endl;
00333                 
00334                 if ( useAdjacency_==2 && clusters_deltaR < 0.4-0.0004*InputObjectPt)
00335                 {
00336                         if ( verbose_ )cout<<"clusters too close. modified adj. break."<<endl;
00337                         break;
00338                 } 
00339 
00340                 // Check if clusters are adjacent in the calorimeter. 
00341                 if ( useAdjacency_==3 &&  adjacentCells(ja,jb,cell_particles,theClusterSequence,nCellMin) ){                  
00342                         if ( verbose_ )cout<<"clusters too close in the calorimeter. calorimeter adj. break."<<endl;
00343                         break;         // the clusters are "adjacent" in the calorimeter => shouldn't have decomposed
00344                 }
00345                                 
00346                 if ( verbose_ )cout<<"clusters pass distance cut"<<endl;
00347                 
00349                 
00350                 if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
00351                 
00352                 if (ja.perp() < ptHard && jb.perp() < ptHard){
00353                         if ( verbose_ )cout<<"two soft clusters. dead end"<<endl;
00354                         break;         // broke into two soft clusters, dead end
00355                 }
00356                 
00357                 if (ja.perp() > ptHard && jb.perp() > ptHard){
00358                         if ( verbose_ )cout<<"two hard clusters. done"<<endl;
00359                         return true;   // broke into two hard clusters, we're done!
00360                 }
00361                 
00362                 else if (ja.perp() > jb.perp()) {                              // broke into one hard and one soft, ditch the soft one and try again
00363                         if ( verbose_ )cout<<"ja hard jb soft. try to split hard. j = ja"<<endl; 
00364                         j = ja;
00365                         vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(jb);
00366                         leftovers.insert(leftovers.end(),particles.begin(),particles.end());
00367                 }
00368                 else {
00369                         if ( verbose_ )cout<<"ja hard jb soft. try to split hard. j = jb"<<endl; 
00370                         j = jb;
00371                         vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(ja);
00372                         leftovers.insert(leftovers.end(),particles.begin(),particles.end());
00373                 }
00374         }
00375         
00376         if ( verbose_ )cout<<"did not decluster."<<endl;  // did not decluster into hard subjets
00377         
00378         ja.reset(0,0,0,0);
00379         jb.reset(0,0,0,0);
00380         leftovers.clear();
00381         return false;
00382 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoJets/JetAlgorithms/src/CATopJetAlgorithm.cc
