FastJet  3.2.1
ConeSplitMerge.hpp
1 #ifndef D0RunIIconeJets_CONESPLITMERGE
2 #define D0RunIIconeJets_CONESPLITMERGE
3 // ---------------------------------------------------------------------------
4 // ConeSplitMerge.hpp
5 //
6 // Created: 28-JUL-2000 Francois Touze
7 //
8 // Purpose: Implements the pT ordered jet split/merge algorithm for the
9 // Improved Legacy Cone Algorithm split/merge algo.
10 //
11 // Modified:
12 // 31-JUL-2000 Laurent Duflot
13 // + introduce support for additional informations (ConeJetInfo)
14 // 1-AUG-2000 Laurent Duflot
15 // + jet merge criteria was wrong, now calculate shared_ET and compare to
16 // neighbour jet pT.
17 // + split was incomplete: shared items not really removed from jets.
18 // 4-Aug-2000 Laurent Duflot
19 // + use item methods y() and phi() rather than p4vec() and then P2y and
20 // P2phi
21 // 7-Aug-2000 Laurent Duflot
22 // + force the list to be organized by decreasing ET and, for identical ET,
23 // by decreasing seedET. Identical protojets can be found by eg nearby
24 // seeds. The seedET ordering is such that for identical jets, the one
25 // with the highest seedET is kept, which is what we want for efficiency
26 // calculations.
27 // + avoid unecessary copies of lists by using reference when possible
28 // 9-Aug-2000 Laurent Duflot
29 // + save initial jet ET before split/merge
30 // 1-May-2007 Lars Sonnenschein
31 // extracted from D0 software framework and modified to remove subsequent dependencies
32 //
33 //
34 // This file is distributed with FastJet under the terms of the GNU
35 // General Public License (v2). Permission to do so has been granted
36 // by Lars Sonnenschein and the D0 collaboration (see COPYING for
37 // details)
38 //
39 // History of changes in FastJet compared tothe original version of
40 // ConeSplitMerge.hpp
41 //
42 // 2011-12-13 Gregory Soyez <soyez@fastjet.fr>
43 //
44 // * added license information
45 //
46 // 2009-01-17 Gregory Soyez <soyez@fastjet.fr>
47 //
48 // * put the code in the fastjet::d0 namespace
49 //
50 // 2007-12-14 Gavin Salam <salam@lpthe.jussieu.fr>
51 //
52 // * replaced make_pair by std::make_pair
53 //
54 // ---------------------------------------------------------------------------
55 
56 
57 #include <iostream>
58 #include <map>
59 #include <utility>
60 #include <vector>
61 #include "ProtoJet.hpp"
62 
63 //using namespace D0RunIIconeJets_CONEJETINFO;
64 
65 #include <fastjet/internal/base.hh>
66 
67 FASTJET_BEGIN_NAMESPACE
68 
69 namespace d0{
70 
71 //
72 // this class is used to order ProtoJets by decreasing ET and seed ET
73 template <class Item>
74 class ProtoJet_ET_seedET_order
75 {
76 public:
77  bool operator()(const ProtoJet<Item> & first, const ProtoJet<Item> & second) const
78  {
79  if ( first.pT() > second.pT() ) return true;
80  else
81  if ( first.pT() < second.pT() ) return false;
82  else return ( first.info().seedET() > second.info().seedET() );
83  }
84 };
85 
86 
87 template <class Item>
88 class ConeSplitMerge {
89 
90 public :
91 
92  typedef std::multimap<ProtoJet<Item>,float,ProtoJet_ET_seedET_order<Item> > PJMMAP;
93 
94  ConeSplitMerge();
95  ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector);
96  ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist);
97  ~ConeSplitMerge() {;}
98  void split_merge(std::vector<ProtoJet<Item> >& ecv,float s, float pT_min_leading_protojet, float pT_min_second_protojet, int MergeMax, float pT_min_noMergeMax);
99 
100 private :
101  PJMMAP _members;
102 };
103 ///////////////////////////////////////////////////////////////////////////////
104 template<class Item>
105 ConeSplitMerge<Item>::ConeSplitMerge():_members() {;}
106 
107 template<class Item>
108 ConeSplitMerge<Item>::ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector)
109 {
110  // sort proto_jets in Et descending order
111  typename std::vector<ProtoJet<Item> >::const_iterator jt;
112  for(jt = jvector.begin(); jt != jvector.end(); ++jt)
113  {
114  // this is supposed to be a stable cone, declare so
115  ProtoJet<Item> jet(*jt);
116  jet.NowStable();
117  _members.insert(std::make_pair(jet,jet.pT()));
118  }
119 }
120 
121 template<class Item>
122 ConeSplitMerge<Item>::ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist)
123 {
124  //_max_nb_items =-1;
125  // sort proto_jets in Et descending order
126  typename std::list<ProtoJet<Item> >::const_iterator jt;
127  for(jt = jlist.begin(); jt != jlist.end(); ++jt)
128  {
129  // this is supposed to be a stable cone, declare so
130  ProtoJet<Item> jet(*jt);
131  jet.NowStable();
132  _members.insert(std::make_pair(jet,jet.pT()));
133  }
134 }
135 
136 template<class Item>
137 void ConeSplitMerge<Item>::split_merge(std::vector<ProtoJet<Item> >& jcv,
138  float shared_ET_fraction,
139  float pT_min_leading_protojet, float pT_min_second_protojet,
140  int MergeMax, float pT_min_noMergeMax)
141 {
142  while(!_members.empty())
143  {
144  /*
145  {
146  std::cout << std::endl;
147  std::cout << " --------------- list of protojets ------------------ " <<std::endl;
148  for ( PJMMAP::iterator it = _members.begin();
149  it != _members.end(); ++it)
150  {
151  std::cout << " pT y phi " << (*it).pT() << " " << (*it).y() << " " << (*it).phi() << " " << (*it).info().seedET() << " " << (*it).info().nbMerge() << " " << (*it).info().nbSplit() << std::endl;
152  }
153  std::cout << " ----------------------------------------------------- " <<std::endl;
154  }
155  */
156 
157 
158  // select highest Et jet
159  typename PJMMAP::iterator itmax= _members.begin();
160  ProtoJet<Item> imax((*itmax).first);
161  const std::list<const Item*>& ilist(imax.LItems());
162 
163  _members.erase(itmax);
164 
165  // does jet share items?
166  bool share= false;
167  float shared_ET = 0.;
168  typename PJMMAP::iterator jtmax;
169  typename PJMMAP::iterator jt;
170  for(jt = _members.begin(); jt != _members.end(); ++jt)
171  {
172  const std::list<const Item*>& jlist((*jt).first.LItems());
173  typename std::list<const Item*>::const_iterator tk;
174  int count;
175  for(tk = ilist.begin(), count = 0; tk != ilist.end();
176  ++tk, ++count)
177  {
178  typename std::list<const Item*>::const_iterator where=
179  find(jlist.begin(),jlist.end(),*tk);
180  if(where != jlist.end())
181  {
182  share= true;
183  shared_ET += (*tk)->pT();
184  }
185  }
186  if(share)
187  {
188  jtmax = jt;
189  break;
190  }
191  }
192 
193  if(!share) {
194  // add jet to the final jet list
195  jcv.push_back(imax);
196  //std::cout << " final jet " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl;
197  }
198  else
199  {
200  // find highest Et neighbor
201  ProtoJet<Item> jmax((*jtmax).first);
202 
203  // drop neighbor
204  _members.erase(jtmax);
205 
206 
207  //std::cout << " split or merge ? " << imax.pT() << " " << shared_ET << " " << jmax.pT() << " " << s << " " << (jmax.pT())*s << std::endl;
208 
209  // merge
210  if( shared_ET > (jmax.pT())*shared_ET_fraction
211  && (imax.pT() > pT_min_leading_protojet || jmax.pT() > pT_min_second_protojet)
212  && (imax.info().nbMerge() < MergeMax || imax.pT() > pT_min_noMergeMax))
213  {
214  // add neighbor's items to imax
215  const std::list<const Item*>& jlist(jmax.LItems());
216  typename std::list<const Item*>::const_iterator tk;
217  typename std::list<const Item*>::const_iterator iend= ilist.end();
218  bool same = true; // is jmax just the same as imax ?
219  for(tk = jlist.begin(); tk != jlist.end(); ++tk)
220  {
221  typename std::list<const Item*>::const_iterator where=
222  find(ilist.begin(),iend,*tk);
223  if(where == iend)
224  {
225  imax.addItem(*tk);
226  same = false;
227  }
228  }
229  if ( ! same )
230  {
231  // recalculate
232  //float old_pT = imax.pT();
233 
234  imax.updateJet();
235  imax.merged();
236  //std::cout << " jet merge :: " << old_pT << " " << jmax.pT() << " " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl;
237  }
238  }
239 
240  //split and assign removed shared cells from lowest pT protojet
241  else if(shared_ET > (jmax.pT())*shared_ET_fraction)
242  {
243  // here we need to pull the lists, because there are items to remove
244  std::list<const Item*> ilist(imax.LItems());
245  std::list<const Item*> jlist(jmax.LItems());
246 
247  typename std::list<const Item*>::iterator tk;
248  for(tk = jlist.begin(); tk != jlist.end(); )
249  {
250  typename std::list<const Item*>::iterator where=
251  find(ilist.begin(),ilist.end(),*tk);
252  if(where != ilist.end()) {
253  tk = jlist.erase(tk);
254  }
255  else ++tk;
256  }
257 
258  jmax.erase();
259 
260  for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
261  it != jlist.end(); ++it) jmax.addItem(*it);
262 
263  // recalculated jet quantities
264  jmax.updateJet();
265  jmax.splitted();
266  //std::cout << " jet split 1 :: " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl;
267  _members.insert(std::make_pair(jmax,jmax.pT()));
268  }
269 
270  // split and assign shared cells to nearest protojet
271  else
272  {
273  // here we need to pull the lists, because there are items to remove
274  std::list<const Item*> ilist(imax.LItems());
275  std::list<const Item*> jlist(jmax.LItems());
276 
277  typename std::list<const Item*>::iterator tk;
278  for(tk = jlist.begin(); tk != jlist.end(); )
279  {
280  typename std::list<const Item*>::iterator where=
281  find(ilist.begin(),ilist.end(),*tk);
282  if(where != ilist.end()) {
283  float yk = (*tk)->y();
284  float phik = (*tk)->phi();
285  float di= RD2(imax.y(),imax.phi(),yk,phik);
286  float dj= RD2(jmax.y(),jmax.phi(),yk,phik);
287  if(dj > di)
288  {
289  tk = jlist.erase(tk);
290  //std::cout << " shared item " << (*tk)->pT() << " removed from neighbour jet " << std::endl;
291  }
292  else
293  {
294  ilist.erase(where);
295  ++tk;
296  //std::cout << " shared item " << (*tk)->pT() << " removed from leading jet " << std::endl;
297  }
298  }
299  else ++tk;
300  }
301  // recalculate jets imax and jmax
302 
303  // first erase all items
304  imax.erase();
305  // put items that are left
306  for ( typename std::list<const Item*>::const_iterator it = ilist.begin();
307  it != ilist.end(); ++it) imax.addItem(*it);
308  // recalculated jet quantities
309  imax.updateJet();
310  imax.splitted();
311  //std::cout << " jet split 2 :: " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl;
312 
313 
314  // first erase all items
315  jmax.erase();
316  // put items that are left
317  for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
318  it != jlist.end(); ++it) jmax.addItem(*it);
319  // recalculated jet quantities
320  jmax.updateJet();
321  jmax.splitted();
322  //std::cout << " jet split " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl;
323 
324  _members.insert(std::make_pair(jmax,jmax.pT()));
325  }
326  _members.insert(std::make_pair(imax,imax.pT()));
327  }
328  } // while
329 }
330 ///////////////////////////////////////////////////////////////////////////////
331 
332 } // namespace d0
333 
334 FASTJET_END_NAMESPACE
335 
336 #endif