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