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