FastJet 3.0.0
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 #include <iostream>
00035 #include <map>
00036 #include <utility>
00037 #include <vector>
00038 #include "ProtoJet.hpp"
00039 
00040 //using namespace D0RunIIconeJets_CONEJETINFO;
00041 
00042 #include <fastjet/internal/base.hh>
00043 
00044 FASTJET_BEGIN_NAMESPACE
00045 
00046 namespace d0{
00047 
00048 //
00049 // this class is used to order ProtoJets by decreasing ET and seed ET
00050 template <class Item>
00051 class ProtoJet_ET_seedET_order
00052 {
00053 public:
00054   bool operator()(const ProtoJet<Item> & first, const ProtoJet<Item> & second)
00055   {
00056     if ( first.pT() > second.pT() ) return true;
00057     else
00058       if ( first.pT() < second.pT() ) return false;
00059       else return ( first.info().seedET() > second.info().seedET() );
00060   }
00061 };
00062 
00063 
00064 template <class Item>
00065 class ConeSplitMerge {
00066 
00067 public :
00068 
00069   typedef std::multimap<ProtoJet<Item>,float,ProtoJet_ET_seedET_order<Item> > PJMMAP;
00070   
00071   ConeSplitMerge();
00072   ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector);
00073   ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist);
00074   ~ConeSplitMerge() {;}
00075   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);
00076 
00077 private :
00078   PJMMAP _members;
00079 };
00080 ///////////////////////////////////////////////////////////////////////////////
00081 template<class Item>
00082 ConeSplitMerge<Item>::ConeSplitMerge():_members() {;}
00083 
00084 template<class Item>
00085 ConeSplitMerge<Item>::ConeSplitMerge(const std::vector<ProtoJet<Item> >& jvector) 
00086 {
00087   // sort proto_jets in Et descending order
00088   typename std::vector<ProtoJet<Item> >::const_iterator jt;
00089   for(jt = jvector.begin(); jt != jvector.end(); ++jt) 
00090   {
00091     // this is supposed to be a stable cone, declare so
00092     ProtoJet<Item> jet(*jt);
00093     jet.NowStable();
00094     _members.insert(std::make_pair(jet,jet.pT()));
00095   }
00096 }
00097 
00098 template<class Item>
00099 ConeSplitMerge<Item>::ConeSplitMerge(const std::list<ProtoJet<Item> >& jlist) 
00100 {
00101   //_max_nb_items =-1;
00102   // sort proto_jets in Et descending order
00103   typename std::list<ProtoJet<Item> >::const_iterator jt;
00104   for(jt = jlist.begin(); jt != jlist.end(); ++jt) 
00105   {
00106     // this is supposed to be a stable cone, declare so
00107     ProtoJet<Item> jet(*jt);
00108     jet.NowStable();
00109     _members.insert(std::make_pair(jet,jet.pT()));
00110   }
00111 }
00112 
00113 template<class Item>
00114 void ConeSplitMerge<Item>::split_merge(std::vector<ProtoJet<Item> >& jcv,
00115                                        float shared_ET_fraction,
00116                                        float pT_min_leading_protojet, float pT_min_second_protojet,
00117                                        int MergeMax, float pT_min_noMergeMax) 
00118 {
00119   while(!_members.empty()) 
00120   {
00121     /*
00122     {
00123       std::cout << std::endl;
00124       std::cout << " ---------------  list of protojets ------------------ " <<std::endl;
00125       for ( PJMMAP::iterator it = _members.begin();
00126             it != _members.end(); ++it)
00127       {
00128         std::cout << " pT y phi " << (*it).pT() << " " << (*it).y() << " " << (*it).phi() << " " << (*it).info().seedET() <<  " " << (*it).info().nbMerge() << " " << (*it).info().nbSplit() << std::endl;
00129       }
00130       std::cout << " ----------------------------------------------------- " <<std::endl;
00131     }
00132     */
00133 
00134 
00135     // select highest Et jet
00136     typename PJMMAP::iterator itmax= _members.begin();
00137     ProtoJet<Item> imax((*itmax).first);
00138     const std::list<const Item*>& ilist(imax.LItems());
00139 
00140     _members.erase(itmax);
00141  
00142     // does jet share items?
00143     bool share= false;
00144     float shared_ET = 0.;
00145     typename PJMMAP::iterator jtmax;
00146     typename PJMMAP::iterator jt;
00147     for(jt = _members.begin(); jt != _members.end(); ++jt) 
00148     {
00149       const std::list<const Item*>& jlist((*jt).first.LItems());
00150       typename std::list<const Item*>::const_iterator tk;
00151       int count;
00152       for(tk = ilist.begin(), count = 0; tk != ilist.end(); 
00153           ++tk, ++count) 
00154       {
00155         typename std::list<const Item*>::const_iterator where= 
00156           find(jlist.begin(),jlist.end(),*tk);   
00157         if(where != jlist.end()) 
00158         {
00159           share= true;
00160           shared_ET += (*tk)->pT();
00161         }
00162       }
00163       if(share) 
00164       {
00165         jtmax = jt;
00166         break;
00167       }
00168     }
00169     
00170     if(!share) {
00171       // add jet to the final jet list
00172       jcv.push_back(imax);
00173       //std::cout << " final jet " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
00174     }
00175     else 
00176     {
00177       // find highest Et neighbor
00178       ProtoJet<Item> jmax((*jtmax).first);
00179 
00180       // drop neighbor
00181       _members.erase(jtmax);
00182 
00183 
00184       //std::cout << " split or merge ? " << imax.pT() << " " << shared_ET << " " << jmax.pT() << " " << s << " " << (jmax.pT())*s << std::endl;
00185 
00186       // merge
00187       if( shared_ET > (jmax.pT())*shared_ET_fraction 
00188           && (imax.pT() > pT_min_leading_protojet || jmax.pT() > pT_min_second_protojet)
00189           && (imax.info().nbMerge() < MergeMax || imax.pT() > pT_min_noMergeMax))
00190         {
00191           // add neighbor's items to imax
00192           const std::list<const Item*>& jlist(jmax.LItems());
00193           typename std::list<const Item*>::const_iterator tk;
00194           typename std::list<const Item*>::const_iterator iend= ilist.end();
00195           bool same = true; // is jmax just the same as imax ? 
00196           for(tk = jlist.begin(); tk != jlist.end(); ++tk) 
00197             {
00198               typename std::list<const Item*>::const_iterator where= 
00199                 find(ilist.begin(),iend,*tk);   
00200               if(where == iend) 
00201                 {
00202                   imax.addItem(*tk);
00203                   same = false;
00204                 }
00205             }
00206           if ( ! same ) 
00207             {
00208               // recalculate
00209               //float old_pT = imax.pT();
00210               
00211               imax.updateJet();
00212               imax.merged();
00213               //std::cout << " jet merge :: " << old_pT << " " << jmax.pT() << " " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
00214             }
00215         }
00216       
00217       //split and assign removed shared cells from lowest pT protojet
00218       else if(shared_ET > (jmax.pT())*shared_ET_fraction)
00219       {
00220         // here we need to pull the lists, because there are items to remove                                                                           
00221         std::list<const Item*> ilist(imax.LItems());
00222         std::list<const Item*> jlist(jmax.LItems());
00223 
00224         typename std::list<const Item*>::iterator tk;
00225         for(tk = jlist.begin(); tk != jlist.end(); )
00226           {
00227             typename std::list<const Item*>::iterator where=
00228               find(ilist.begin(),ilist.end(),*tk);
00229             if(where != ilist.end()) {
00230               tk = jlist.erase(tk);
00231             }
00232             else ++tk;
00233           }
00234         
00235         jmax.erase();
00236 
00237         for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
00238               it != jlist.end(); ++it) jmax.addItem(*it);
00239 
00240         // recalculated jet quantities 
00241         jmax.updateJet();
00242         jmax.splitted();
00243         //std::cout << " jet split 1 :: " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl;                         
00244         _members.insert(std::make_pair(jmax,jmax.pT()));
00245       }
00246 
00247       // split and assign shared cells to nearest protojet
00248       else 
00249       {
00250         // here we need to pull the lists, because there are items to remove
00251         std::list<const Item*> ilist(imax.LItems());
00252         std::list<const Item*> jlist(jmax.LItems());
00253         
00254         typename std::list<const Item*>::iterator tk;
00255         for(tk = jlist.begin(); tk != jlist.end(); ) 
00256         {
00257           typename std::list<const Item*>::iterator where= 
00258             find(ilist.begin(),ilist.end(),*tk);
00259           if(where != ilist.end()) {
00260             float yk   = (*tk)->y();
00261             float phik = (*tk)->phi();
00262             float di= RD2(imax.y(),imax.phi(),yk,phik);
00263             float dj= RD2(jmax.y(),jmax.phi(),yk,phik);
00264             if(dj > di) 
00265             {
00266               tk = jlist.erase(tk);
00267               //std::cout << " shared item " << (*tk)->pT() << " removed from neighbour jet " << std::endl;
00268             }
00269             else
00270             {
00271               ilist.erase(where);
00272               ++tk;
00273               //std::cout << " shared item " << (*tk)->pT() << " removed from leading jet " << std::endl;
00274             }
00275           }
00276           else ++tk;
00277         }
00278         // recalculate jets imax and jmax
00279         
00280         // first erase all items
00281         imax.erase();
00282         // put items that are left
00283         for ( typename std::list<const Item*>::const_iterator it = ilist.begin();
00284               it != ilist.end(); ++it) imax.addItem(*it);
00285         // recalculated jet quantities
00286         imax.updateJet();
00287         imax.splitted();
00288         //std::cout << " jet split 2 :: " << imax.pT() << " "<< imax.info().nbMerge() << " " << imax.info().nbSplit() << std::endl; 
00289 
00290 
00291         // first erase all items
00292         jmax.erase();
00293         // put items that are left
00294         for ( typename std::list<const Item*>::const_iterator it = jlist.begin();
00295               it != jlist.end(); ++it) jmax.addItem(*it);
00296         // recalculated jet quantities
00297         jmax.updateJet();
00298         jmax.splitted();
00299         //std::cout << " jet split " << jmax.pT() << " "<< jmax.info().nbMerge() << " " << jmax.info().nbSplit() << std::endl; 
00300 
00301         _members.insert(std::make_pair(jmax,jmax.pT()));
00302       }
00303       _members.insert(std::make_pair(imax,imax.pT()));
00304     }
00305   } // while
00306 }
00307 ///////////////////////////////////////////////////////////////////////////////
00308 
00309 }  // namespace d0
00310 
00311 FASTJET_END_NAMESPACE
00312 
00313 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends