00001 #ifndef D0RunIIconeJets_ILCONEALGORITHM
00002 #define D0RunIIconeJets_ILCONEALGORITHM
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00044 #include <vector>
00045 #include <list>
00046 #include <utility>
00047 #include <map>
00048 #include <algorithm>
00049 #include <iostream>
00050
00051
00052
00053
00054 #include "ProtoJet.hpp"
00055
00056 #include "ConeSplitMerge.hpp"
00057
00058
00059 #include "inline_maths.h"
00060
00062 #include <fastjet/internal/base.hh>
00063
00064 FASTJET_BEGIN_NAMESPACE
00065
00066 namespace d0{
00067
00068 using namespace inline_maths;
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 #ifdef ILCA_USE_ORDERED_LIST
00087
00088 template <class Item>
00089 class rapidity_order
00090 {
00091 public:
00092 bool operator()(const Item* first, const Item* second)
00093 {
00094 return (first->y() < second->y());
00095 }
00096 bool operator()(float const & first, const Item* second)
00097 {
00098 return (first < second->y());
00099 }
00100 bool operator()(const Item* first, float const& second)
00101 {
00102 return (first->y() < second);
00103 }
00104 };
00105 #endif
00106
00107
00108
00109 template <class Item>
00110 class ILConeAlgorithm
00111 {
00112
00113 public:
00114
00115
00116
00117 ILConeAlgorithm():
00118 _CONE_RADIUS(0.),
00119 _MIN_JET_ET(99999.),
00120 _ET_MIN_RATIO(0.5),
00121 _FAR_RATIO(0.5),
00122 _SPLIT_RATIO(0.5),
00123 _DUPLICATE_DR(0.005),
00124 _DUPLICATE_DPT(0.01),
00125 _SEARCH_CONE(0.5),
00126 _PT_MIN_LEADING_PROTOJET(0),
00127 _PT_MIN_SECOND_PROTOJET(0),
00128 _MERGE_MAX(10000),
00129 _PT_MIN_noMERGE_MAX(0)
00130 {;}
00131
00132
00133 ILConeAlgorithm(float cone_radius, float min_jet_Et, float split_ratio,
00134 float far_ratio=0.5, float Et_min_ratio=0.5,
00135 bool kill_duplicate=true, float duplicate_dR=0.005,
00136 float duplicate_dPT=0.01, float search_factor=1.0,
00137 float pT_min_leading_protojet=0., float pT_min_second_protojet=0.,
00138 int merge_max=10000, float pT_min_nomerge=0.) :
00139
00140 _CONE_RADIUS(cone_radius),
00141
00142 _MIN_JET_ET(min_jet_Et),
00143
00144 _ET_MIN_RATIO(Et_min_ratio),
00145
00146 _FAR_RATIO(far_ratio),
00147
00148 _SPLIT_RATIO(split_ratio),
00149 _DUPLICATE_DR(duplicate_dR),
00150 _DUPLICATE_DPT(duplicate_dPT),
00151 _SEARCH_CONE(cone_radius/search_factor),
00152
00153
00154 _KILL_DUPLICATE(kill_duplicate),
00155 _PT_MIN_LEADING_PROTOJET(pT_min_leading_protojet),
00156 _PT_MIN_SECOND_PROTOJET(pT_min_second_protojet),
00157 _MERGE_MAX(merge_max),
00158 _PT_MIN_noMERGE_MAX(pT_min_nomerge)
00159 {;}
00160
00161
00162 ~ILConeAlgorithm() {;}
00163
00164
00165
00166 void makeClusters(
00167
00168 std::list<Item> &jets,
00169 std::list<const Item*>& itemlist,
00170
00172
00173
00174 const float Item_ET_Threshold);
00175
00176
00177 std::vector<ProtoJet<Item> > ilcv;
00178
00179 private:
00180
00181 float _CONE_RADIUS;
00182 float _MIN_JET_ET;
00183 float _ET_MIN_RATIO;
00184 float _FAR_RATIO;
00185 float _SPLIT_RATIO;
00186 float _DUPLICATE_DR;
00187 float _DUPLICATE_DPT;
00188 float _SEARCH_CONE;
00189
00190 bool _KILL_DUPLICATE;
00191
00192 float _PT_MIN_LEADING_PROTOJET;
00193 float _PT_MIN_SECOND_PROTOJET;
00194 int _MERGE_MAX;
00195 float _PT_MIN_noMERGE_MAX;
00196
00197
00198
00199
00200 class TemporaryJet : public ProtoJet<Item>
00201 {
00202
00203 public :
00204
00205 TemporaryJet(float seedET) : ProtoJet<Item>(seedET) {;}
00206
00207 TemporaryJet(float seedET,float y,float phi) :
00208 ProtoJet<Item>(seedET,y,phi) {;}
00209
00210 ~TemporaryJet() {;}
00211
00212 float dist(TemporaryJet& jet) const
00213 {
00214 return RDelta(this->_y,this->_phi,jet.y(),jet.phi());
00215 }
00216
00217 void midpoint(const TemporaryJet& jet,float & y, float & phi) const
00218 {
00219
00220
00221 float pTsum = this->_pT + jet.pT();
00222 y = (this->_y*this->_pT + jet.y()*jet.pT())/pTsum;
00223
00224 phi = (this->_phi*this->_pT + jet.phi()*jet.pT())/pTsum;
00225
00226
00227
00228
00229 if ( fabs(phi-this->_phi)>2.0 ) {
00230 phi = fmod( this->_phi+PI, TWOPI);
00231 if (phi < 0.0) phi += TWOPI;
00232 phi -= PI;
00233
00234 float temp=fmod( jet.phi()+PI, TWOPI);
00235 if (temp < 0.0) temp += TWOPI;
00236 temp -= PI;
00237
00238 phi = (phi*this->_pT + temp*jet.pT()) /pTsum;
00239 }
00240
00241 if ( phi < 0. ) phi += TWOPI;
00242 }
00243
00244
00246 #ifdef ILCA_USE_MMAP
00247 bool is_stable(const std::multimap<float,const Item*>& items,
00248 float radius, float min_ET, int max_iterations=50)
00249 #else
00250 bool is_stable(const std::list<const Item*>& itemlist, float radius,
00251 float min_ET, int max_iterations=50)
00252 #endif
00253
00254 {
00255 float radius2 = radius*radius;
00256 float Rcut= 1.E-06;
00257
00258
00259
00260 bool stable= true;
00261 int trial= 0;
00262 float Yst;
00263 float PHIst;
00264 do {
00265 trial++;
00266
00267 Yst = this->_y;
00268 PHIst= this->_phi;
00269
00270 this->erase();
00271
00272 this->setJet(Yst,PHIst,0.0);
00273
00274
00275 #ifdef ILCA_USE_ORDERED_LIST
00276 std::list<const Item*>::const_iterator lower =
00277 lower_bound(itemlist.begin(),itemlist.end(),Yst-radius,
00278 rapidity_order<Item>());
00279 std::list<const Item*>::const_iterator upper =
00280 upper_bound(itemlist.begin(),itemlist.end(),Yst+radius,
00281 rapidity_order<Item>());
00282 for(std::list<const Item*>::const_iterator tk = lower; tk != upper; ++tk) {
00283
00284 if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
00285 {
00286 addItem(*tk);
00287 }
00288 }
00289 #else
00290 #ifdef ILCA_USE_MMAP
00291
00292 for ( std::multimap<float,const Item*>::const_iterator
00293 tk = items.lower_bound(Yst-radius);
00294 tk != items.upper_bound(Yst+radius); ++tk )
00295 {
00296
00297 if(RD2(((*tk).second)->y(),((*tk).second)->phi(),Yst,PHIst) <= radius2)
00298 {
00299 addItem((*tk).second);
00300 }
00301 }
00302
00303 #else
00304
00305
00306 for(typename std::list<const Item*>::const_iterator tk = itemlist.begin(); tk != itemlist.end(); ++tk)
00307 {
00308
00309 if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
00310 {
00311
00312 addItem(*tk);
00313 }
00314 }
00315 #endif
00316 #endif
00317
00318
00319 this->updateJet();
00320
00321
00322 if(this->_pT < min_ET )
00323 {
00324 stable= false;
00325 break;
00326 }
00327
00328 } while(RD2(this->_y,this->_phi,Yst,PHIst) >= Rcut && trial <= max_iterations);
00329
00330 return stable;
00331 }
00332
00333 private :
00334
00335 };
00336 };
00338
00339
00340 template <class Item>
00341 void ILConeAlgorithm <Item>::
00342
00343 makeClusters(
00344 std::list<Item> &jets,
00345 std::list<const Item*>& ilist,
00346
00348
00349
00350 const float Item_ET_Threshold)
00351 {
00352
00353 for ( typename std::list<const Item*>::iterator it = ilist.begin();
00354
00355 it != ilist.end(); )
00356 {
00357 if ( (*it)->pT() < Item_ET_Threshold )
00358 {
00359 it = ilist.erase(it);
00360 }
00361 else ++it;
00362 }
00363
00364
00365
00366
00367
00369 std::vector<const Item*> ecv;
00370 for ( typename std::list<const Item*>::iterator it = ilist.begin();
00371 it != ilist.end(); it++) {
00372 ecv.push_back(*it);
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382 float far_def = _FAR_RATIO*_CONE_RADIUS * _FAR_RATIO*_CONE_RADIUS;
00383
00384
00385 float ratio= _MIN_JET_ET*_ET_MIN_RATIO;
00386
00387
00388 #ifdef ILCA_USE_ORDERED_LIST
00389
00390 ilist.sort(rapidity_order<Item>());
00391 #else
00392 #ifdef ILCA_USE_MMAP
00393
00394 std::multimap<float,const Item*> items;
00395 std::list<const Item*>::const_iterator it;
00396 for(it = ilist.begin(); it != ilist.end(); ++it)
00397 {
00398 pair<float,const Item*> p((*it)->y(),*it);
00399 items.insert(p);
00400 }
00401 #endif
00402 #endif
00403
00404 std::vector<ProtoJet<Item> > mcoll;
00405 std::vector<TemporaryJet> scoll;
00406
00407
00408
00409
00410 typename std::vector<const Item*>::iterator jclu;
00411 for(jclu = ecv.begin(); jclu != ecv.end(); ++jclu)
00412 {
00413
00414 const Item* ptr= *jclu;
00415 float p[4];
00416 ptr->p4vec(p);
00417 float Yst = P2y(p);
00418 float PHIst= P2phi(p);
00419
00420
00421 bool is_far= true;
00422
00423 for(unsigned int i = 0; i < scoll.size(); ++i)
00424 {
00425 float y = scoll[i].y();
00426 float phi= scoll[i].phi();
00427 if(RD2(Yst,PHIst,y,phi) < far_def)
00428 {
00429 is_far= false;
00430 break;
00431 }
00432 }
00433
00434
00435 if(is_far)
00436 {
00437 TemporaryJet jet(ptr->pT(),Yst,PHIst);
00438
00439
00440
00441
00442
00443 #ifdef ILCA_USE_MMAP
00444 if(jet.is_stable(items,_CONE_RADIUS,ratio,0) && jet.is_stable(items,_SEARCH_CONE,3.0))
00445 #else
00446 if(jet.is_stable(ilist,_CONE_RADIUS,ratio,0) && jet.is_stable(ilist,_SEARCH_CONE,3.0))
00447 #endif
00448 {
00449
00450
00451
00452
00453 #ifdef ILCA_USE_MMAP
00454 jet.is_stable(items,_CONE_RADIUS,ratio,0) ;
00455 #else
00456 jet.is_stable(ilist,_CONE_RADIUS,ratio,0) ;
00457 #endif
00458
00459
00460 if ( _KILL_DUPLICATE )
00461 {
00462
00463 float distmax = 999.; int imax = -1;
00464 for(unsigned int i = 0; i < scoll.size(); ++i)
00465 {
00466 float dist = jet.dist(scoll[i]);
00467 if ( dist < distmax )
00468 {
00469 distmax = dist;
00470 imax = i;
00471 }
00472 }
00473 if ( distmax > _DUPLICATE_DR ||
00474 fabs((jet.pT()-scoll[imax].pT())/scoll[imax].pT())>_DUPLICATE_DPT )
00475 {
00476 scoll.push_back(jet);
00477 mcoll.push_back(jet);
00478
00479 }
00480
00481
00482
00483
00484
00485
00486
00487 }
00488 else
00489 {
00490 scoll.push_back(jet);
00491 mcoll.push_back(jet);
00492
00493 }
00494
00495 }
00496 }
00497 }
00498
00499
00500 for(unsigned int i = 0; i < scoll.size(); ++i)
00501 {
00502 for(unsigned int k = i+1; k < scoll.size(); ++k)
00503 {
00504 float djet= scoll[i].dist(scoll[k]);
00505 if(djet > _CONE_RADIUS && djet < 2.*_CONE_RADIUS)
00506 {
00507 float y_mid,phi_mid;
00508 scoll[i].midpoint(scoll[k],y_mid,phi_mid);
00509 TemporaryJet jet(-999999.,y_mid,phi_mid);
00510
00511
00512
00513 #ifdef ILCA_USE_MMAP
00514 if(jet.is_stable(items,_CONE_RADIUS,ratio))
00515 #else
00516 if(jet.is_stable(ilist,_CONE_RADIUS,ratio))
00517 #endif
00518 {
00519 mcoll.push_back(jet);
00520
00521 }
00522 }
00523 }
00524 }
00525
00526
00527
00528 ConeSplitMerge<Item> pjets(mcoll);
00529 ilcv.clear();
00530 pjets.split_merge(ilcv,_SPLIT_RATIO, _PT_MIN_LEADING_PROTOJET, _PT_MIN_SECOND_PROTOJET, _MERGE_MAX, _PT_MIN_noMERGE_MAX);
00531
00532
00533 for(unsigned int i = 0; i < ilcv.size(); ++i)
00534 {
00535 if ( ilcv[i].pT() > _MIN_JET_ET )
00536 {
00537
00538 Item ptrclu;
00539
00540
00541 std::list<const Item*> tlist=ilcv[i].LItems();
00542 typename std::list<const Item*>::iterator tk;
00543 for(tk = tlist.begin(); tk != tlist.end(); ++tk)
00544 {
00545
00546 float pk[4];
00547 (*tk)->p4vec(pk);
00548
00549
00550
00551
00552 ptrclu.Add(**tk);
00553 }
00554
00555
00556
00557
00558 jets.push_back(ptrclu);
00559 }
00560 }
00561 }
00562
00563 }
00564
00565 FASTJET_END_NAMESPACE
00566
00567 #endif
00568
00569