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