FastJet 3.0.0
|
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