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