Changeset 273e668 in git for external/fastjet/plugins
- Timestamp:
- Oct 15, 2014, 10:55:55 AM (10 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 35b9204, b25d4cf
- Parents:
- f14bd6a
- git-author:
- Pavel Demin <pavel.demin@…> (10/10/14 08:56:40)
- git-committer:
- Pavel Demin <pavel.demin@…> (10/15/14 10:55:55)
- Location:
- external/fastjet/plugins
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
external/fastjet/plugins/D0RunICone/fastjet/D0RunIConePlugin.hh
rf14bd6a r273e668 32 32 //FJENDHEADER 33 33 34 #include "fastjet/internal/base.hh" // namespace macros (include explicitly to help Doxygen) 34 35 #include "D0RunIBaseConePlugin.hh" 35 36 -
external/fastjet/plugins/D0RunICone/fastjet/D0RunIpre96ConePlugin.hh
rf14bd6a r273e668 32 32 //FJENDHEADER 33 33 34 #include "fastjet/internal/base.hh" // namespace macros (include explicitly to help Doxygen) 34 35 #include "D0RunIBaseConePlugin.hh" 35 36 -
external/fastjet/plugins/SISCone/SISConePlugin.cc
rf14bd6a r273e668 21 21 four_vector.E); 22 22 } 23 24 //====================================================================== 25 // wrap-up around siscone's user-defined scales 26 namespace siscone_plugin_internal{ 27 /// @ingroup internal 28 /// \class SISConeUserScale 29 /// class that makes the transition between the internal SISCone 30 /// user-defined scale choice (using SISCone's Cjet) and 31 /// user-defined scale choices in the plugn above (using FastJet's 32 /// PseudoJets) 33 class SISConeUserScale : public siscone::Csplit_merge::Cuser_scale_base{ 34 public: 35 /// ctor takes the "fastjet-style" user-defined scale as well as a 36 /// reference to the current cluster sequence (to access the 37 /// particles if needed) 38 SISConeUserScale(const SISConePlugin::UserScaleBase *user_scale, 39 const ClusterSequence &cs) 40 : _user_scale(user_scale), _cs(cs){} 41 42 /// returns the scale associated to a given jet 43 virtual double operator()(const siscone::Cjet &jet) const{ 44 return _user_scale->result(_build_PJ_from_Cjet(jet)); 45 } 46 47 /// returns true id the scasle associated to jet a is larger than 48 /// the scale associated to jet b 49 virtual bool is_larger(const siscone::Cjet &a, const siscone::Cjet &b) const{ 50 return _user_scale->is_larger(_build_PJ_from_Cjet(a), _build_PJ_from_Cjet(b)); 51 } 52 53 private: 54 /// constructs a PseudoJet from a siscone::Cjet 55 /// 56 /// Note that it is tempting to overload the PseudoJet ctor. This 57 /// would not work because down the line we need to access the 58 /// original PseudoJet through the ClusterSequence and therefore 59 /// the PseudoJet structure needs to be aware of the 60 /// ClusterSequence. 61 PseudoJet _build_PJ_from_Cjet(const siscone::Cjet &jet) const{ 62 PseudoJet j(jet.v.px, jet.v.py, jet.v.pz, jet.v.E); 63 j.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>( 64 new SISConePlugin::UserScaleBaseStructureType<siscone::Cjet>(jet,_cs))); 65 return j; 66 } 67 68 const SISConePlugin::UserScaleBase *_user_scale; 69 const ClusterSequence & _cs; 70 }; 71 } 72 73 // end of the internal material 74 //====================================================================== 23 75 24 76 … … 45 97 desc << "SISCone jet algorithm with " ; 46 98 desc << "cone_radius = " << cone_radius () << ", "; 47 desc << "overlap_threshold = " << overlap_threshold () << ", "; 99 if (_progressive_removal) 100 desc << "progressive-removal mode, "; 101 else 102 desc << "overlap_threshold = " << overlap_threshold () << ", "; 48 103 desc << "n_pass_max = " << n_pass_max () << ", "; 49 104 desc << "protojet_ptmin = " << protojet_ptmin() << ", "; 50 desc << sm_scale_string << ", "; 51 desc << "caching turned " << (caching() ? on : off); 52 desc << ", SM stop scale = " << _split_merge_stopping_scale; 105 if (_progressive_removal && _user_scale) { 106 desc << "using a user-defined scale for ordering of stable cones"; 107 string user_scale_desc = _user_scale->description(); 108 if (user_scale_desc != "") desc << " (" << user_scale_desc << ")"; 109 } else { 110 desc << sm_scale_string; 111 } 112 if (!_progressive_removal){ 113 desc << ", caching turned " << (caching() ? on : off); 114 desc << ", SM stop scale = " << _split_merge_stopping_scale; 115 } 53 116 54 117 // add a note to the description if we use the pt-weighted splitting … … 85 148 bool new_siscone = true; // by default we'll be running it 86 149 87 if (caching() ) {150 if (caching() && !_progressive_removal) { 88 151 89 152 // Establish if we have a cached run with the same R, npass and … … 138 201 // run the jet finding 139 202 //cout << "plg sms: " << split_merge_scale() << endl; 140 siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(), 141 n_pass_max(), protojet_or_ghost_ptmin(), 142 Esplit_merge_scale(split_merge_scale())); 203 if (_progressive_removal){ 204 // handle the optional user-defined scale choice 205 SharedPtr<siscone_plugin_internal::SISConeUserScale> internal_scale; 206 if (_user_scale){ 207 internal_scale.reset(new siscone_plugin_internal::SISConeUserScale(_user_scale, clust_seq)); 208 siscone->set_user_scale(internal_scale.get()); 209 } 210 siscone->compute_jets_progressive_removal(siscone_momenta, cone_radius(), 211 n_pass_max(), protojet_or_ghost_ptmin(), 212 Esplit_merge_scale(split_merge_scale())); 213 } else { 214 siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(), 215 n_pass_max(), protojet_or_ghost_ptmin(), 216 Esplit_merge_scale(split_merge_scale())); 217 } 143 218 } else { 144 219 // rerun the jet finding … … 155 230 SISConeExtras * extras = new SISConeExtras(n); 156 231 232 // the ordering in which the inclusive jets are transfered here is 233 // deliberate and ensures that when a user asks for 234 // inclusive_jets(), they are provided in the order in which SISCone 235 // created them. 157 236 for (int ijet = njet-1; ijet >= 0; ijet--) { 158 237 const Cjet & jet = siscone->jets[ijet]; // shorthand … … 220 299 } 221 300 301 // //====================================================================== 302 // // material to handle user-defined scales 303 // 304 // //-------------------------------------------------- 305 // // SISCone structure type 306 // 307 // // the textual descripotion 308 // std::string SISConePlugin::UserScaleBase::StructureType::description() const{ 309 // return "PseudoJet wrapping a siscone::Cjet stable cone"; 310 // } 311 // 312 // // retrieve the constituents 313 // // 314 // // if you simply need to iterate over the constituents, it will be 315 // // faster to access them via constituent(i) 316 // vector<PseudoJet> SISConePlugin::UserScaleBase::StructureType::constituents(const PseudoJet &) const{ 317 // vector<PseudoJet> constits; 318 // constits.reserve(_jet.n); 319 // for (unsigned int i=0; i<(unsigned int)_jet.n;i++) 320 // constits.push_back(constituent(i)); 321 // return constits; 322 // } 323 // 324 // // returns the number of constituents 325 // unsigned int SISConePlugin::UserScaleBase::StructureType::size() const{ 326 // return _jet.n; 327 // } 328 // 329 // // returns the index (in the original particle list) of the ith 330 // // constituent 331 // int SISConePlugin::UserScaleBase::StructureType::constituent_index(unsigned int i) const{ 332 // return _jet.contents[i]; 333 // } 334 // 335 // // returns the ith constituent (as a PseusoJet) 336 // const PseudoJet & SISConePlugin::UserScaleBase::StructureType::constituent(unsigned int i) const{ 337 // return _cs.jets()[_jet.contents[i]]; 338 // } 339 // 340 // // returns the scalar pt of this stable cone 341 // double SISConePlugin::UserScaleBase::StructureType::pt_tilde() const{ 342 // return _jet.pt_tilde; 343 // } 344 // 345 // // returns the sm_var2 (signed ordering variable squared) for this stable cone 346 // double SISConePlugin::UserScaleBase::StructureType::ordering_var2() const{ 347 // return _jet.sm_var2; 348 // } 349 350 222 351 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh -
external/fastjet/plugins/SISCone/config.h
rf14bd6a r273e668 49 49 50 50 /* Define to the full name and version of this package. */ 51 #define PACKAGE_STRING "SISCone 2.0.6"51 #define PACKAGE_STRING "SISCone 3.0.0" 52 52 53 53 /* Define to the one symbol short name of this package. */ 54 54 #define PACKAGE_TARNAME "siscone" 55 55 56 /* Define to the home page for this package. */ 57 /* #undef PACKAGE_URL */ 58 56 59 /* Define to the version of this package. */ 57 #define PACKAGE_VERSION " 2.0.6"60 #define PACKAGE_VERSION "3.0.0" 58 61 59 62 /* Define to 1 if you have the ANSI C header files. */ … … 61 64 62 65 /* Version number of package */ 63 #define VERSION " 2.0.6"66 #define VERSION "3.0.0" -
external/fastjet/plugins/SISCone/fastjet/SISConeBasePlugin.hh
rf14bd6a r273e668 43 43 SISConeBasePlugin (){ 44 44 _use_jet_def_recombiner = false; 45 set_progressive_removal(false); 45 46 } 46 47 … … 49 50 *this = plugin; 50 51 } 52 53 /// set whether to use SISCone with progressive removal instead of 54 /// the default split_merge step. 55 /// 56 /// If progressive removal is enabled, the following SISCone 57 /// variables are not used: 58 /// 59 /// - overlap_threshold 60 /// - caching 61 /// - split_merge_stopping_scale 62 /// 63 /// The split_merge_scale choice is reinterpreted as the ordering 64 /// variable for progressive removal. It is also possible for the 65 /// user to supply his/her own function for the scale that orders 66 /// progressive removal, with set_user_scale(...) 67 void set_progressive_removal(bool progressive_removal_in=true){ 68 _progressive_removal = progressive_removal_in; 69 } 70 71 /// returns true if progressive_removal is enabled 72 bool progressive_removal() const{ return _progressive_removal;} 51 73 52 74 /// the cone radius … … 105 127 } 106 128 129 // user-defined scale for progressive removal 130 //------------------------------------------------------------ 131 132 /// \class UserScaleBase 133 /// base class for user-defined ordering of stable cones (used for 134 /// prorgessive removal) 135 /// 136 /// derived classes have to implement the () operator that returns 137 /// the scale associated with a given jet. 138 /// 139 /// It is also highly recommended to implement the is_larger() 140 /// method whenever possible, in order to avoid rounding issues 141 /// known to lead to possible infrared unsafeties. 142 /// 143 /// The jets that are passed to this class will carry the structure 144 /// of type SISConePlugin::StructureType which allows to retreive 145 /// easily the following information: 146 /// 147 /// vector<PseudoJet> constituents = jet.constituents(); 148 /// unsigned int n_constituents = jet.structure_of<SISConePlugin::UserScaleBase>().size(); 149 /// int index = jet.structure_of<SISConePlugin::UserScaleBase>().constituent_index(index i); 150 /// const PseudoJet & p = jet.structure_of<SISConePlugin::UserScaleBase>().constituent(index i); 151 /// double scalar_pt = jet.structure_of<SISConePlugin::UserScaleBase>().pt_tilde(); 152 /// 153 /// see SISConePlugin::StructureType below for further details 154 class UserScaleBase : public FunctionOfPseudoJet<double>{ 155 public: 156 /// empty virtual dtor 157 virtual ~UserScaleBase(){} 158 159 /// returns the scale associated with a given jet 160 /// 161 /// "progressive removal" iteratively removes the stable cone with 162 /// the largest scale 163 virtual double result(const PseudoJet & jet) const = 0; 164 165 /// returns true when the scale associated with jet a is larger than 166 /// the scale associated with jet b 167 /// 168 /// By default this does a simple direct comparison but it can be 169 /// overloaded for higher precision [recommended if possible] 170 virtual bool is_larger(const PseudoJet & a, const PseudoJet & b) const; 171 172 class StructureType; // defined below 173 }; 174 175 // template class derived from UserScaleBase::StryctureType that 176 // works for both SISCone jet classes 177 // implemented below 178 template<class Tjet> 179 class UserScaleBaseStructureType; 180 181 /// set a user-defined scale for stable-cone ordering in 182 /// progressive removal 183 void set_user_scale(const UserScaleBase *user_scale_in){ _user_scale = user_scale_in;} 184 185 /// returns the user-defined scale in use (0 if none) 186 const UserScaleBase * user_scale() const{ return _user_scale;} 187 188 107 189 // the things that one MUST overload required by base class 108 190 //--------------------------------------------------------- … … 120 202 double _split_merge_stopping_scale; 121 203 bool _use_jet_def_recombiner; 204 bool _progressive_removal; 122 205 123 206 mutable double _ghost_sep_scale; … … 126 209 /// call the re-clustering itself 127 210 virtual void reset_stored_plugin() const =0; 211 212 const UserScaleBase * _user_scale; 128 213 129 214 }; … … 185 270 inline SISConeBaseExtras::~SISConeBaseExtras(){} 186 271 272 //---------------------------------------------------------------------- 273 // implementation of the structure type associated with the UserScaleBase class 274 275 /// \class SISConeBasePlugin::UserScaleBase::StructureType 276 /// the structure that allows to store the information contained 277 /// into a siscone::Cjet (built internally in SISCone from a stable 278 /// cone) into a PseudoJet 279 class SISConeBasePlugin::UserScaleBase::StructureType : public PseudoJetStructureBase { 280 public: 281 /// base ctor (constructed from a ClusterSequence tin order to have 282 /// access to the initial particles 283 StructureType(const ClusterSequence &cs) 284 : _cs(cs){} 285 286 /// empty virtual dtor 287 virtual ~StructureType(){} 288 289 //-------------------------------------------------- 290 // members inherited from the base class 291 /// the textual descripotion 292 virtual std::string description() const{ 293 return "PseudoJet wrapping a siscone jet from a stable cone"; 294 } 295 296 /// this structure has constituents 297 virtual bool has_constituents() const {return true;} 298 299 /// retrieve the constituents 300 /// 301 /// if you simply need to iterate over the constituents, it will be 302 /// faster to access them via constituent(i) 303 virtual std::vector<PseudoJet> constituents(const PseudoJet & /*reference*/) const{ 304 std::vector<PseudoJet> constits; 305 constits.reserve(size()); 306 for (unsigned int i=0; i<size();i++) 307 constits.push_back(constituent(i)); 308 return constits; 309 } 310 311 //-------------------------------------------------- 312 // additional information relevant for this structure 313 314 /// returns the number of constituents 315 virtual unsigned int size() const = 0; 316 317 /// returns the index (in the original particle list) of the ith 318 /// constituent 319 virtual int constituent_index(unsigned int i) const = 0; 320 321 /// returns the ith constituent (as a PseusoJet) 322 const PseudoJet & constituent(unsigned int i) const{ 323 return _cs.jets()[constituent_index(i)]; 324 } 325 326 // /// returns the scalar pt of this stable cone 327 // virtual double pt_tilde() const = 0; 328 329 /// returns the sm_var2 (signed ordering variable squared) for this stable cone 330 virtual double ordering_var2() const = 0; 331 332 protected: 333 const ClusterSequence &_cs; ///< a reference to the CS (for access to the particles) 334 }; 335 336 337 ///@ingroup internal 338 /// template class derived from UserScaleBase::StryctureType that 339 /// works for both SISCone jet classes 340 /// implemented below 341 template<class Tjet> 342 class SISConeBasePlugin::UserScaleBaseStructureType : public UserScaleBase::StructureType{ 343 public: 344 UserScaleBaseStructureType(const Tjet &jet, const ClusterSequence &cs) 345 : UserScaleBase::StructureType(cs), _jet(jet){} 346 347 /// empty virtual dtor 348 virtual ~UserScaleBaseStructureType(){} 349 350 //-------------------------------------------------- 351 // additional information relevant for this structure 352 353 /// returns the number of constituents 354 virtual unsigned int size() const{ 355 return _jet.n; 356 } 357 358 /// returns the index (in the original particle list) of the ith 359 /// constituent 360 virtual int constituent_index(unsigned int i) const{ 361 return _jet.contents[i]; 362 } 363 364 // /// returns the scalar pt of this stable cone 365 // virtual double pt_tilde() const{ 366 // return _jet.pt_tilde; 367 // } 368 369 /// returns the sm_var2 (signed ordering variable squared) for this stable cone 370 virtual double ordering_var2() const{ 371 return _jet.sm_var2; 372 } 373 374 protected: 375 const Tjet &_jet; ///< a reference to the internal jet in SISCone 376 }; 377 187 378 188 379 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh -
external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh
rf14bd6a r273e668 7 7 namespace siscone{ 8 8 class Csiscone; 9 class Cjet; 9 10 } 10 11 … … 74 75 /// enum for the different split-merge scale choices; 75 76 /// Note that order _must_ be the same as in siscone 76 enum SplitMergeScale {SM_pt, ///< transverse momentum (E-scheme), IR unsafe 77 SM_Et, ///< transverse energy (E-scheme), not long. boost invariant 78 ///< original run-II choice [may not be implemented] 79 SM_mt, ///< transverse mass (E-scheme), IR safe except 80 ///< in decays of two identical narrow heavy particles 81 SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should 82 ///< be IR safe in all cases 77 enum SplitMergeScale { 78 SM_pt, ///< transverse momentum (E-scheme), IR unsafe 79 SM_Et, ///< transverse energy (E-scheme), not long. boost invariant 80 ///< original run-II choice [may not be implemented] 81 SM_mt, ///< transverse mass (E-scheme), IR safe except 82 ///< in decays of two identical narrow heavy particles 83 SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should 84 ///< be IR safe in all cases 83 85 }; 84 86 … … 108 110 _split_merge_stopping_scale = split_merge_stopping_scale_in; 109 111 _ghost_sep_scale = 0.0; 110 _use_pt_weighted_splitting = false;} 112 _use_pt_weighted_splitting = false; 113 _user_scale = 0;} 111 114 112 115 … … 125 128 _split_merge_stopping_scale = 0.0; 126 129 _split_merge_scale = split_merge_on_transverse_mass_in ? SM_mt : SM_pttilde; 127 _ghost_sep_scale = 0.0;} 130 _ghost_sep_scale = 0.0; 131 _user_scale = 0;} 128 132 129 133 /// backwards compatible constructor for the SISCone Plugin class … … 141 145 _split_merge_stopping_scale = 0.0; 142 146 _ghost_sep_scale = 0.0; 143 _use_pt_weighted_splitting = false;} 147 _use_pt_weighted_splitting = false; 148 _user_scale = 0;} 144 149 145 150 /// minimum pt for a protojet to be considered in the split-merge step … … 191 196 }; 192 197 198 199 /////\class SISConePlugin::UserScaleBase::StructureType 200 ///// the structure that allows to store the information contained 201 ///// into a siscone::Cjet (built internally in SISCone from a stable 202 ///// cone) into a PseudoJet 203 //class SISConePlugin::UserScaleBase::StructureType : public PseudoJetStructureBase { 204 //public: 205 // StructureType(const siscone::Cjet & jet, const ClusterSequence &cs) 206 // : _jet(jet), _cs(cs){} 207 // 208 // //-------------------------------------------------- 209 // // members inherited from the base class 210 // /// the textual descripotion 211 // virtual std::string description() const; 212 // 213 // /// this structure has constituents 214 // virtual bool has_constituents() const {return true;} 215 // 216 // /// retrieve the constituents 217 // /// 218 // /// if you simply need to iterate over the constituents, it will be 219 // /// faster to access them via constituent(i) 220 // virtual std::vector<PseudoJet> constituents(const PseudoJet & /*reference*/) const; 221 // 222 // //-------------------------------------------------- 223 // // additional information relevant for this structure 224 // 225 // /// returns the number of constituents 226 // unsigned int size() const; 227 // 228 // /// returns the index (in the original particle list) of the ith 229 // /// constituent 230 // int constituent_index(unsigned int i) const; 231 // 232 // /// returns the ith constituent (as a PseusoJet) 233 // const PseudoJet & constituent(unsigned int i) const; 234 // 235 // /// returns the scalar pt of this stable cone 236 // double pt_tilde() const; 237 // 238 // /// returns the sm_var2 (signed ordering variable squared) for this stable cone 239 // double ordering_var2() const; 240 // 241 //protected: 242 // const siscone::Cjet &_jet; ///< a dreference to the internal jet in SISCone 243 // const ClusterSequence &_cs; ///< a reference to the CS (for access to the particles) 244 //}; 193 245 194 246 //====================================================================== -
external/fastjet/plugins/SISCone/siscone.cc
rf14bd6a r273e668 21 21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 22 22 // // 23 // $Revision:: 3 41 $//24 // $Date:: 201 3-04-08 12:21:06 +0200 (Mon, 08 Apr 2013) $//23 // $Revision:: 371 $// 24 // $Date:: 2014-09-09 10:05:32 +0200 (Tue, 09 Sep 2014) $// 25 25 /////////////////////////////////////////////////////////////////////////////// 26 26 … … 29 29 //#else 30 30 //#define PACKAGE_NAME "SISCone" 31 //#define VERSION " 2.0.6"31 //#define VERSION "3.0.0" 32 32 //#warning "No config.h file available, using preset values" 33 33 //#endif … … 87 87 int _n_pass_max, double _ptmin, 88 88 Esplit_merge_scale _split_merge_scale){ 89 // initialise random number generator 90 if (!init_done){ 91 // initialise random number generator 92 ranlux_init(); 93 94 // do not do this again 95 init_done=true; 96 97 // print the banner 98 if (_banner_ostr != 0){ 99 (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl; 100 (*_banner_ostr) << "# SISCone version " << setw(28) << left << siscone_version() << "o" << endl; 101 (*_banner_ostr) << "# http://projects.hepforge.org/siscone o" << endl; 102 (*_banner_ostr) << "# o" << endl; 103 (*_banner_ostr) << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm o" << endl; 104 (*_banner_ostr) << "# SISCone was written by Gavin Salam and Gregory Soyez o" << endl; 105 (*_banner_ostr) << "# It is released under the terms of the GNU General Public License o" << endl; 106 (*_banner_ostr) << "# o" << endl; 107 (*_banner_ostr) << "# A description of the algorithm is available in the publication o" << endl; 108 (*_banner_ostr) << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)]. o" << endl; 109 (*_banner_ostr) << "# Please cite it if you use SISCone. o" << endl; 110 (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl; 111 (*_banner_ostr) << endl; 112 113 _banner_ostr->flush(); 114 } 115 } 89 _initialise_if_needed(); 116 90 117 91 // run some general safety tests (NB: f will be checked in split-merge) … … 174 148 } 175 149 150 151 /* 152 * compute the jets from a given particle set doing multiple passes 153 * such pass N looks for jets among all particles not put into jets 154 * during previous passes. 155 * - _particles list of particles 156 * - _radius cone radius 157 * - _n_pass_max maximum number of runs 158 * - _ptmin minimum pT of the protojets 159 * - _ordering_scale the ordering scale to decide which stable 160 * cone is removed 161 * return the number of jets found. 162 **********************************************************************/ 163 int Csiscone::compute_jets_progressive_removal(vector<Cmomentum> &_particles, double _radius, 164 int _n_pass_max, double _ptmin, 165 Esplit_merge_scale _ordering_scale){ 166 _initialise_if_needed(); 167 168 // run some general safety tests (NB: f will be checked in split-merge) 169 if (_radius <= 0.0 || _radius >= 0.5*M_PI) { 170 ostringstream message; 171 message << "Illegal value for cone radius, R = " << _radius 172 << " (legal values are 0<R<pi/2)"; 173 throw Csiscone_error(message.str()); 174 } 175 176 ptcomparison.split_merge_scale = _ordering_scale; 177 partial_clear(); // make sure some things are initialised properly 178 179 // init the split_merge algorithm with the initial list of particles 180 // this initialises particle list p_left of remaining particles to deal with 181 // 182 // this stores the "processed" particles in p_uncol_hard 183 init_particles(_particles); 184 jets.clear(); 185 186 bool unclustered_left; 187 rerun_allowed = false; 188 protocones_list.clear(); 189 190 do{ 191 //cout << n_left << " particle left" << endl; 192 193 // initialise stable_cone finder 194 // here we use the list of remaining particles 195 // AFTER COLLINEAR CLUSTERING !!!!!! 196 Cstable_cones::init(p_uncol_hard); 197 198 // get stable cones (stored in 'protocones') 199 unclustered_left = get_stable_cones(_radius); 200 201 // add the hardest stable cone to the list of jets 202 if (add_hardest_protocone_to_jets(&protocones, R2, _ptmin)) break; 203 204 _n_pass_max--; 205 } while ((unclustered_left) && (n_left>0) && (_n_pass_max!=0)); 206 207 // split & merge 208 return jets.size(); 209 } 210 211 176 212 /* 177 213 * recompute the jets with a different overlap parameter. … … 205 241 } 206 242 243 // ensure things are initialised 244 void Csiscone::_initialise_if_needed(){ 245 // initialise random number generator 246 if (init_done) return; 247 248 // initialise random number generator 249 ranlux_init(); 250 251 // do not do this again 252 init_done=true; 253 254 // print the banner 255 if (_banner_ostr != 0){ 256 (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl; 257 (*_banner_ostr) << "# SISCone version " << setw(28) << left << siscone_version() << "o" << endl; 258 (*_banner_ostr) << "# http://projects.hepforge.org/siscone o" << endl; 259 (*_banner_ostr) << "# o" << endl; 260 (*_banner_ostr) << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm o" << endl; 261 (*_banner_ostr) << "# SISCone was written by Gavin Salam and Gregory Soyez o" << endl; 262 (*_banner_ostr) << "# It is released under the terms of the GNU General Public License o" << endl; 263 (*_banner_ostr) << "# o" << endl; 264 (*_banner_ostr) << "# A description of the algorithm is available in the publication o" << endl; 265 (*_banner_ostr) << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)]. o" << endl; 266 (*_banner_ostr) << "# Please cite it if you use SISCone. o" << endl; 267 (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl; 268 (*_banner_ostr) << endl; 269 270 _banner_ostr->flush(); 271 } 272 } 207 273 208 274 // finally, a bunch of functions to access to -
external/fastjet/plugins/SISCone/siscone.h
rf14bd6a r273e668 22 22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 23 23 // // 24 // $Revision:: 3 25$//25 // $Date:: 201 1-11-25 12:41:17 +0100 (Fri, 25 Nov 2011) $//24 // $Revision:: 369 $// 25 // $Date:: 2014-09-04 16:57:55 +0200 (Thu, 04 Sep 2014) $// 26 26 /////////////////////////////////////////////////////////////////////////////// 27 27 … … 79 79 80 80 /** 81 * compute the jets from a given particle set. 82 * We are doing multiple passes such pass n_pass looks for jets among 83 * all particles not put into jets during previous passes. 84 * By default the number of passes is infinite (0). 85 * \param _particles list of particles 86 * \param _radius cone radius 87 * \param _n_pass_max maximum number of passes (0=full search) 88 * \param _ptmin minimum pT of the protojets 89 * \param _ordering_scale the ordering scale to decide which stable 90 * cone is removed 91 * 92 * Note that the Csplit_merge::SM_var2_hardest_cut_off cut is not 93 * used in the progressive removal variant. 94 * 95 * \return the number of jets found. 96 */ 97 int compute_jets_progressive_removal(std::vector<Cmomentum> &_particles, double _radius, 98 int _n_pass_max=0, double _ptmin=0.0, 99 Esplit_merge_scale _ordering_scale=SM_pttilde); 100 101 /** 81 102 * recompute the jets with a different overlap parameter. 82 103 * we use the same particles and R as in the preceeding call. … … 93 114 Esplit_merge_scale _split_merge_scale=SM_pttilde); 94 115 95 /// list of protocones found pass-by-pass 116 /// list of protocones found pass-by-pass (not filled by compute_jets_progressive_removal) 96 117 std::vector<std::vector<Cmomentum> > protocones_list; 97 118 … … 125 146 bool rerun_allowed; ///< is recompute_jets allowed ? 126 147 static std::ostream * _banner_ostr; ///< stream to use for banners 148 149 /// ensure things are initialised 150 void _initialise_if_needed(); 151 127 152 }; 128 153 -
external/fastjet/plugins/SISCone/split_merge.cc
rf14bd6a r273e668 21 21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 22 22 // // 23 // $Revision:: 3 22$//24 // $Date:: 201 1-11-15 10:12:36 +0100 (Tue, 15 Nov 2011) $//23 // $Revision:: 370 $// 24 // $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014) $// 25 25 /////////////////////////////////////////////////////////////////////////////// 26 26 … … 28 28 #include "siscone_error.h" 29 29 #include "momentum.h" 30 #include <math.h>31 30 #include <limits> // for max 32 31 #include <iostream> … … 34 33 #include <sstream> 35 34 #include <cassert> 35 #include <cmath> 36 36 37 37 namespace siscone{ … … 229 229 #endif 230 230 #endif 231 _user_scale = NULL; 231 232 indices = NULL; 232 233 … … 237 238 238 239 // no hardest cut (col-unsafe) 239 SM_var2_hardest_cut_off = - 1.0;240 SM_var2_hardest_cut_off = -numeric_limits<double>::max(); 240 241 241 242 // no pt cutoff for the particles to put in p_uncol_hard … … 555 556 } 556 557 558 559 /* 560 * remove the hardest protocone and declare it a a jet 561 * - protocones list of protocones (initial jet candidates) 562 * - R2 cone radius (squared) 563 * - ptmin minimal pT allowed for jets 564 * return 0 on success, 1 on error 565 * 566 * The list of remaining particles (and the uncollinear-hard ones) 567 * is updated. 568 */ 569 int Csplit_merge::add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin){ 570 571 int i; 572 Cmomentum *c; 573 Cmomentum *v; 574 double eta, phi; 575 double dx, dy; 576 double R; 577 Cjet jet, jet_candidate; 578 bool found_jet = false; 579 580 if (protocones->size()==0) 581 return 1; 582 583 pt_min2 = ptmin*ptmin; 584 R = sqrt(R2); 585 586 // browse protocones 587 // for each of them, build the list of particles in them 588 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){ 589 // initialise variables 590 c = &(*p_it); 591 592 // note: cones have been tested => their (eta,phi) coordinates are computed 593 eta = c->eta; 594 phi = c->phi; 595 596 // NOTE: this is common to this method and add_protocones, so it 597 // could be moved into a 'build_jet_from_protocone' method 598 // 599 // browse particles to create cone contents 600 jet_candidate.v = Cmomentum(); 601 jet_candidate.pt_tilde=0; 602 jet_candidate.contents.clear(); 603 for (i=0;i<n_left;i++){ 604 v = &(p_remain[i]); 605 // for security, avoid including particles with infinite rapidity) 606 // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING 607 //if (fabs(v->pz)!=v->E){ 608 dx = eta - v->eta; 609 dy = fabs(phi - v->phi); 610 if (dy>M_PI) 611 dy -= twopi; 612 if (dx*dx+dy*dy<R2){ 613 jet_candidate.contents.push_back(v->parent_index); 614 jet_candidate.v+= *v; 615 jet_candidate.pt_tilde+= pt[v->parent_index]; 616 v->index=0; 617 } 618 } 619 jet_candidate.n=jet_candidate.contents.size(); 620 621 // set the momentum in protocones 622 // (it was only known through eta and phi up to now) 623 *c = jet_candidate.v; 624 c->eta = eta; // restore exact original coords 625 c->phi = phi; // to avoid rounding error inconsistencies 626 627 // set the jet range 628 jet_candidate.range=Ceta_phi_range(eta,phi,R); 629 630 // check that the protojet has large enough pt 631 if (jet_candidate.v.perp2()<pt_min2) 632 continue; 633 634 // assign the split-merge (or progressive-removal) squared scale variable 635 if (_user_scale) { 636 // sm_var2 is the signed square of the user scale returned 637 // for the jet candidate 638 jet_candidate.sm_var2 = (*_user_scale)(jet_candidate); 639 jet_candidate.sm_var2 *= abs(jet_candidate.sm_var2); 640 } else { 641 jet_candidate.sm_var2 = get_sm_var2(jet_candidate.v, jet_candidate.pt_tilde); 642 } 643 644 // now check if it is possibly the hardest 645 if ((! found_jet) || 646 (_user_scale ? _user_scale->is_larger(jet_candidate, jet) 647 : ptcomparison(jet_candidate, jet))){ 648 jet = jet_candidate; 649 found_jet = true; 650 } 651 } 652 653 // make sure at least one of the jets has passed the selection 654 if (!found_jet) return 1; 655 656 // add the jet to the list of jets 657 jets.push_back(jet); 658 jets[jets.size()-1].v.build_etaphi(); 659 660 #ifdef DEBUG_SPLIT_MERGE 661 cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:"; 662 #endif 663 664 // update the list of what particles are left 665 int p_remain_index = 0; 666 int contents_index = 0; 667 //sort(next_jet.contents.begin(),next_jet.contents.end()); 668 for (int index=0;index<n_left;index++){ 669 if ((contents_index<(int) jet.contents.size()) && 670 (p_remain[index].parent_index == jet.contents[contents_index])){ 671 // this particle belongs to the newly found jet 672 // set pass in initial list 673 particles[p_remain[index].parent_index].index = n_pass; 674 #ifdef DEBUG_SPLIT_MERGE 675 cout << " " << jet.contents[contents_index]; 676 #endif 677 contents_index++; 678 } else { 679 // this particle still has to be clustered 680 p_remain[p_remain_index] = p_remain[index]; 681 p_remain[p_remain_index].parent_index = p_remain[index].parent_index; 682 p_remain[p_remain_index].index=1; 683 p_remain_index++; 684 } 685 } 686 p_remain.resize(n_left-jet.contents.size()); 687 n_left = p_remain.size(); 688 jets[jets.size()-1].pass = particles[jet.contents[0]].index; 689 690 // increase the pass index 691 n_pass++; 692 693 #ifdef DEBUG_SPLIT_MERGE 694 cout << endl; 695 #endif 696 697 // male sure the list of uncol_hard particles (used for the next 698 // stable cone finding) is updated [could probably be more 699 // efficient] 700 merge_collinear_and_remove_soft(); 701 702 return 0; 703 } 557 704 558 705 /* -
external/fastjet/plugins/SISCone/split_merge.h
rf14bd6a r273e668 22 22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 23 23 // // 24 // $Revision:: 268$//25 // $Date:: 20 09-03-12 21:24:16 +0100 (Thu, 12 Mar 2009) $//24 // $Revision:: 367 $// 25 // $Date:: 2014-09-04 15:57:37 +0200 (Thu, 04 Sep 2014) $// 26 26 /////////////////////////////////////////////////////////////////////////////// 27 27 … … 141 141 /// the split-merge process i.e. the variable we use for 142 142 /// 1. ordering jet candidates; 143 /// 2. computing t e overlap fraction of two candidates.143 /// 2. computing the overlap fraction of two candidates. 144 144 /// The default value uses pttile (p-scheme pt). Other alternatives are 145 145 /// pt, mt=sqrt(pt^2+m^2)=sqrt(E^2-pz^2) or Et. … … 151 151 /// the default value i.e. to use pt only for the purpose of 152 152 /// investigating the IR issue 153 /// - using Et is safe but do not respect boost invariance153 /// - using Et is safe but does not respect boost invariance 154 154 /// - using mt solves the IR unsafety issues with the pt variable 155 155 /// for QCD jets but the IR unsafety remains for nack-to-back … … 234 234 int full_clear(); 235 235 236 /////////////////////////////////////// 237 // user-defined stable-cone ordering // 238 /////////////////////////////////////// 239 240 /// \class Cuser_scale_base 241 /// base class for user-defined ordering of stable cones 242 /// 243 /// derived classes have to implement the () operator that returns 244 /// the scale associated with a given jet. 245 class Cuser_scale_base{ 246 public: 247 /// empty virtual dtor 248 virtual ~Cuser_scale_base(){} 249 250 /// the scale associated with a given jet 251 /// 252 /// "progressive removal" iteratively removes the stable cone with 253 /// the largest scale 254 virtual double operator()(const Cjet & jet) const = 0; 255 256 /// returns true when the scale associated with jet a is larger than 257 /// the scale associated with jet b 258 /// 259 /// By default this does a simple direct comparison but it can be 260 /// overloaded for higher precision [recommended if possible] 261 /// 262 /// This function assumes that a.sm_var2 and b.sm_var2 have been 263 /// correctly initialised with the signed squared output of 264 /// operator(), as is by default the case when is_larger is called 265 /// from within siscone. 266 virtual bool is_larger(const Cjet & a, const Cjet & b) const{ 267 return (a.sm_var2 > b.sm_var2); 268 } 269 }; 270 271 /// associate a user-defined scale to order the stable cones 272 /// 273 /// Note that this is only used in "progressive-removal mode", 274 /// e.g. in add_hardest_protocone_to_jets(). 275 void set_user_scale(const Cuser_scale_base * user_scale_in){ 276 _user_scale = user_scale_in; 277 } 278 279 /// return the user-defined scale (NULL if none) 280 const Cuser_scale_base * user_scale() const { return _user_scale; } 281 236 282 237 283 ///////////////////////////////// … … 256 302 */ 257 303 int add_protocones(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0); 304 305 /** 306 * remove the hardest protocone and declare it a jet 307 * \param protocones list of protocones (initial jet candidates) 308 * \param R2 cone radius (squared) 309 * \param ptmin minimal pT allowed for jets 310 * \return 0 on success, 1 on error 311 * 312 * The list of remaining particles (and the uncollinear-hard ones) 313 * is updated. 314 */ 315 int add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0); 258 316 259 317 /** … … 314 372 Csplit_merge_ptcomparison ptcomparison; 315 373 316 /// stop split--merge when the SM_var of the hardest protojet 317 /// is below this cut-off. 374 /// stop split--merge or progressive-removal when the squared SM_var 375 /// of the hardest protojet is below this cut-off. Note that this is 376 /// a signed square (ie SM_var*|SM_var|) to be able to handle 377 /// negative values. 378 /// 379 /// Note that the cut-off is set on the variable squared. 380 double SM_var2_hardest_cut_off; 381 382 /// pt cutoff for the particles to put in p_uncol_hard 383 /// this is meant to allow removing soft particles in the 384 /// stable-cone search. 385 /// 318 386 /// This is not collinear-safe so you should not use this 319 387 /// variable unless you really know what you are doing 320 388 /// Note that the cut-off is set on the variable squared. 321 double SM_var2_hardest_cut_off;322 323 /// pt cutoff for the particles to put in p_uncol_hard324 /// this is meant to allow removing soft particles in the325 /// stable-cone search.326 389 double stable_cone_soft_pt2_cutoff; 327 390 … … 390 453 bool use_pt_weighted_splitting; 391 454 455 /// use a user-defined scale to order the stable cones and jet 456 /// candidates 457 const Cuser_scale_base *_user_scale; 458 392 459 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES 393 460 /// checkxor for the candidates (to avoid having twice the same contents)
Note:
See TracChangeset
for help on using the changeset viewer.