Fork me on GitHub

source: svn/trunk/external/fastjet/plugins/D0RunIICone/ConeSplitMerge.hpp@ 1369

Last change on this file since 1369 was 1332, checked in by Pavel Demin, 11 years ago

upgrade fastjet to 3.0.6
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 10.6 KB
RevLine 
[859]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
67FASTJET_BEGIN_NAMESPACE
68
69namespace d0{
70
71//
72// this class is used to order ProtoJets by decreasing ET and seed ET
73template <class Item>
74class ProtoJet_ET_seedET_order
75{
76public:
[1332]77 bool operator()(const ProtoJet<Item> & first, const ProtoJet<Item> & second) const
[859]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
87template <class Item>
88class ConeSplitMerge {
89
90public :
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
100private :
101 PJMMAP _members;
102};
103///////////////////////////////////////////////////////////////////////////////
104template<class Item>
105ConeSplitMerge<Item>::ConeSplitMerge():_members() {;}
106
107template<class Item>
108ConeSplitMerge<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
121template<class Item>
122ConeSplitMerge<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
136template<class Item>
137void 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
334FASTJET_END_NAMESPACE
335
336#endif
Note: See TracBrowser for help on using the repository browser.