1 #ifndef D0RunIIconeJets_ILCONEALGORITHM 2 #define D0RunIIconeJets_ILCONEALGORITHM 89 #include "ProtoJet.hpp" 91 #include "ConeSplitMerge.hpp" 94 #include "inline_maths.h" 97 #include <fastjet/internal/base.hh> 99 FASTJET_BEGIN_NAMESPACE
103 using namespace inline_maths;
121 #ifdef ILCA_USE_ORDERED_LIST 123 template <
class Item>
127 bool operator()(
const Item* first,
const Item* second)
129 return (first->y() < second->y());
131 bool operator()(
float const & first,
const Item* second)
133 return (first < second->y());
135 bool operator()(
const Item* first,
float const& second)
137 return (first->y() < second);
144 template <
class Item>
145 class ILConeAlgorithm
158 _DUPLICATE_DR(0.005),
159 _DUPLICATE_DPT(0.01),
161 _PT_MIN_LEADING_PROTOJET(0),
162 _PT_MIN_SECOND_PROTOJET(0),
164 _PT_MIN_noMERGE_MAX(0)
168 ILConeAlgorithm(
float cone_radius,
float min_jet_Et,
float split_ratio,
169 float far_ratio=0.5,
float Et_min_ratio=0.5,
170 bool kill_duplicate=
true,
float duplicate_dR=0.005,
171 float duplicate_dPT=0.01,
float search_factor=1.0,
172 float pT_min_leading_protojet=0.,
float pT_min_second_protojet=0.,
173 int merge_max=10000,
float pT_min_nomerge=0.) :
175 _CONE_RADIUS(cone_radius),
177 _MIN_JET_ET(min_jet_Et),
179 _ET_MIN_RATIO(Et_min_ratio),
181 _FAR_RATIO(far_ratio),
183 _SPLIT_RATIO(split_ratio),
184 _DUPLICATE_DR(duplicate_dR),
185 _DUPLICATE_DPT(duplicate_dPT),
186 _SEARCH_CONE(cone_radius/search_factor),
189 _KILL_DUPLICATE(kill_duplicate),
190 _PT_MIN_LEADING_PROTOJET(pT_min_leading_protojet),
191 _PT_MIN_SECOND_PROTOJET(pT_min_second_protojet),
192 _MERGE_MAX(merge_max),
193 _PT_MIN_noMERGE_MAX(pT_min_nomerge)
197 ~ILConeAlgorithm() {;}
203 std::list<Item> &jets,
204 std::list<const Item*>& itemlist,
209 const float Item_ET_Threshold);
212 std::vector<ProtoJet<Item> > ilcv;
222 float _DUPLICATE_DPT;
225 bool _KILL_DUPLICATE;
227 float _PT_MIN_LEADING_PROTOJET;
228 float _PT_MIN_SECOND_PROTOJET;
230 float _PT_MIN_noMERGE_MAX;
235 class TemporaryJet :
public ProtoJet<Item>
240 TemporaryJet(
float seedET) : ProtoJet<Item>(seedET) {;}
242 TemporaryJet(
float seedET,
float y_in,
float phi_in) :
243 ProtoJet<Item>(seedET,y_in,phi_in) {;}
247 float dist(TemporaryJet& jet)
const 249 return RDelta(this->_y,this->_phi,jet.y(),jet.phi());
252 void midpoint(
const TemporaryJet& jet,
float & y_out,
float & phi_out)
const 256 float pTsum = this->_pT + jet.pT();
257 y_out = (this->_y*this->_pT + jet.y()*jet.pT())/pTsum;
259 phi_out = (this->_phi*this->_pT + jet.phi()*jet.pT())/pTsum;
264 if ( fabs(phi_out-this->_phi)>2.0 ) {
265 phi_out = fmod( this->_phi+PI, TWOPI);
266 if (phi_out < 0.0) phi_out += TWOPI;
269 float temp=fmod( jet.phi()+PI, TWOPI);
270 if (temp < 0.0) temp += TWOPI;
273 phi_out = (phi_out*this->_pT + temp*jet.pT()) /pTsum;
276 if ( phi_out < 0. ) phi_out += TWOPI;
282 bool is_stable(
const std::multimap<float,const Item*>& items,
283 float radius,
float min_ET,
int max_iterations=50)
285 bool is_stable(
const std::list<const Item*>& itemlist,
float radius,
286 float min_ET,
int max_iterations=50)
290 float radius2 = radius*radius;
307 this->setJet(Yst,PHIst,0.0);
310 #ifdef ILCA_USE_ORDERED_LIST 311 std::list<const Item*>::const_iterator lower =
312 lower_bound(itemlist.begin(),itemlist.end(),Yst-radius,
313 rapidity_order<Item>());
314 std::list<const Item*>::const_iterator upper =
315 upper_bound(itemlist.begin(),itemlist.end(),Yst+radius,
316 rapidity_order<Item>());
317 for(std::list<const Item*>::const_iterator tk = lower; tk != upper; ++tk) {
319 if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
327 for ( std::multimap<float,const Item*>::const_iterator
328 tk = items.lower_bound(Yst-radius);
329 tk != items.upper_bound(Yst+radius); ++tk )
332 if(RD2(((*tk).second)->y(),((*tk).second)->phi(),Yst,PHIst) <= radius2)
334 this->addItem((*tk).second);
341 for(
typename std::list<const Item*>::const_iterator tk = itemlist.begin(); tk != itemlist.end(); ++tk)
344 if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
357 if(this->_pT < min_ET )
363 }
while(RD2(this->_y,this->_phi,Yst,PHIst) >= Rcut && trial <= max_iterations);
375 template <
class Item>
376 void ILConeAlgorithm <Item>::
379 std::list<Item> &jets,
380 std::list<const Item*>& ilist,
385 const float Item_ET_Threshold)
388 for (
typename std::list<const Item*>::iterator it = ilist.begin();
392 if ( (*it)->pT() < Item_ET_Threshold )
394 it = ilist.erase(it);
404 std::vector<const Item*> ecv;
405 for (
typename std::list<const Item*>::iterator it = ilist.begin();
406 it != ilist.end(); it++) {
417 float far_def = _FAR_RATIO*_CONE_RADIUS * _FAR_RATIO*_CONE_RADIUS;
420 float ratio= _MIN_JET_ET*_ET_MIN_RATIO;
423 #ifdef ILCA_USE_ORDERED_LIST 425 ilist.sort(rapidity_order<Item>());
429 std::multimap<float,const Item*> items;
430 std::list<const Item*>::const_iterator it;
431 for(it = ilist.begin(); it != ilist.end(); ++it)
433 pair<float,const Item*> p((*it)->y(),*it);
439 std::vector<ProtoJet<Item> > mcoll;
440 std::vector<TemporaryJet> scoll;
445 typename std::vector<const Item*>::iterator jclu;
446 for(jclu = ecv.begin(); jclu != ecv.end(); ++jclu)
449 const Item* ptr= *jclu;
453 float PHIst= P2phi(p);
458 for(
unsigned int i = 0; i < scoll.size(); ++i)
460 float y = scoll[i].y();
461 float phi= scoll[i].phi();
462 if(RD2(Yst,PHIst,y,phi) < far_def)
472 TemporaryJet jet(ptr->pT(),Yst,PHIst);
479 if(jet.is_stable(items,_CONE_RADIUS,ratio,0) && jet.is_stable(items,_SEARCH_CONE,3.0))
481 if(jet.is_stable(ilist,_CONE_RADIUS,ratio,0) && jet.is_stable(ilist,_SEARCH_CONE,3.0))
489 jet.is_stable(items,_CONE_RADIUS,ratio,0) ;
491 jet.is_stable(ilist,_CONE_RADIUS,ratio,0) ;
495 if ( _KILL_DUPLICATE )
498 float distmax = 999.;
int imax = -1;
499 for(
unsigned int i = 0; i < scoll.size(); ++i)
501 float dist = jet.dist(scoll[i]);
502 if ( dist < distmax )
508 if ( distmax > _DUPLICATE_DR ||
509 fabs((jet.pT()-scoll[imax].pT())/scoll[imax].pT())>_DUPLICATE_DPT )
511 scoll.push_back(jet);
512 mcoll.push_back(jet);
525 scoll.push_back(jet);
526 mcoll.push_back(jet);
535 for(
unsigned int i = 0; i < scoll.size(); ++i)
537 for(
unsigned int k = i+1; k < scoll.size(); ++k)
539 float djet= scoll[i].dist(scoll[k]);
540 if(djet > _CONE_RADIUS && djet < 2.*_CONE_RADIUS)
543 scoll[i].midpoint(scoll[k],y_mid,phi_mid);
544 TemporaryJet jet(-999999.,y_mid,phi_mid);
549 if(jet.is_stable(items,_CONE_RADIUS,ratio))
551 if(jet.is_stable(ilist,_CONE_RADIUS,ratio))
554 mcoll.push_back(jet);
563 ConeSplitMerge<Item> pjets(mcoll);
565 pjets.split_merge(ilcv,_SPLIT_RATIO, _PT_MIN_LEADING_PROTOJET, _PT_MIN_SECOND_PROTOJET, _MERGE_MAX, _PT_MIN_noMERGE_MAX);
568 for(
unsigned int i = 0; i < ilcv.size(); ++i)
570 if ( ilcv[i].pT() > _MIN_JET_ET )
576 std::list<const Item*> tlist=ilcv[i].LItems();
577 typename std::list<const Item*>::iterator tk;
578 for(tk = tlist.begin(); tk != tlist.end(); ++tk)
593 jets.push_back(ptrclu);
600 FASTJET_END_NAMESPACE