FastJet 3.0.2
ConeSplitMerge.hpp
00001 #ifndef  D0RunIIconeJets_CONESPLITMERGE
00002 #define  D0RunIIconeJets_CONESPLITMERGE
00003 // ---------------------------------------------------------------------------
00004 // ConeSplitMerge.hpp
00005 //
00006 // Created: 28-JUL-2000 Francois Touze
00007 //
00008 // Purpose: Implements the pT ordered jet split/merge algorithm for the 
00009 //   Improved Legacy Cone Algorithm split/merge algo.
00010 //
00011 // Modified:
00012 //   31-JUL-2000  Laurent Duflot
00013 //     + introduce support for additional informations (ConeJetInfo)
00014 //    1-AUG-2000  Laurent Duflot
00015 //     + jet merge criteria was wrong, now calculate shared_ET and compare to 
00016 //       neighbour jet pT.
00017 //     + split was incomplete: shared items not really removed from jets.
00018 //    4-Aug-2000  Laurent Duflot
00019 //     + use item methods y() and phi() rather than p4vec() and then P2y and 
00020 //       P2phi
00021 //    7-Aug-2000  Laurent Duflot 
00022 //     + force the list to be organized by decreasing ET and, for identical ET,
00023 //       by decreasing seedET. Identical protojets can be found by eg nearby
00024 //       seeds. The seedET ordering is such that for identical jets, the one
00025 //       with the highest seedET is kept, which is what we want for efficiency
00026 //       calculations.
00027 //     + avoid unecessary copies of lists by using reference when possible
00028 //    9-Aug-2000  Laurent Duflot
00029 //     + save initial jet ET before split/merge
00030 //    1-May-2007 Lars Sonnenschein
00031 //    extracted from D0 software framework and modified to remove subsequent dependencies
00032 //
00033 //
00034 // This file is distributed with FastJet under the terms of the GNU
00035 // General Public License (v2). Permission to do so has been granted
00036 // by Lars Sonnenschein and the D0 collaboration (see COPYING for
00037 // details)
00038 //
00039 // History of changes in FastJet compared tothe original version of
00040 // ConeSplitMerge.hpp
00041 //
00042 // 2011-12-13  Gregory Soyez  <soyez@fastjet.fr>
00043 // 
00044 //        * added license information
00045 //
00046 // 2009-01-17  Gregory Soyez  <soyez@fastjet.fr>
00047 //
00048 //        * put the code in the fastjet::d0 namespace
00049 //
00050 // 2007-12-14  Gavin Salam  <salam@lpthe.jussieu.fr>
00051 // 
00052 //        * replaced make_pair by std::make_pair
00053 //
00054 // ---------------------------------------------------------------------------
00055 
00056 
00057 #include <iostream>
00058 #include <map>
00059 #include <utility>
00060 #include <vector>
00061 #include "ProtoJet.hpp"
00062 
00063 //using namespace D0RunIIconeJets_CONEJETINFO;
00064 
00065 #include <fastjet/internal/base.hh>
00066 
00067 FASTJET_BEGIN_NAMESPACE
00068 
00069 namespace d0{
00070 
00071 //
00072 // this class is used to order ProtoJets by decreasing ET and seed ET
00073 template <class Item>
00074 class ProtoJet_ET_seedET_order
00075 {
00076 public:
00077   bool operator()(const ProtoJet<Item> & first, const ProtoJet<Item> & second)
00078   {
00079     if ( first.pT() > second.pT() ) return true;
00080     else
00081       if ( first.pT() < second.pT() ) return false;
00082       else return ( first.info().seedET() > second.info().seedET() );
00083   }
00084 };
00085 
00086 
00087 template <class Item>
00088 class ConeSplitMerge {
00089 
00090 public :
00091 
00092   typedef std::multimap<ProtoJet<Item>,float,ProtoJet_ET_seedET_order<Item> > PJMMAP;
00093   
00094   ConeSplitMerge();
00095   ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector);
00096   ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist);
00097   ~ConeSplitMerge() {;}
00098   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);
00099 
00100 private :
00101   PJMMAP _members;
00102 };
00103 ///////////////////////////////////////////////////////////////////////////////
00104 template<class Item>
00105 ConeSplitMerge<Item>::ConeSplitMerge():_members() {;}
00106 
00107 template<class Item>
00108 ConeSplitMerge<Item>::ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector) 
00109 {
00110   // sort proto_jets in Et descending order
00111   typename std::vector<ProtoJet<Item> >::const_iterator jt;
00112   for(jt = jvector.begin(); jt != jvector.end(); ++jt) 
00113   {
00114     // this is supposed to be a stable cone, declare so
00115     ProtoJet<Item> jet(*jt);
00116     jet.NowStable();
00117     _members.insert(std::make_pair(jet,jet.pT()));
00118   }
00119 }
00120 
00121 template<class Item>
00122 ConeSplitMerge<Item>::ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist) 
00123 {
00124   //_max_nb_items =-1;
00125   // sort proto_jets in Et descending order
00126   typename std::list<ProtoJet<Item> >::const_iterator jt;
00127   for(jt = jlist.begin(); jt != jlist.end(); ++jt) 
00128   {
00129     // this is supposed to be a stable cone, declare so
00130     ProtoJet<Item> jet(*jt);
00131     jet.NowStable();
00132     _members.insert(std::make_pair(jet,jet.pT()));
00133   }
00134 }
00135 
00136 template<class Item>
00137 void ConeSplitMerge<Item>::split_merge(std::vector<ProtoJet<Item> >& jcv,
00138                                        float shared_ET_fraction,
00139                                        float pT_min_leading_protojet, float pT_min_second_protojet,
00140                                        int MergeMax, float pT_min_noMergeMax) 
00141 {
00142   while(!_members.empty()) 
00143   {
00144     /*
00145     {
00146       std::cout << std::endl;
00147       std::cout << " ---------------  list of protojets ------------------ " <<std::endl;
00148       for ( PJMMAP::iterator it = _members.begin();
00149             it != _members.end(); ++it)
00150       {
00151         std::cout << " pT y phi " << (*it).pT() << " " << (*it).y() << " " << (*it).phi() << " " << (*it).info().seedET() <<  " " << (*it).info().nbMerge() << " " << (*it).info().nbSplit() << std::endl;
00152       }
00153       std::cout << " ----------------------------------------------------- " <<std::endl;
00154     }
00155     */
00156 
00157 
00158     // select highest Et jet
00159     typename PJMMAP::iterator itmax= _members.begin();
00160     ProtoJet<Item> imax((*itmax).first);
00161     const std::list<const Item*>& ilist(imax.LItems());
00162 
00163     _members.erase(itmax);
00164  
00165     // does jet share items?
00166     bool share= false;
00167     float shared_ET = 0.;
00168     typename PJMMAP::iterator jtmax;
00169     typename PJMMAP::iterator jt;
00170     for(jt = _members.begin(); jt != _members.end(); ++jt) 
00171     {
00172       const std::list<const Item*>& jlist((*jt).first.LItems());
00173       typename std::list<const Item*>::const_iterator tk;
00174       int count;
00175       for(tk = ilist.begin(), count = 0; tk != ilist.end(); 
00176           ++tk, ++count) 
00177       {
00178         typename std::list<const Item*>::const_iterator where= 
00179           find(jlist.begin(),jlist.end(),*tk);   
00180         if(where != jlist.end()) 
00181         {
00182           share= true;
00183           shared_ET += (*tk)->pT();
00184         }
00185       }
00186       if(share) 
00187       {
00188         jtmax = jt;
00189         break;
00190       }
00191     }
00192     
00193     if(!share) {
00194       // add jet to the final jet list
00195       jcv.push_back(imax);
00196       //std::cout << " final jet " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
00197     }
00198     else 
00199     {
00200       // find highest Et neighbor
00201       ProtoJet<Item> jmax((*jtmax).first);
00202 
00203       // drop neighbor
00204       _members.erase(jtmax);
00205 
00206 
00207       //std::cout << " split or merge ? " << imax.pT() << " " << shared_ET << " " << jmax.pT() << " " << s << " " << (jmax.pT())*s << std::endl;
00208 
00209       // merge
00210       if( shared_ET > (jmax.pT())*shared_ET_fraction 
00211           && (imax.pT() > pT_min_leading_protojet || jmax.pT() > pT_min_second_protojet)
00212           && (imax.info().nbMerge() < MergeMax || imax.pT() > pT_min_noMergeMax))
00213         {
00214           // add neighbor's items to imax
00215           const std::list<const Item*>& jlist(jmax.LItems());
00216           typename std::list<const Item*>::const_iterator tk;
00217           typename std::list<const Item*>::const_iterator iend= ilist.end();
00218           bool same = true; // is jmax just the same as imax ? 
00219           for(tk = jlist.begin(); tk != jlist.end(); ++tk) 
00220             {
00221               typename std::list<const Item*>::const_iterator where= 
00222                 find(ilist.begin(),iend,*tk);   
00223               if(where == iend) 
00224                 {
00225                   imax.addItem(*tk);
00226                   same = false;
00227                 }
00228             }
00229           if ( ! same ) 
00230             {
00231               // recalculate
00232               //float old_pT = imax.pT();
00233               
00234               imax.updateJet();
00235               imax.merged();
00236               //std::cout << " jet merge :: " << old_pT << " " << jmax.pT() << " " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
00237             }
00238         }
00239       
00240       //split and assign removed shared cells from lowest pT protojet
00241       else if(shared_ET > (jmax.pT())*shared_ET_fraction)
00242       {
00243         // here we need to pull the lists, because there are items to remove                                                                           
00244         std::list<const Item*> ilist(imax.LItems());
00245         std::list<const Item*> jlist(jmax.LItems());
00246 
00247         typename std::list<const Item*>::iterator tk;
00248         for(tk = jlist.begin(); tk != jlist.end(); )
00249           {
00250             typename std::list<const Item*>::iterator where=
00251               find(ilist.begin(),ilist.end(),*tk);
00252             if(where != ilist.end()) {
00253               tk = jlist.erase(tk);
00254             }
00255             else ++tk;
00256           }
00257         
00258         jmax.erase();
00259 
00260         for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
00261               it != jlist.end(); ++it) jmax.addItem(*it);
00262 
00263         // recalculated jet quantities 
00264         jmax.updateJet();
00265         jmax.splitted();
00266         //std::cout << " jet split 1 :: " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl;                         
00267         _members.insert(std::make_pair(jmax,jmax.pT()));
00268       }
00269 
00270       // split and assign shared cells to nearest protojet
00271       else 
00272       {
00273         // here we need to pull the lists, because there are items to remove
00274         std::list<const Item*> ilist(imax.LItems());
00275         std::list<const Item*> jlist(jmax.LItems());
00276         
00277         typename std::list<const Item*>::iterator tk;
00278         for(tk = jlist.begin(); tk != jlist.end(); ) 
00279         {
00280           typename std::list<const Item*>::iterator where= 
00281             find(ilist.begin(),ilist.end(),*tk);
00282           if(where != ilist.end()) {
00283             float yk   = (*tk)->y();
00284             float phik = (*tk)->phi();
00285             float di= RD2(imax.y(),imax.phi(),yk,phik);
00286             float dj= RD2(jmax.y(),jmax.phi(),yk,phik);
00287             if(dj > di) 
00288             {
00289               tk = jlist.erase(tk);
00290               //std::cout << " shared item " << (*tk)->pT() << " removed from neighbour jet " << std::endl;
00291             }
00292             else
00293             {
00294               ilist.erase(where);
00295               ++tk;
00296               //std::cout << " shared item " << (*tk)->pT() << " removed from leading jet " << std::endl;
00297             }
00298           }
00299           else ++tk;
00300         }
00301         // recalculate jets imax and jmax
00302         
00303         // first erase all items
00304         imax.erase();
00305         // put items that are left
00306         for ( typename std::list<const Item*>::const_iterator it = ilist.begin();
00307               it != ilist.end(); ++it) imax.addItem(*it);
00308         // recalculated jet quantities
00309         imax.updateJet();
00310         imax.splitted();
00311         //std::cout << " jet split 2 :: " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
00312 
00313 
00314         // first erase all items
00315         jmax.erase();
00316         // put items that are left
00317         for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
00318               it != jlist.end(); ++it) jmax.addItem(*it);
00319         // recalculated jet quantities
00320         jmax.updateJet();
00321         jmax.splitted();
00322         //std::cout << " jet split " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl; 
00323 
00324         _members.insert(std::make_pair(jmax,jmax.pT()));
00325       }
00326       _members.insert(std::make_pair(imax,imax.pT()));
00327     }
00328   } // while
00329 }
00330 ///////////////////////////////////////////////////////////////////////////////
00331 
00332 }  // namespace d0
00333 
00334 FASTJET_END_NAMESPACE
00335 
00336 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends