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) |
---|
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 |
---|