FastJet 3.0.0
JetConeFinderTool.cc
00001 //*******************************************************************************
00002 // Filename : JetConeFinderTool.cc 
00003 // Author   : Ambreesh Gupta
00004 // Created  : Nov, 2000
00005 //
00006 // Jan 2004: Use CLHEP units. Use phi = (-pi,pi].
00007 //
00008 // File taken from SpartyJet v2.20.0
00009 // Modifications:
00010 //   removed the string name in the ctor
00011 //   removed the Message m_log
00012 //   replaced StatusCode by an int (as defined previously)
00013 //   cleaned the comments
00014 //*******************************************************************************
00015 
00016 #include <iostream>
00017 
00018 #include "JetConeFinderTool.hh"
00019 #include "Jet.hh"
00020 #include "JetDistances.hh"
00021 #include "CommonUtils.hh"
00022 
00023 #include <vector>
00024 #include <math.h>
00025 
00026 #include <fastjet/internal/base.hh>
00027 
00028 FASTJET_BEGIN_NAMESPACE
00029 
00030 namespace atlas { 
00031 
00032 // set the default energy scale
00033   double GeV = 1.0; //1000;
00034 
00035 JetConeFinderTool::JetConeFinderTool() :m_coneR(0.7)
00036   , m_ptcut(0.0*GeV)
00037   , m_eps(0.05)
00038   , m_seedPt(2.0*GeV)
00039   , m_etaMax(5.0)
00040 {}
00041 
00042 JetConeFinderTool::~JetConeFinderTool()
00043 {}
00044 
00045 /////////////////////////////////////////////////////////////////////////////////
00046 //Execution                                                                     /
00047 /////////////////////////////////////////////////////////////////////////////////
00048 int JetConeFinderTool::execute(jetcollection_t & theJets)
00049 {
00050   sort_jet_list<JetSorter_Et>(theJets);
00051 
00052   m_pjetV = &theJets;
00053 
00054   if(theJets.size()==0) return 0;
00055 
00056   // Initiale ctr/dctr counter for object counting.
00057   m_ctr = 0;
00058   m_dctr = 0;
00059 
00060   //////////////////////
00061   // Reconstruct Jets // 
00062   //////////////////////
00063   this->reconstruct();
00064 
00065   //////////////////////////
00066   // ReFill JetCollection //
00067   //////////////////////////
00068   clear_list(theJets);
00069   jetcollection_t::iterator it = m_jetOV->begin();
00070   jetcollection_t::iterator itE = m_jetOV->end();
00071   for(; it!=itE; ++it){
00072     theJets.push_back(*it);
00073   }
00074 
00075 
00076   delete m_jetOV;
00077   return 1;
00078 }
00079 
00080 ///////////////////////////////////////////////////////////////////////////////
00081 // Reconstruction algorithm specific methods                                  /
00082 ///////////////////////////////////////////////////////////////////////////////
00083 
00084 void
00085 JetConeFinderTool::reconstruct()
00086 {
00087   m_jetOV = new jetcollection_t();
00088 
00089   jetcollection_t::iterator tItr;
00090   jetcollection_t::iterator tItr_begin = m_pjetV->begin();
00091   jetcollection_t::iterator tItr_end   = m_pjetV->end();
00092 
00093   // order towers in pt
00094 
00095   for ( tItr=tItr_begin; tItr!=tItr_end; ++tItr ) {    
00096 
00097     // Seed Cut
00098     double tEt = (*tItr)->et();
00099     if ( tEt < m_seedPt ) break;
00100 
00101     // Tower eta, phi
00102     double etaT = (*tItr)->eta();
00103     double phiT = (*tItr)->phi();    
00104     
00105     // Iteration logic
00106     bool stable = false;
00107     bool inGeom = true;
00108     
00109     Jet* preJet;
00110     
00111     int count = 1;
00112     do { // Iteration Loop
00113 
00114       // Make cone  
00115       preJet = calc_cone(etaT,phiT);
00116       double etaC = preJet->eta();
00117       double phiC = preJet->phi();
00118       
00119       double deta = fabs(etaT - etaC);
00120       double dphi = fabs(JetDistances::deltaPhi(phiT,phiC));
00121 
00122       // Is Stable ?
00123       if ( deta < m_eps && dphi < m_eps ) 
00124         stable = true;
00125       
00126       // In Geometry ?
00127       if ( fabs(etaC) > m_etaMax ) 
00128         inGeom = false; 
00129 
00130       etaT = etaC;
00131       phiT = phiC;
00132 
00133       if ( !stable && inGeom ) {
00134         delete preJet;
00135         m_dctr +=1;
00136       }
00137       ++count;
00138 
00139     }while ( !stable && inGeom && count < 10  );      
00140   
00141     if ( count > 9 && (!stable && inGeom) ) continue;  // FIXME 9 ?
00142 
00143     // If iteration was succesfull -- check if this is a new jet and 
00144     // add it to OV.
00145 
00146     if ( stable && inGeom ) {
00147       jetcollection_t::iterator pItr   = m_jetOV->begin();
00148       jetcollection_t::iterator pItrE  = m_jetOV->end();
00149     
00150       bool newJet = true;
00151       double etaT = preJet->eta();
00152       double phiT = preJet->phi();
00153 
00154       for ( ; pItr != pItrE ; ++pItr ) {
00155         double etaC = (*pItr)->eta();
00156         double phiC = (*pItr)->phi();
00157 
00158         double deta = fabs(etaT - etaC);
00159         double dphi = fabs(JetDistances::deltaPhi(phiT,phiC));
00160 
00161         if ( deta < 0.05  && dphi < 0.05 ) { 
00162           // Debugging done by Gregory Soyez:
00163           //
00164           // Becase of the cut on the Et difference imposed on the
00165           // ordering (line 80 of Jet.hh), the ordering of the input
00166           // particles in Et is not robust agains different
00167           // implementations of sort (which has an undefined behaviour
00168           // if 2 particles are equal). A consequence of this is that
00169           // stable cone search will consider these 2 seeds in an
00170           // undefined order. If the 2 resulting stable cones are too
00171           // close (deta<0.05, dphi<0.05) one will be accepted and the
00172           // other rejected. Which one depends on the ordering and is
00173           // thus undefined. If the 2 stable cones do not have the
00174           // same number of constituents this could affect the result
00175           // of the clustering.
00176           //
00177           // The line below helps debugging these cases by printing
00178           // the rejected stable cones
00179           //std::cout << "rejecting " << etaT << " " << phiT << " " << preJet->et() << (*tItr)->eta() << " " << (*tItr)->phi() << " " << (*tItr)->et() << std::endl;
00180           newJet = false;
00181           break;
00182         }
00183       }
00184       if ( newJet ) {
00185         m_jetOV->push_back( preJet );
00186         // Debugging done by Gregory Soyez:
00187         //
00188         // Becase of the cut on the Et difference imposed on the
00189         // ordering (line 80 of Jet.hh), the ordering of the input
00190         // particles in Et is not robust agains different
00191         // implementations of sort (which has an undefined behaviour
00192         // if 2 particles are equal). A consequence of this is that
00193         // stable cone search will consider these 2 seeds in an
00194         // undefined order. If the 2 resulting stable cones are too
00195         // close (deta<0.05, dphi<0.05) one will be accepted and the
00196         // other rejected. Which one depends on the ordering and is
00197         // thus undefined. If the 2 stable cones do not have the
00198         // same number of constituents this could affect the result
00199         // of the clustering.
00200         //
00201         // The line below helps debugging these cases by printing
00202         // the accepted stable cones
00203         //std::cout << "accepting " << etaT << " " << phiT << " " << preJet->et() << (*tItr)->eta() << " " << (*tItr)->phi() << " " << (*tItr)->et() << std::endl;
00204       }
00205       else {
00206         delete preJet;
00207         m_dctr +=1;
00208       }
00209     }
00210     else {
00211       delete preJet;
00212       m_dctr +=1;
00213     }
00214   }    
00215 }
00216 
00217 Jet* JetConeFinderTool::calc_cone(double eta, double phi)
00218 {
00219   // Create a new Jet   
00220   Jet* j = new Jet();
00221   m_ctr +=1; 
00222 
00223   // Add all ProtoJet within m_coneR to this Jet  
00224   jetcollection_t::iterator itr  = m_pjetV->begin();
00225   jetcollection_t::iterator itrE = m_pjetV->end();
00226 
00227   for ( ; itr!=itrE; ++itr ) {
00228     double dR = JetDistances::deltaR(eta,phi,(*itr)->eta(),(*itr)->phi());
00229     if ( dR < m_coneR ) {
00230       j->addJet( (*itr) );
00231     }
00232   }    
00233 
00234   return j;
00235 }
00236 
00237 
00238 
00239 }  // namespace atlas
00240 
00241 FASTJET_END_NAMESPACE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends