FastJet 3.0.2
|
00001 #ifndef D0RunIIconeJets_CONESPLITMERGE 00002 #define D0RunIIconeJets_CONESPLITMERGE 00003 // --------------------------------------------------------------------------- 00004 // ConeSplitMerge.hpp 00005 // 00006 // Created: 28-JUL-2000 Francois Touze 00007 // 00008 // Purpose: Implements the pT ordered jet split/merge algorithm for the 00009 // Improved Legacy Cone Algorithm split/merge algo. 00010 // 00011 // Modified: 00012 // 31-JUL-2000 Laurent Duflot 00013 // + introduce support for additional informations (ConeJetInfo) 00014 // 1-AUG-2000 Laurent Duflot 00015 // + jet merge criteria was wrong, now calculate shared_ET and compare to 00016 // neighbour jet pT. 00017 // + split was incomplete: shared items not really removed from jets. 00018 // 4-Aug-2000 Laurent Duflot 00019 // + use item methods y() and phi() rather than p4vec() and then P2y and 00020 // P2phi 00021 // 7-Aug-2000 Laurent Duflot 00022 // + force the list to be organized by decreasing ET and, for identical ET, 00023 // by decreasing seedET. Identical protojets can be found by eg nearby 00024 // seeds. The seedET ordering is such that for identical jets, the one 00025 // with the highest seedET is kept, which is what we want for efficiency 00026 // calculations. 00027 // + avoid unecessary copies of lists by using reference when possible 00028 // 9-Aug-2000 Laurent Duflot 00029 // + save initial jet ET before split/merge 00030 // 1-May-2007 Lars Sonnenschein 00031 // extracted from D0 software framework and modified to remove subsequent dependencies 00032 // 00033 // 00034 // This file is distributed with FastJet under the terms of the GNU 00035 // General Public License (v2). Permission to do so has been granted 00036 // by Lars Sonnenschein and the D0 collaboration (see COPYING for 00037 // details) 00038 // 00039 // History of changes in FastJet compared tothe original version of 00040 // ConeSplitMerge.hpp 00041 // 00042 // 2011-12-13 Gregory Soyez <soyez@fastjet.fr> 00043 // 00044 // * added license information 00045 // 00046 // 2009-01-17 Gregory Soyez <soyez@fastjet.fr> 00047 // 00048 // * put the code in the fastjet::d0 namespace 00049 // 00050 // 2007-12-14 Gavin Salam <salam@lpthe.jussieu.fr> 00051 // 00052 // * replaced make_pair by std::make_pair 00053 // 00054 // --------------------------------------------------------------------------- 00055 00056 00057 #include <iostream> 00058 #include <map> 00059 #include <utility> 00060 #include <vector> 00061 #include "ProtoJet.hpp" 00062 00063 //using namespace D0RunIIconeJets_CONEJETINFO; 00064 00065 #include <fastjet/internal/base.hh> 00066 00067 FASTJET_BEGIN_NAMESPACE 00068 00069 namespace d0{ 00070 00071 // 00072 // this class is used to order ProtoJets by decreasing ET and seed ET 00073 template <class Item> 00074 class ProtoJet_ET_seedET_order 00075 { 00076 public: 00077 bool operator()(const ProtoJet<Item> & first, const ProtoJet<Item> & second) 00078 { 00079 if ( first.pT() > second.pT() ) return true; 00080 else 00081 if ( first.pT() < second.pT() ) return false; 00082 else return ( first.info().seedET() > second.info().seedET() ); 00083 } 00084 }; 00085 00086 00087 template <class Item> 00088 class ConeSplitMerge { 00089 00090 public : 00091 00092 typedef std::multimap<ProtoJet<Item>,float,ProtoJet_ET_seedET_order<Item> > PJMMAP; 00093 00094 ConeSplitMerge(); 00095 ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector); 00096 ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist); 00097 ~ConeSplitMerge() {;} 00098 void split_merge(std::vector<ProtoJet<Item> >& ecv,float s, float pT_min_leading_protojet, float pT_min_second_protojet, int MergeMax, float pT_min_noMergeMax); 00099 00100 private : 00101 PJMMAP _members; 00102 }; 00103 /////////////////////////////////////////////////////////////////////////////// 00104 template<class Item> 00105 ConeSplitMerge<Item>::ConeSplitMerge():_members() {;} 00106 00107 template<class Item> 00108 ConeSplitMerge<Item>::ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector) 00109 { 00110 // sort proto_jets in Et descending order 00111 typename std::vector<ProtoJet<Item> >::const_iterator jt; 00112 for(jt = jvector.begin(); jt != jvector.end(); ++jt) 00113 { 00114 // this is supposed to be a stable cone, declare so 00115 ProtoJet<Item> jet(*jt); 00116 jet.NowStable(); 00117 _members.insert(std::make_pair(jet,jet.pT())); 00118 } 00119 } 00120 00121 template<class Item> 00122 ConeSplitMerge<Item>::ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist) 00123 { 00124 //_max_nb_items =-1; 00125 // sort proto_jets in Et descending order 00126 typename std::list<ProtoJet<Item> >::const_iterator jt; 00127 for(jt = jlist.begin(); jt != jlist.end(); ++jt) 00128 { 00129 // this is supposed to be a stable cone, declare so 00130 ProtoJet<Item> jet(*jt); 00131 jet.NowStable(); 00132 _members.insert(std::make_pair(jet,jet.pT())); 00133 } 00134 } 00135 00136 template<class Item> 00137 void ConeSplitMerge<Item>::split_merge(std::vector<ProtoJet<Item> >& jcv, 00138 float shared_ET_fraction, 00139 float pT_min_leading_protojet, float pT_min_second_protojet, 00140 int MergeMax, float pT_min_noMergeMax) 00141 { 00142 while(!_members.empty()) 00143 { 00144 /* 00145 { 00146 std::cout << std::endl; 00147 std::cout << " --------------- list of protojets ------------------ " <<std::endl; 00148 for ( PJMMAP::iterator it = _members.begin(); 00149 it != _members.end(); ++it) 00150 { 00151 std::cout << " pT y phi " << (*it).pT() << " " << (*it).y() << " " << (*it).phi() << " " << (*it).info().seedET() << " " << (*it).info().nbMerge() << " " << (*it).info().nbSplit() << std::endl; 00152 } 00153 std::cout << " ----------------------------------------------------- " <<std::endl; 00154 } 00155 */ 00156 00157 00158 // select highest Et jet 00159 typename PJMMAP::iterator itmax= _members.begin(); 00160 ProtoJet<Item> imax((*itmax).first); 00161 const std::list<const Item*>& ilist(imax.LItems()); 00162 00163 _members.erase(itmax); 00164 00165 // does jet share items? 00166 bool share= false; 00167 float shared_ET = 0.; 00168 typename PJMMAP::iterator jtmax; 00169 typename PJMMAP::iterator jt; 00170 for(jt = _members.begin(); jt != _members.end(); ++jt) 00171 { 00172 const std::list<const Item*>& jlist((*jt).first.LItems()); 00173 typename std::list<const Item*>::const_iterator tk; 00174 int count; 00175 for(tk = ilist.begin(), count = 0; tk != ilist.end(); 00176 ++tk, ++count) 00177 { 00178 typename std::list<const Item*>::const_iterator where= 00179 find(jlist.begin(),jlist.end(),*tk); 00180 if(where != jlist.end()) 00181 { 00182 share= true; 00183 shared_ET += (*tk)->pT(); 00184 } 00185 } 00186 if(share) 00187 { 00188 jtmax = jt; 00189 break; 00190 } 00191 } 00192 00193 if(!share) { 00194 // add jet to the final jet list 00195 jcv.push_back(imax); 00196 //std::cout << " final jet " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 00197 } 00198 else 00199 { 00200 // find highest Et neighbor 00201 ProtoJet<Item> jmax((*jtmax).first); 00202 00203 // drop neighbor 00204 _members.erase(jtmax); 00205 00206 00207 //std::cout << " split or merge ? " << imax.pT() << " " << shared_ET << " " << jmax.pT() << " " << s << " " << (jmax.pT())*s << std::endl; 00208 00209 // merge 00210 if( shared_ET > (jmax.pT())*shared_ET_fraction 00211 && (imax.pT() > pT_min_leading_protojet || jmax.pT() > pT_min_second_protojet) 00212 && (imax.info().nbMerge() < MergeMax || imax.pT() > pT_min_noMergeMax)) 00213 { 00214 // add neighbor's items to imax 00215 const std::list<const Item*>& jlist(jmax.LItems()); 00216 typename std::list<const Item*>::const_iterator tk; 00217 typename std::list<const Item*>::const_iterator iend= ilist.end(); 00218 bool same = true; // is jmax just the same as imax ? 00219 for(tk = jlist.begin(); tk != jlist.end(); ++tk) 00220 { 00221 typename std::list<const Item*>::const_iterator where= 00222 find(ilist.begin(),iend,*tk); 00223 if(where == iend) 00224 { 00225 imax.addItem(*tk); 00226 same = false; 00227 } 00228 } 00229 if ( ! same ) 00230 { 00231 // recalculate 00232 //float old_pT = imax.pT(); 00233 00234 imax.updateJet(); 00235 imax.merged(); 00236 //std::cout << " jet merge :: " << old_pT << " " << jmax.pT() << " " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 00237 } 00238 } 00239 00240 //split and assign removed shared cells from lowest pT protojet 00241 else if(shared_ET > (jmax.pT())*shared_ET_fraction) 00242 { 00243 // here we need to pull the lists, because there are items to remove 00244 std::list<const Item*> ilist(imax.LItems()); 00245 std::list<const Item*> jlist(jmax.LItems()); 00246 00247 typename std::list<const Item*>::iterator tk; 00248 for(tk = jlist.begin(); tk != jlist.end(); ) 00249 { 00250 typename std::list<const Item*>::iterator where= 00251 find(ilist.begin(),ilist.end(),*tk); 00252 if(where != ilist.end()) { 00253 tk = jlist.erase(tk); 00254 } 00255 else ++tk; 00256 } 00257 00258 jmax.erase(); 00259 00260 for ( typename std::list<const Item*>::const_iterator it = jlist.begin(); 00261 it != jlist.end(); ++it) jmax.addItem(*it); 00262 00263 // recalculated jet quantities 00264 jmax.updateJet(); 00265 jmax.splitted(); 00266 //std::cout << " jet split 1 :: " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl; 00267 _members.insert(std::make_pair(jmax,jmax.pT())); 00268 } 00269 00270 // split and assign shared cells to nearest protojet 00271 else 00272 { 00273 // here we need to pull the lists, because there are items to remove 00274 std::list<const Item*> ilist(imax.LItems()); 00275 std::list<const Item*> jlist(jmax.LItems()); 00276 00277 typename std::list<const Item*>::iterator tk; 00278 for(tk = jlist.begin(); tk != jlist.end(); ) 00279 { 00280 typename std::list<const Item*>::iterator where= 00281 find(ilist.begin(),ilist.end(),*tk); 00282 if(where != ilist.end()) { 00283 float yk = (*tk)->y(); 00284 float phik = (*tk)->phi(); 00285 float di= RD2(imax.y(),imax.phi(),yk,phik); 00286 float dj= RD2(jmax.y(),jmax.phi(),yk,phik); 00287 if(dj > di) 00288 { 00289 tk = jlist.erase(tk); 00290 //std::cout << " shared item " << (*tk)->pT() << " removed from neighbour jet " << std::endl; 00291 } 00292 else 00293 { 00294 ilist.erase(where); 00295 ++tk; 00296 //std::cout << " shared item " << (*tk)->pT() << " removed from leading jet " << std::endl; 00297 } 00298 } 00299 else ++tk; 00300 } 00301 // recalculate jets imax and jmax 00302 00303 // first erase all items 00304 imax.erase(); 00305 // put items that are left 00306 for ( typename std::list<const Item*>::const_iterator it = ilist.begin(); 00307 it != ilist.end(); ++it) imax.addItem(*it); 00308 // recalculated jet quantities 00309 imax.updateJet(); 00310 imax.splitted(); 00311 //std::cout << " jet split 2 :: " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 00312 00313 00314 // first erase all items 00315 jmax.erase(); 00316 // put items that are left 00317 for ( typename std::list<const Item*>::const_iterator it = jlist.begin(); 00318 it != jlist.end(); ++it) jmax.addItem(*it); 00319 // recalculated jet quantities 00320 jmax.updateJet(); 00321 jmax.splitted(); 00322 //std::cout << " jet split " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl; 00323 00324 _members.insert(std::make_pair(jmax,jmax.pT())); 00325 } 00326 _members.insert(std::make_pair(imax,imax.pT())); 00327 } 00328 } // while 00329 } 00330 /////////////////////////////////////////////////////////////////////////////// 00331 00332 } // namespace d0 00333 00334 FASTJET_END_NAMESPACE 00335 00336 #endif