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:
operator ()(const ProtoJet<Item> & first,const ProtoJet<Item> & second) const77 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);
~ConeSplitMerge()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>
ConeSplitMerge()105 ConeSplitMerge<Item>::ConeSplitMerge():_members() {;}
106
107 template<class Item>
ConeSplitMerge(const std::vector<ProtoJet<Item>> & jvector)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>
ConeSplitMerge(const std::list<ProtoJet<Item>> & jlist)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>
split_merge(std::vector<ProtoJet<Item>> & jcv,float shared_ET_fraction,float pT_min_leading_protojet,float pT_min_second_protojet,int MergeMax,float pT_min_noMergeMax)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
337