1 #ifndef D0RunIIconeJets_ILCONEALGORITHM
2 #define D0RunIIconeJets_ILCONEALGORITHM
82 #include "ProtoJet.hpp"
84 #include "ConeSplitMerge.hpp"
87 #include "inline_maths.h"
90 #include <fastjet/internal/base.hh>
92 FASTJET_BEGIN_NAMESPACE
96 using namespace inline_maths;
114 #ifdef ILCA_USE_ORDERED_LIST
116 template <
class Item>
120 bool operator()(
const Item* first,
const Item* second)
122 return (first->y() < second->y());
124 bool operator()(
float const & first,
const Item* second)
126 return (first < second->y());
128 bool operator()(
const Item* first,
float const& second)
130 return (first->y() < second);
137 template <
class Item>
138 class ILConeAlgorithm
151 _DUPLICATE_DR(0.005),
152 _DUPLICATE_DPT(0.01),
154 _PT_MIN_LEADING_PROTOJET(0),
155 _PT_MIN_SECOND_PROTOJET(0),
157 _PT_MIN_noMERGE_MAX(0)
161 ILConeAlgorithm(
float cone_radius,
float min_jet_Et,
float split_ratio,
162 float far_ratio=0.5,
float Et_min_ratio=0.5,
163 bool kill_duplicate=
true,
float duplicate_dR=0.005,
164 float duplicate_dPT=0.01,
float search_factor=1.0,
165 float pT_min_leading_protojet=0.,
float pT_min_second_protojet=0.,
166 int merge_max=10000,
float pT_min_nomerge=0.) :
168 _CONE_RADIUS(cone_radius),
170 _MIN_JET_ET(min_jet_Et),
172 _ET_MIN_RATIO(Et_min_ratio),
174 _FAR_RATIO(far_ratio),
176 _SPLIT_RATIO(split_ratio),
177 _DUPLICATE_DR(duplicate_dR),
178 _DUPLICATE_DPT(duplicate_dPT),
179 _SEARCH_CONE(cone_radius/search_factor),
182 _KILL_DUPLICATE(kill_duplicate),
183 _PT_MIN_LEADING_PROTOJET(pT_min_leading_protojet),
184 _PT_MIN_SECOND_PROTOJET(pT_min_second_protojet),
185 _MERGE_MAX(merge_max),
186 _PT_MIN_noMERGE_MAX(pT_min_nomerge)
190 ~ILConeAlgorithm() {;}
196 std::list<Item> &jets,
197 std::list<const Item*>& itemlist,
202 const float Item_ET_Threshold);
205 std::vector<ProtoJet<Item> > ilcv;
215 float _DUPLICATE_DPT;
218 bool _KILL_DUPLICATE;
220 float _PT_MIN_LEADING_PROTOJET;
221 float _PT_MIN_SECOND_PROTOJET;
223 float _PT_MIN_noMERGE_MAX;
228 class TemporaryJet :
public ProtoJet<Item>
233 TemporaryJet(
float seedET) : ProtoJet<Item>(seedET) {;}
235 TemporaryJet(
float seedET,
float y_in,
float phi_in) :
236 ProtoJet<Item>(seedET,y_in,phi_in) {;}
240 float dist(TemporaryJet& jet)
const
242 return RDelta(this->_y,this->_phi,jet.y(),jet.phi());
245 void midpoint(
const TemporaryJet& jet,
float & y_out,
float & phi_out)
const
249 float pTsum = this->_pT + jet.pT();
250 y_out = (this->_y*this->_pT + jet.y()*jet.pT())/pTsum;
252 phi_out = (this->_phi*this->_pT + jet.phi()*jet.pT())/pTsum;
257 if ( fabs(phi_out-this->_phi)>2.0 ) {
258 phi_out = fmod( this->_phi+PI, TWOPI);
259 if (phi_out < 0.0) phi_out += TWOPI;
262 float temp=fmod( jet.phi()+PI, TWOPI);
263 if (temp < 0.0) temp += TWOPI;
266 phi_out = (phi_out*this->_pT + temp*jet.pT()) /pTsum;
269 if ( phi_out < 0. ) phi_out += TWOPI;
275 bool is_stable(
const std::multimap<float,const Item*>& items,
276 float radius,
float min_ET,
int max_iterations=50)
278 bool is_stable(
const std::list<const Item*>& itemlist,
float radius,
279 float min_ET,
int max_iterations=50)
283 float radius2 = radius*radius;
300 this->setJet(Yst,PHIst,0.0);
303 #ifdef ILCA_USE_ORDERED_LIST
304 std::list<const Item*>::const_iterator lower =
305 lower_bound(itemlist.begin(),itemlist.end(),Yst-radius,
306 rapidity_order<Item>());
307 std::list<const Item*>::const_iterator upper =
308 upper_bound(itemlist.begin(),itemlist.end(),Yst+radius,
309 rapidity_order<Item>());
310 for(std::list<const Item*>::const_iterator tk = lower; tk != upper; ++tk) {
312 if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
320 for ( std::multimap<float,const Item*>::const_iterator
321 tk = items.lower_bound(Yst-radius);
322 tk != items.upper_bound(Yst+radius); ++tk )
325 if(RD2(((*tk).second)->y(),((*tk).second)->phi(),Yst,PHIst) <= radius2)
327 addItem((*tk).second);
334 for(
typename std::list<const Item*>::const_iterator tk = itemlist.begin(); tk != itemlist.end(); ++tk)
337 if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
350 if(this->_pT < min_ET )
356 }
while(RD2(this->_y,this->_phi,Yst,PHIst) >= Rcut && trial <= max_iterations);
368 template <
class Item>
369 void ILConeAlgorithm <Item>::
372 std::list<Item> &jets,
373 std::list<const Item*>& ilist,
378 const float Item_ET_Threshold)
381 for (
typename std::list<const Item*>::iterator it = ilist.begin();
385 if ( (*it)->pT() < Item_ET_Threshold )
387 it = ilist.erase(it);
397 std::vector<const Item*> ecv;
398 for (
typename std::list<const Item*>::iterator it = ilist.begin();
399 it != ilist.end(); it++) {
410 float far_def = _FAR_RATIO*_CONE_RADIUS * _FAR_RATIO*_CONE_RADIUS;
413 float ratio= _MIN_JET_ET*_ET_MIN_RATIO;
416 #ifdef ILCA_USE_ORDERED_LIST
418 ilist.sort(rapidity_order<Item>());
422 std::multimap<float,const Item*> items;
423 std::list<const Item*>::const_iterator it;
424 for(it = ilist.begin(); it != ilist.end(); ++it)
426 pair<float,const Item*> p((*it)->y(),*it);
432 std::vector<ProtoJet<Item> > mcoll;
433 std::vector<TemporaryJet> scoll;
438 typename std::vector<const Item*>::iterator jclu;
439 for(jclu = ecv.begin(); jclu != ecv.end(); ++jclu)
442 const Item* ptr= *jclu;
446 float PHIst= P2phi(p);
451 for(
unsigned int i = 0; i < scoll.size(); ++i)
453 float y = scoll[i].y();
454 float phi= scoll[i].phi();
455 if(RD2(Yst,PHIst,y,phi) < far_def)
465 TemporaryJet jet(ptr->pT(),Yst,PHIst);
472 if(jet.is_stable(items,_CONE_RADIUS,ratio,0) && jet.is_stable(items,_SEARCH_CONE,3.0))
474 if(jet.is_stable(ilist,_CONE_RADIUS,ratio,0) && jet.is_stable(ilist,_SEARCH_CONE,3.0))
482 jet.is_stable(items,_CONE_RADIUS,ratio,0) ;
484 jet.is_stable(ilist,_CONE_RADIUS,ratio,0) ;
488 if ( _KILL_DUPLICATE )
491 float distmax = 999.;
int imax = -1;
492 for(
unsigned int i = 0; i < scoll.size(); ++i)
494 float dist = jet.dist(scoll[i]);
495 if ( dist < distmax )
501 if ( distmax > _DUPLICATE_DR ||
502 fabs((jet.pT()-scoll[imax].pT())/scoll[imax].pT())>_DUPLICATE_DPT )
504 scoll.push_back(jet);
505 mcoll.push_back(jet);
518 scoll.push_back(jet);
519 mcoll.push_back(jet);
528 for(
unsigned int i = 0; i < scoll.size(); ++i)
530 for(
unsigned int k = i+1; k < scoll.size(); ++k)
532 float djet= scoll[i].dist(scoll[k]);
533 if(djet > _CONE_RADIUS && djet < 2.*_CONE_RADIUS)
536 scoll[i].midpoint(scoll[k],y_mid,phi_mid);
537 TemporaryJet jet(-999999.,y_mid,phi_mid);
542 if(jet.is_stable(items,_CONE_RADIUS,ratio))
544 if(jet.is_stable(ilist,_CONE_RADIUS,ratio))
547 mcoll.push_back(jet);
556 ConeSplitMerge<Item> pjets(mcoll);
558 pjets.split_merge(ilcv,_SPLIT_RATIO, _PT_MIN_LEADING_PROTOJET, _PT_MIN_SECOND_PROTOJET, _MERGE_MAX, _PT_MIN_noMERGE_MAX);
561 for(
unsigned int i = 0; i < ilcv.size(); ++i)
563 if ( ilcv[i].pT() > _MIN_JET_ET )
569 std::list<const Item*> tlist=ilcv[i].LItems();
570 typename std::list<const Item*>::iterator tk;
571 for(tk = tlist.begin(); tk != tlist.end(); ++tk)
586 jets.push_back(ptrclu);
593 FASTJET_END_NAMESPACE