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