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