FastJet  3.0.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ILConeAlgorithm.hpp
1 #ifndef D0RunIIconeJets_ILCONEALGORITHM
2 #define D0RunIIconeJets_ILCONEALGORITHM
3 // ---------------------------------------------------------------------------
4 // ILConeAlgorithm.hpp
5 //
6 // Created: 28-JUL-2000 Francois Touze (+ Laurent Duflot)
7 //
8 // Purpose: Implements the Improved Legacy Cone Algorithm
9 //
10 // Modified:
11 // 31-JUL-2000 Laurent Duflot
12 // + introduce support for additional informations (ConeJetInfo)
13 // 1-AUG-2000 Laurent Duflot
14 // + seedET for midpoint jet is -999999., consistent with seedET ordering
15 // in ConeSplitMerge: final jets with seedET=-999999. will be midpoint
16 // jets which were actually different from stable cones.
17 // 4-Aug-2000 Laurent Duflot
18 // + remove unecessary copy in TemporaryJet::is_stable()
19 // 11-Aug-2000 Laurent Duflot
20 // + remove using namespace std
21 // + add threshold on item. The input list to makeClusters() IS modified
22 // 20-June-2002 John Krane
23 // + remove bug in midpoint calculation based on pT weight
24 // + started to add new midpoint calculation based on 4-vectors,
25 // but not enough info held by this class
26 // 24-June-2002 John Krane
27 // + modify is_stable() to not iterate if desired
28 // (to expand search cone w/out moving it)
29 // + added search cone logic for initial jets but not midpoint jets as per
30 // agreement with CDF
31 // 19-July-2002 John Krane
32 // + _SEARCH_CONE size factor now provided by calreco/CalClusterReco.cpp
33 // + (default = 1.0 = like original ILCone behavior)
34 // 10-Oct-2002 John Krane
35 // + Check min Pt of cone with full size first, then iterate with search cone
36 // 07-Dec-2002 John Krane
37 // + speed up the midpoint phi-wrap solution
38 // 01-May-2007 Lars Sonnenschein
39 // extracted from D0 software framework and modified to remove subsequent dependencies
40 //
41 //
42 // This file is distributed with FastJet under the terms of the GNU
43 // General Public License (v2). Permission to do so has been granted
44 // by Lars Sonnenschein and the D0 collaboration (see COPYING for
45 // details)
46 //
47 // History of changes in FastJet compared tothe original version of
48 // ProtoJet.hpp
49 //
50 // 2011-12-13 Gregory Soyez <soyez@fastjet.fr>
51 //
52 // * added license information
53 //
54 // 2011-11-14 Gregory Soyez <soyez@fastjet.fr>
55 //
56 // * changed the name of a few parameters to avoid a gcc
57 // -Wshadow warning
58 //
59 // 2009-01-17 Gregory Soyez <soyez@fastjet.fr>
60 //
61 // * put the code in the fastjet::d0 namespace
62 //
63 // 2007-12-14 Gavin Salam <salam@lpthe.jussieu.fr>
64 //
65 // * moved the 'std::vector<ProtoJet<Item> > ilcv' structure
66 // containing the info about the final jets from a local
67 // variable to a class variable (for integration in the
68 // FastJet plugin core)
69 //
70 // ---------------------------------------------------------------------------
71 
72 #include <vector>
73 #include <list>
74 #include <utility> // defines pair<,>
75 #include <map>
76 #include <algorithm>
77 #include <iostream>
78 
79 
80 //#include "energycluster/EnergyClusterReco.hpp"
81 //#include "energycluster/ProtoJet.hpp"
82 #include "ProtoJet.hpp"
83 //#include "energycluster/ConeSplitMerge.hpp"
84 #include "ConeSplitMerge.hpp"
85 //#include "energycluster/ConeJetInfoChunk.hpp"
86 
87 #include "inline_maths.h"
88 
89 ///////////////////////////////////////////////////////////////////////////////
90 #include <fastjet/internal/base.hh>
91 
92 FASTJET_BEGIN_NAMESPACE
93 
94 namespace d0{
95 
96 using namespace inline_maths;
97 
98 /*
99  NB: Some attempt at optimizing the code has been made by ordering the object
100  along rapidity but this did not improve the speed. There are traces of
101  these atemps in the code, that will be cleaned up in the future.
102  */
103 
104 // at most one of those !
105 // order the input list and use lower_bound() and upper_bound() to restrict the
106 // on item to those that could possibly be in the cone.
107 //#define ILCA_USE_ORDERED_LIST
108 
109 // idem but use an intermediate multimap in hope that lower_bound() and
110 // upper_bound() are faster in this case.
111 //#define ILCA_USE_MMAP
112 
113 
114 #ifdef ILCA_USE_ORDERED_LIST
115 // this class is used to order the item list along rapidity
116 template <class Item>
117 class rapidity_order
118 {
119 public:
120  bool operator()(const Item* first, const Item* second)
121  {
122  return (first->y() < second->y());
123  }
124  bool operator()(float const & first, const Item* second)
125  {
126  return (first < second->y());
127  }
128  bool operator()(const Item* first, float const& second)
129  {
130  return (first->y() < second);
131  }
132 };
133 #endif
134 
135 
136 //template <class Item,class ItemAddress,class IChunk>
137 template <class Item>
138 class ILConeAlgorithm
139 {
140 
141 public:
142 
143  // default constructor (default parameters are crazy: you should not use that
144  // constructor !)
145  ILConeAlgorithm():
146  _CONE_RADIUS(0.),
147  _MIN_JET_ET(99999.),
148  _ET_MIN_RATIO(0.5),
149  _FAR_RATIO(0.5),
150  _SPLIT_RATIO(0.5),
151  _DUPLICATE_DR(0.005),
152  _DUPLICATE_DPT(0.01),
153  _SEARCH_CONE(0.5),
154  _PT_MIN_LEADING_PROTOJET(0),
155  _PT_MIN_SECOND_PROTOJET(0),
156  _MERGE_MAX(10000),
157  _PT_MIN_noMERGE_MAX(0)
158  {;}
159 
160  // full constructor
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.) :
167  // cone radius
168  _CONE_RADIUS(cone_radius),
169  // minimum jet ET
170  _MIN_JET_ET(min_jet_Et),
171  // stable cones must have ET > _ET_MIN_RATIO*_MIN_JET_ET at any iteration
172  _ET_MIN_RATIO(Et_min_ratio),
173  // precluster at least _FAR_RATIO*_CONE_RADIUS away from stable cones
174  _FAR_RATIO(far_ratio),
175  // split or merge criterium
176  _SPLIT_RATIO(split_ratio),
177  _DUPLICATE_DR(duplicate_dR),
178  _DUPLICATE_DPT(duplicate_dPT),
179  _SEARCH_CONE(cone_radius/search_factor),
180  // kill stable cone if within _DUPLICATE_DR and delta(pT)<_DUPLICATE_DPT
181  // of another stable cone.
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)
187  {;}
188 
189  //destructor
190  ~ILConeAlgorithm() {;}
191 
192  // make jet clusters using improved legacy cone algorithm
193  //void makeClusters(const EnergyClusterReco* r,
194  void makeClusters(
195  // the input list modified (ordered)
196  std::list<Item> &jets,
197  std::list<const Item*>& itemlist,
198  //float zvertex,
199  ////std::list<const Item*>& itemlist);
200  //const EnergyClusterCollection<ItemAddress>& preclu,
201  //IChunk* chunkptr, ConeJetInfoChunk* infochunkptr,
202  const float Item_ET_Threshold);
203 
204  // this will hold the final jets + contents
205  std::vector<ProtoJet<Item> > ilcv;
206 
207 private:
208 
209  float _CONE_RADIUS;
210  float _MIN_JET_ET;
211  float _ET_MIN_RATIO;
212  float _FAR_RATIO;
213  float _SPLIT_RATIO;
214  float _DUPLICATE_DR;
215  float _DUPLICATE_DPT;
216  float _SEARCH_CONE;
217 
218  bool _KILL_DUPLICATE;
219 
220  float _PT_MIN_LEADING_PROTOJET;
221  float _PT_MIN_SECOND_PROTOJET;
222  int _MERGE_MAX;
223  float _PT_MIN_noMERGE_MAX;
224 
225  // private class
226  // This is a ProtoJet with additional methods dist(), midpoint() and
227  // is_stable()
228  class TemporaryJet : public ProtoJet<Item>
229  {
230 
231  public :
232 
233  TemporaryJet(float seedET) : ProtoJet<Item>(seedET) {;}
234 
235  TemporaryJet(float seedET,float y_in,float phi_in) :
236  ProtoJet<Item>(seedET,y_in,phi_in) {;}
237 
238  ~TemporaryJet() {;}
239 
240  float dist(TemporaryJet& jet) const
241  {
242  return RDelta(this->_y,this->_phi,jet.y(),jet.phi());
243  }
244 
245  void midpoint(const TemporaryJet& jet,float & y_out, float & phi_out) const
246  {
247  // Midpoint should probably be computed w/4-vectors but don't
248  // have that info. Preserving Pt-weighted calculation - JPK
249  float pTsum = this->_pT + jet.pT();
250  y_out = (this->_y*this->_pT + jet.y()*jet.pT())/pTsum;
251 
252  phi_out = (this->_phi*this->_pT + jet.phi()*jet.pT())/pTsum;
253  // careful with phi-wrap area: convert from [0,2pi] to [-pi,pi]
254  //ls: original D0 code, as of 23/Mar/2007
255  //if ( abs(phi-this->_phi)>2.0 ) { // assumes cones R=1.14 or smaller, merge within 2R only
256  //ls: abs bug fixed 26/Mar/2007
257  if ( fabs(phi_out-this->_phi)>2.0 ) { // assumes cones R=1.14 or smaller, merge within 2R only
258  phi_out = fmod( this->_phi+PI, TWOPI);
259  if (phi_out < 0.0) phi_out += TWOPI;
260  phi_out -= PI;
261 
262  float temp=fmod( jet.phi()+PI, TWOPI);
263  if (temp < 0.0) temp += TWOPI;
264  temp -= PI;
265 
266  phi_out = (phi_out*this->_pT + temp*jet.pT()) /pTsum;
267  }
268 
269  if ( phi_out < 0. ) phi_out += TWOPI;
270  }
271 
272 
273 ////////////////////////////////////////
274 #ifdef ILCA_USE_MMAP
275  bool is_stable(const std::multimap<float,const Item*>& items,
276  float radius, float min_ET, int max_iterations=50)
277 #else
278  bool is_stable(const std::list<const Item*>& itemlist, float radius,
279  float min_ET, int max_iterations=50)
280 #endif
281  // Note: max_iterations = 0 will just recompute the jet using the specified cone
282  {
283  float radius2 = radius*radius;
284  float Rcut= 1.E-06;
285 
286 
287  // ?? if(_Increase_Delta_R) Rcut= 1.E-04;
288  bool stable= true;
289  int trial= 0;
290  float Yst;
291  float PHIst;
292  do {
293  trial++;
294  //std::cout << " trial " << trial << " " << _y << " " << _phi << std::endl;
295  Yst = this->_y;
296  PHIst= this->_phi;
297  //cout << "is_stable beginning do loop: this->_pT=" << this->_pT << " this->_y=" << this->_y << " this->_phi=" << this->_phi << endl;
298  this->erase();
299 
300  this->setJet(Yst,PHIst,0.0);
301 
302 
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) {
311  //std::cout << " is_stable: item y=" << (*tk)->y() << " phi=" << (*tk)->phi() << " RD2=" << RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) << " " << Yst-radius << " " << Yst+radius << endl;
312  if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
313  {
314  addItem(*tk);
315  }
316  }
317 #else
318 #ifdef ILCA_USE_MMAP
319  // need to loop only on the subset with Yst-R < y < Yst+R
320  for ( std::multimap<float,const Item*>::const_iterator
321  tk = items.lower_bound(Yst-radius);
322  tk != items.upper_bound(Yst+radius); ++tk )
323  {
324  //std::cout << " item " << (*tk)->y() << " " << (*tk)->phi() << " " << RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) << " " << Yst-radius << " " << Yst+radius << endl;
325  if(RD2(((*tk).second)->y(),((*tk).second)->phi(),Yst,PHIst) <= radius2)
326  {
327  addItem((*tk).second);
328  }
329  }
330 
331 #else
332 
333  //cout << " is_stable: itemlist.size()=" << itemlist.size() << endl;
334  for(typename std::list<const Item*>::const_iterator tk = itemlist.begin(); tk != itemlist.end(); ++tk)
335  {
336  //std::cout << " is_stable: item (*tk)->y()=" << (*tk)->y() << " (*tk)->phi=" << (*tk)->phi() << " RD2=" << RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) << " Yst-rad=" << Yst-radius << " Yst+rad=" << Yst+radius << endl;
337  if(RD2((*tk)->y(),(*tk)->phi(),Yst,PHIst) <= radius2)
338  {
339  //cout << "add item to *tk" << endl;
340  this->addItem(*tk);
341  }
342  }
343 #endif
344 #endif
345 
346  //std::cout << "is_stable, before update: jet this->_y=" << this->_y << " _phi=" << this->_phi << " _pT=" << this->_pT << " min_ET=" << min_ET << std::endl;
347  this->updateJet();
348  //std::cout << "is_stable, after update: jet this->_y=" << this->_y << " _phi=" << this->_phi << " _pT=" << this->_pT << " min_ET=" << min_ET << std::endl;
349 
350  if(this->_pT < min_ET )
351  {
352  stable= false;
353  break;
354  }
355  //cout << "is_stable end while loop: this->_pT=" << this->_pT << " this->_y=" << this->_y << " this->_phi=" << this->_phi << endl;
356  } while(RD2(this->_y,this->_phi,Yst,PHIst) >= Rcut && trial <= max_iterations);
357  //std::cout << " trial stable " << trial << " " << stable << std::endl;
358  return stable;
359  }
360 
361  private :
362 
363  };
364 };
365 ///////////////////////////////////////////////////////////////////////////////
366 //template <class Item,class ItemAddress,class IChunk>
367 //void ILConeAlgorithm <Item,ItemAddress,IChunk>::
368 template <class Item>
369 void ILConeAlgorithm <Item>::
370 //makeClusters(const EnergyClusterReco* r,
371 makeClusters(
372  std::list<Item> &jets,
373  std::list<const Item*>& ilist,
374  //float Zvertex,
375  ////std::list<const Item*>& ilist)
376  //const EnergyClusterCollection<ItemAddress>& preclu,
377  //IChunk* chunkptr, ConeJetInfoChunk* infochunkptr,
378  const float Item_ET_Threshold)
379 {
380  // remove items below threshold
381  for ( typename std::list<const Item*>::iterator it = ilist.begin();
382 
383  it != ilist.end(); )
384  {
385  if ( (*it)->pT() < Item_ET_Threshold )
386  {
387  it = ilist.erase(it);
388  }
389  else ++it;
390  }
391 
392  // create an energy cluster collection for jets
393  //EnergyClusterCollection<ItemAddress>* ptrcol;
394  //Item* ptrcol;
395  //r->createClusterCollection(chunkptr,ptrcol);
396  ////std::vector<const EnergyCluster<ItemAddress>*> ecv;
397  std::vector<const Item*> ecv;
398  for ( typename std::list<const Item*>::iterator it = ilist.begin();
399  it != ilist.end(); it++) {
400  ecv.push_back(*it);
401  }
402 
403 
404  //preclu.getClusters(ecv);
405  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% need to fill here vector ecv
406 
407  //std::cout << " number of seed clusters: " << ecv.size() << std::endl;
408 
409  // skip precluster close to jets
410  float far_def = _FAR_RATIO*_CONE_RADIUS * _FAR_RATIO*_CONE_RADIUS;
411 
412  // skip if jet Et is below some value
413  float ratio= _MIN_JET_ET*_ET_MIN_RATIO;
414 
415 
416 #ifdef ILCA_USE_ORDERED_LIST
417  // sort the list in rapidity order
418  ilist.sort(rapidity_order<Item>());
419 #else
420 #ifdef ILCA_USE_MMAP
421  // create a y ordered list of items
422  std::multimap<float,const Item*> items;
423  std::list<const Item*>::const_iterator it;
424  for(it = ilist.begin(); it != ilist.end(); ++it)
425  {
426  pair<float,const Item*> p((*it)->y(),*it);
427  items.insert(p);
428  }
429 #endif
430 #endif
431 
432  std::vector<ProtoJet<Item> > mcoll;
433  std::vector<TemporaryJet> scoll;
434 
435 
436  // find stable jets around seeds
437  //typename std::vector<const EnergyCluster<ItemAddress>* >::iterator jclu;
438  typename std::vector<const Item*>::iterator jclu;
439  for(jclu = ecv.begin(); jclu != ecv.end(); ++jclu)
440  {
441  //const EnergyCluster<ItemAddress>* ptr= *jclu;
442  const Item* ptr= *jclu;
443  float p[4];
444  ptr->p4vec(p);
445  float Yst = P2y(p);
446  float PHIst= P2phi(p);
447 
448  // don't keep preclusters close to jet
449  bool is_far= true;
450  // ?? if(_Kill_Far_Clusters) {
451  for(unsigned int i = 0; i < scoll.size(); ++i)
452  {
453  float y = scoll[i].y();
454  float phi= scoll[i].phi();
455  if(RD2(Yst,PHIst,y,phi) < far_def)
456  {
457  is_far= false;
458  break;
459  }
460  }
461  // ?? }
462 
463  if(is_far)
464  {
465  TemporaryJet jet(ptr->pT(),Yst,PHIst);
466  //cout << "temporary jet: pt=" << ptr->pT() << " y=" << Yst << " phi=" << PHIst << endl;
467 
468  // Search cones are smaller, so contain less jet Et
469  // Don't throw out too many little jets during search phase!
470  // Strategy: check Et first with full cone, then search with low-GeV min_et thresh
471 #ifdef ILCA_USE_MMAP
472  if(jet.is_stable(items,_CONE_RADIUS,ratio,0) && jet.is_stable(items,_SEARCH_CONE,3.0))
473 #else
474  if(jet.is_stable(ilist,_CONE_RADIUS,ratio,0) && jet.is_stable(ilist,_SEARCH_CONE,3.0))
475 #endif
476  {
477 
478  //cout << "temporary jet is stable" << endl;
479 
480 // jpk Resize the found jets
481 #ifdef ILCA_USE_MMAP
482  jet.is_stable(items,_CONE_RADIUS,ratio,0) ;
483 #else
484  jet.is_stable(ilist,_CONE_RADIUS,ratio,0) ;
485 #endif
486  //cout << "found jet has been resized if any" << endl;
487 
488  if ( _KILL_DUPLICATE )
489  {
490  // check if we are not finding the same jet again
491  float distmax = 999.; int imax = -1;
492  for(unsigned int i = 0; i < scoll.size(); ++i)
493  {
494  float dist = jet.dist(scoll[i]);
495  if ( dist < distmax )
496  {
497  distmax = dist;
498  imax = i;
499  }
500  }
501  if ( distmax > _DUPLICATE_DR ||
502  fabs((jet.pT()-scoll[imax].pT())/scoll[imax].pT())>_DUPLICATE_DPT )
503  {
504  scoll.push_back(jet);
505  mcoll.push_back(jet);
506  //std::cout << " stable cone " << jet.y() << " " << jet.phi() << " " << jet.pT() << std::endl;
507  }
508  /*
509  else
510  {
511  std::cout << " jet too close to a found jet " << distmax << " " <<
512  fabs((jet.pT()-scoll[imax].pT())/scoll[imax].pT()) << std::endl;
513  }
514  */
515  }
516  else
517  {
518  scoll.push_back(jet);
519  mcoll.push_back(jet);
520  //std::cout << " stable cone " << jet.y() << " " << jet.phi() << " " << jet.pT() << std::endl;
521  }
522 
523  }
524  }
525  }
526 
527  // find stable jets around midpoints
528  for(unsigned int i = 0; i < scoll.size(); ++i)
529  {
530  for(unsigned int k = i+1; k < scoll.size(); ++k)
531  {
532  float djet= scoll[i].dist(scoll[k]);
533  if(djet > _CONE_RADIUS && djet < 2.*_CONE_RADIUS)
534  {
535  float y_mid,phi_mid;
536  scoll[i].midpoint(scoll[k],y_mid,phi_mid);
537  TemporaryJet jet(-999999.,y_mid,phi_mid);
538  //std::cout << " midpoint: " << scoll[i].pT() << " " << scoll[i].info().seedET() << " " << scoll[k].pT() << " " << scoll[k].info().seedET() << " " << y_mid << " " << phi_mid << std::endl;
539 
540 // midpoint jets are full size
541 #ifdef ILCA_USE_MMAP
542  if(jet.is_stable(items,_CONE_RADIUS,ratio))
543 #else
544  if(jet.is_stable(ilist,_CONE_RADIUS,ratio))
545 #endif
546  {
547  mcoll.push_back(jet);
548  //std::cout << " stable midpoint cone " << jet.y() << " " << jet.phi() << " " << jet.pT() << std::endl;
549  }
550  }
551  }
552  }
553 
554 
555  // do a pT ordered splitting/merging
556  ConeSplitMerge<Item> pjets(mcoll);
557  ilcv.clear();
558  pjets.split_merge(ilcv,_SPLIT_RATIO, _PT_MIN_LEADING_PROTOJET, _PT_MIN_SECOND_PROTOJET, _MERGE_MAX, _PT_MIN_noMERGE_MAX);
559 
560 
561  for(unsigned int i = 0; i < ilcv.size(); ++i)
562  {
563  if ( ilcv[i].pT() > _MIN_JET_ET )
564  {
565  //EnergyCluster<ItemAddress>* ptrclu;
566  Item ptrclu;
567  //r->createCluster(ptrcol,ptrclu);
568 
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)
572  {
573  //ItemAddress addr= (*tk)->address();
574  float pk[4];
575  (*tk)->p4vec(pk);
576  //std::cout << (*tk)->index <<" " << (*tk) << std::endl;
577  //float emE= (*tk)->emE();
578  //r->addClusterItem(ptrclu,addr,pk,emE);
579  //ptrclu->Add(*tk);
580  ptrclu.Add(**tk);
581  }
582  // print out
583  //ptrclu->print(cout);
584 
585  //infochunkptr->addInfo(ilcv[i].info());
586  jets.push_back(ptrclu);
587  }
588  }
589 }
590 
591 } // namespace d0
592 
593 FASTJET_END_NAMESPACE
594 
595 #endif
596 
597