EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-02-22 10:55:45
0001 // -*- C++ -*- 0002 /////////////////////////////////////////////////////////////////////////////// 0003 // File: split_merge.h // 0004 // Description: header file for splitting/merging (contains the CJet class) // 0005 // This file is part of the SISCone project. // 0006 // WARNING: this is not the main SISCone trunk but // 0007 // an adaptation to spherical coordinates // 0008 // For more details, see http://projects.hepforge.org/siscone // 0009 // // 0010 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez // 0011 // // 0012 // This program is free software; you can redistribute it and/or modify // 0013 // it under the terms of the GNU General Public License as published by // 0014 // the Free Software Foundation; either version 2 of the License, or // 0015 // (at your option) any later version. // 0016 // // 0017 // This program is distributed in the hope that it will be useful, // 0018 // but WITHOUT ANY WARRANTY; without even the implied warranty of // 0019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // 0020 // GNU General Public License for more details. // 0021 // // 0022 // You should have received a copy of the GNU General Public License // 0023 // along with this program; if not, write to the Free Software // 0024 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 0025 // // 0026 // $Revision:: $// 0027 // $Date:: $// 0028 /////////////////////////////////////////////////////////////////////////////// 0029 0030 #ifndef __SPH_SPLIT_MERGE_H__ 0031 #define __SPH_SPLIT_MERGE_H__ 0032 0033 #include <siscone/defines.h> 0034 #include <siscone/config.h> 0035 #include "geom_2d.h" 0036 #include "momentum.h" 0037 #include <stdio.h> 0038 #include <vector> 0039 #include <set> 0040 #include <memory> 0041 #include <string> 0042 0043 namespace siscone_spherical{ 0044 0045 const int CJET_INEXISTENT_PASS = -2; 0046 0047 /** 0048 * \class CSphjet 0049 * real Jet information. 0050 * 0051 * This class contains information for one single jet. 0052 * That is, first, its momentum carrying information 0053 * about its centre and pT, and second, its particle 0054 * contents 0055 */ 0056 class CSphjet{ 0057 public: 0058 /// default ctor 0059 CSphjet(); 0060 0061 /// default dtor 0062 ~CSphjet(); 0063 0064 CSphmomentum v; ///< jet momentum 0065 double E_tilde; ///< sum of E_i [ 1 +|p_i x p_J|^2/(|p_i|^2 E_J^2)] 0066 int n; ///< number of particles inside 0067 std::vector<int> contents; ///< particle contents (list of indices) 0068 0069 /// ordering variable used for ordering and overlap in the 0070 /// split--merge. This variable is automatically set either to 0071 /// E_tilde or E depending on the siscone parameter. 0072 /// 0073 /// Note that the default behaviour is E_tilde and that other 0074 /// choices may lead to infrared unsafe situations. 0075 /// 0076 /// Note: we use the square of the varible rather than the variable 0077 /// itself 0078 /// 0079 /// Note 2: for the overlap computation, we shall use the jet energy! 0080 double sm_var2; 0081 0082 /// covered range in eta-phi 0083 CSphtheta_phi_range range; 0084 0085 /// pass at which the jet has been found 0086 /// It starts at 0 (first pass), -1 means infinite rapidity 0087 /// (it will be initialised to "CJET_INEXISTENT_PASS" which should 0088 /// never appear after clustering) 0089 int pass; 0090 }; 0091 0092 /// ordering of jets in pt 0093 ///bool jets_pt_less(const CSphjet &j1, const CSphjet &j2); 0094 0095 /// ordering of jets in energy (e.g. used in final jets ordering) 0096 bool jets_E_less(const CSphjet &j1, const CSphjet &j2); 0097 0098 0099 /// the choices of scale variable that can be used in the split-merge 0100 /// step, both for ordering the protojets and for measuing their 0101 /// overlap; E is defined in E-scheme (4-momentum) recombination; 0102 /// Etilde = \sum_{i\in jet} E_i [1+sin^2(theta_{i,jet})] 0103 /// 0104 /// NB: if one changes the order here, one _MUST_ also change the order 0105 /// in the SISCone plugin 0106 enum Esplit_merge_scale { 0107 SM_E, ///< Energy (IR unsafe with momentum conservation) 0108 SM_Etilde ///< sum_{i \in jet} E_i [1+sin^2(theta_iJ)] 0109 }; 0110 0111 /// return the name of the split-merge scale choice 0112 std::string split_merge_scale_name(Esplit_merge_scale sms); 0113 0114 /** 0115 * \class CSphsplit_merge_ptcomparison 0116 * a class that allows us to carry out comparisons of pt of jets, using 0117 * information from exact particle contents where necessary. 0118 */ 0119 class CSphsplit_merge_ptcomparison{ 0120 public: 0121 /// default ctor 0122 CSphsplit_merge_ptcomparison() : 0123 particles(0), split_merge_scale(SM_Etilde){}; 0124 0125 /// return the name corresponding to the SM scale variable 0126 std::string SM_scale_name() const { 0127 return split_merge_scale_name(split_merge_scale);} 0128 0129 std::vector<CSphmomentum> * particles; ///< pointer to the list of particles 0130 std::vector<double> * particles_norm2; ///< pointer to the particles's norm^2 0131 0132 /// comparison of 2 CSphjet 0133 bool operator()(const CSphjet &jet1, const CSphjet &jet2) const; 0134 0135 /** 0136 * get the difference between 2 jets, calculated such that rounding 0137 * errors will not affect the result even if the two jets have 0138 * almost the same content (so that the difference is below the 0139 * rounding errors) 0140 * 0141 * \param j1 first jet 0142 * \param j2 second jet 0143 * \param v jet1-jet2 0144 * \param E_tilde jet1-jet2 E_tilde 0145 */ 0146 void get_difference(const CSphjet &j1, const CSphjet &j2, CSphmomentum *v, double *E_tilde) const; 0147 0148 /// the following parameter controls the variable we're using for 0149 /// the split-merge process i.e. the variable we use for 0150 /// 1. ordering jet candidates; 0151 /// 2. computing the overlap fraction of two candidates. 0152 /// The default value uses Etilde. The other alternative is E 0153 /// NOTE: Modifying the default choice can have nasty effects: 0154 /// - using E is IR-safe for QCD jets but the IR unsafety remains 0155 /// for back-to-back jets of unstable narrow-width particles 0156 /// (e.g. Higgs). 0157 /// Therefore, keeping the default value is strongly advised. 0158 Esplit_merge_scale split_merge_scale; 0159 }; 0160 0161 0162 // iterator types 0163 /// iterator definition for the jet candidates structure 0164 typedef std::multiset<siscone_spherical::CSphjet,CSphsplit_merge_ptcomparison>::iterator cjet_iterator; 0165 0166 /// iterator definition for the jet structure 0167 typedef std::vector<siscone_spherical::CSphjet>::iterator jet_iterator; 0168 0169 0170 0171 /** 0172 * \class CSphsplit_merge 0173 * Class used to split and merge jets. 0174 */ 0175 class CSphsplit_merge{ 0176 public: 0177 /// default ctor 0178 CSphsplit_merge(); 0179 0180 /// default dtor 0181 ~CSphsplit_merge(); 0182 0183 0184 ////////////////////////////// 0185 // initialisation functions // 0186 ////////////////////////////// 0187 0188 /** 0189 * initialisation function 0190 * \param _particles list of particles 0191 * \param protocones list of protocones (initial jet candidates) 0192 * \param R2 cone radius (squared) 0193 * \param Emin minimal energy allowed for jets 0194 * \return 0 on success, 1 on error 0195 */ 0196 int init(std::vector<CSphmomentum> &_particles, std::vector<CSphmomentum> *protocones, double R2, double Emin=0.0); 0197 0198 /** 0199 * initialisation function for particle list 0200 * \param _particles list of particles 0201 * \return 0 on success, 1 on error 0202 */ 0203 int init_particles(std::vector<CSphmomentum> &_particles); 0204 0205 /** 0206 * build initial list of left particles 0207 */ 0208 int init_pleft(); 0209 0210 /** 0211 * use an energy-dependent boundary for splitting 0212 * When called with true, the criterium for splitting two protojets 0213 * will be to compare D1^2/kt1^2 vs. D2^2/kt2^2, the (anti-)kt-weighted 0214 * distance instead of the plain distance D1^2 vs. D2^2. 0215 * This can be set in order to produce more circular hard jets, 0216 * with the same underlying philosophy as for the anti-kt algorithm. 0217 * We thus expect a behaviour closer to the IterativeCone one. 0218 * By default, we use the standard D1^2 vs. D2^2 comparison and this 0219 * function is not called. 0220 */ 0221 inline int set_E_weighted_splitting(bool _use_E_weighted_splitting){ 0222 use_E_weighted_splitting = _use_E_weighted_splitting; 0223 return 0; 0224 } 0225 0226 //////////////////////// 0227 // cleaning functions // 0228 //////////////////////// 0229 0230 /// partial clearance 0231 int partial_clear(); 0232 0233 /// full clearance 0234 int full_clear(); 0235 0236 /////////////////////////////////////// 0237 // user-defined stable-cone ordering // 0238 /////////////////////////////////////// 0239 0240 /// \class Cuser_scale_base 0241 /// base class for user-defined ordering of stable cones 0242 /// 0243 /// derived classes have to implement the () operator that returns 0244 /// the scale associated with a given jet. 0245 class Cuser_scale_base{ 0246 public: 0247 /// empty virtual dtor 0248 virtual ~Cuser_scale_base(){} 0249 0250 /// the scale associated with a given jet 0251 /// 0252 /// "progressive removal" iteratively removes the stable cone with 0253 /// the largest scale 0254 virtual double operator()(const CSphjet & jet) const = 0; 0255 0256 /// returns true when the scale associated with jet a is larger than 0257 /// the scale associated with jet b 0258 /// 0259 /// By default this does a simple direct comparison but it can be 0260 /// overloaded for higher precision [recommended if possible] 0261 /// 0262 /// This function assumes that a.sm_var2 and b.sm_var2 have been 0263 /// correctly initialised with the signed squared output of 0264 /// operator(), as is by default the case when is_larger is called 0265 /// from within siscone. 0266 virtual bool is_larger(const CSphjet & a, const CSphjet & b) const{ 0267 return (a.sm_var2 > b.sm_var2); 0268 } 0269 }; 0270 0271 /// associate a user-defined scale to order the stable cones 0272 /// 0273 /// Note that this is only used in "progressive-removal mode", 0274 /// e.g. in add_hardest_protocone_to_jets(). 0275 void set_user_scale(const Cuser_scale_base * user_scale_in){ 0276 _user_scale = user_scale_in; 0277 } 0278 0279 /// return the user-defined scale (NULL if none) 0280 const Cuser_scale_base * user_scale() const { return _user_scale; } 0281 0282 0283 ///////////////////////////////// 0284 // main parts of the algorithm // 0285 ///////////////////////////////// 0286 0287 /** 0288 * build the list 'p_uncol_hard' from p_remain by clustering 0289 * collinear particles 0290 * note that thins in only used for stable-cone detection 0291 * so the parent_index field is unnecessary 0292 * 0293 * Note that soft particles are not removed here 0294 * This is just a remnant from the trunk version 0295 */ 0296 int merge_collinear_and_remove_soft(); 0297 0298 /** 0299 * add a list of protocones 0300 * \param protocones list of protocones (initial jet candidates) 0301 * \param R2 cone radius (squared) 0302 * \param Emin minimal emergy allowed for jets 0303 * \return 0 on success, 1 on error 0304 */ 0305 int add_protocones(std::vector<CSphmomentum> *protocones, double R2, double Emin=0.0); 0306 0307 /** 0308 * remove the hardest protocone and declare it a a jet 0309 * \param protocones list of protocones (initial jet candidates) 0310 * \param R2 cone radius (squared) 0311 * \param Emin minimal emergy allowed for jets 0312 * \return 0 on success, 1 on error 0313 * 0314 * The list of remaining particles (and the uncollinear-hard ones) 0315 * is updated. 0316 */ 0317 int add_hardest_protocone_to_jets(std::vector<CSphmomentum> *protocones, double R2, double Emin=0.0); 0318 0319 /** 0320 * really do the splitting and merging 0321 * At the end, the vector jets is filled with the jets found. 0322 * the 'contents' field of each jets contains the indices 0323 * of the particles included in that jet. 0324 * \param overlap_tshold threshold for splitting/merging transition 0325 * \param Emin minimal energy allowed for jets 0326 * \return the number of jets is returned 0327 */ 0328 int perform(double overlap_tshold, double Emin=0.0); 0329 0330 0331 ////////////////////////////// 0332 // save and debug functions // 0333 ////////////////////////////// 0334 0335 /// save final jets 0336 /// \param flux stream to save the jet contentss 0337 int save_contents(FILE *flux); 0338 0339 /// show jets/candidates status 0340 int show(); 0341 0342 // particle information 0343 int n; ///< number of particles 0344 std::vector<CSphmomentum> particles; ///< list of particles 0345 std::vector<double> particles_norm2; ///< norm^2 of the particle (3-vect part) 0346 int n_left; ///< numer of particles that does not belong to any jet 0347 std::vector<CSphmomentum> p_remain; ///< list of particles remaining to deal with 0348 std::vector<CSphmomentum> p_uncol_hard; ///< list of particles remaining with collinear clustering 0349 int n_pass; ///< index of the run 0350 0351 /// minimal difference in squared distance between a particle and 0352 /// two overlapping protojets when doing a split (useful when 0353 /// testing approx. collinear safety) 0354 double most_ambiguous_split; 0355 0356 // jets information 0357 std::vector<CSphjet> jets; ///< list of jets 0358 0359 // working entries 0360 int *indices; ///< maximal size array for indices works 0361 int idx_size; ///< number of elements in indices1 0362 0363 /// The following flag indicates that identical protocones 0364 /// are to be merged automatically each time around the split-merge 0365 /// loop and before anything else happens. 0366 /// 0367 /// This flag is only effective if ALLOW_MERGE_IDENTICAL_PROTOCONES 0368 /// is set in 'defines.h' 0369 /// Note that this lead to infrared-unsafety so it is disabled 0370 /// by default 0371 bool merge_identical_protocones; 0372 0373 /// member used for detailed comparisons of pt's 0374 CSphsplit_merge_ptcomparison ptcomparison; 0375 0376 /// stop split--merge or progressive-removal when the squared SM_var 0377 /// of the hardest protojet is below this cut-off. Note that this is 0378 /// a signed square (ie SM_var*|SM_var|) to be able to handle 0379 /// negative values. 0380 /// 0381 /// Note that the cut-off is set on the variable squared. 0382 double SM_var2_hardest_cut_off; 0383 0384 /// Energy cutoff for the particles to put in p_uncol_hard 0385 /// this is meant to allow removing soft particles in the 0386 /// stable-cone search. 0387 /// 0388 /// This is not collinear-safe so you should not use this 0389 /// variable unless you really know what you are doing 0390 /// Note that the cut-off is set on the variable squared. 0391 double stable_cone_soft_E2_cutoff; 0392 0393 private: 0394 /** 0395 * get the overlap between 2 jets 0396 * \param j1 first jet 0397 * \param j2 second jet 0398 * \param v returned overlap^2 (determined by the choice of SM variable) 0399 * \return true if overlapping, false if disjoint 0400 */ 0401 bool get_overlap(const CSphjet &j1, const CSphjet &j2, double *v); 0402 0403 0404 /** 0405 * split the two given jets. 0406 * during this procedure, the jets j1 & j2 are replaced 0407 * by 2 new jets. Common particles are associted to the 0408 * closest initial jet. 0409 * \param it_j1 iterator of the first jet in 'candidates' 0410 * \param it_j2 iterator of the second jet in 'candidates' 0411 * \param j1 first jet (CSphjet instance) 0412 * \param j2 second jet (CSphjet instance) 0413 * \return true on success, false on error 0414 */ 0415 bool split(cjet_iterator &it_j1, cjet_iterator &it_j2); 0416 0417 /** 0418 * merge the two given jet. 0419 * during this procedure, the jets j1 & j2 are replaced 0420 * by 1 single jets containing both of them. 0421 * \param it_j1 iterator of the first jet in 'candidates' 0422 * \param it_j2 iterator of the second jet in 'candidates' 0423 * \return true on success, false on error 0424 */ 0425 bool merge(cjet_iterator &it_j1, cjet_iterator &it_j2); 0426 0427 /** 0428 * Check whether or not a jet has to be inserted in the 0429 * list of protojets. If it has, set its sm_variable and 0430 * insert it to the list of protojets. 0431 * \param jet jet to insert 0432 */ 0433 bool insert(CSphjet &jet); 0434 0435 /** 0436 * given a 4-momentum and its associated pT, return the 0437 * variable tht has to be used for SM 0438 * \param v 4 momentum of the protojet 0439 * \param E_tilde E_tilde of the protojet 0440 */ 0441 double get_sm_var2(CSphmomentum &v, double &E_tilde); 0442 0443 /// compute Etilde for a given jet 0444 void compute_Etilde(CSphjet &j); 0445 0446 // jet information 0447 /// list of jet candidates 0448 #ifdef SISCONE_USES_UNIQUE_PTR_AS_AUTO_PTR 0449 std::unique_ptr<std::multiset<CSphjet,CSphsplit_merge_ptcomparison> > candidates; 0450 #else 0451 std::auto_ptr<std::multiset<CSphjet,CSphsplit_merge_ptcomparison> > candidates; 0452 #endif 0453 0454 /// minimal E 0455 double E_min; 0456 0457 /** 0458 * do we have or not to use the energy-weighted splitting 0459 * (see description for set_E_weighted_splitting) 0460 * This will be false by default 0461 */ 0462 bool use_E_weighted_splitting; 0463 0464 /// use a user-defined scale to order the stable cones and jet 0465 /// candidates 0466 const Cuser_scale_base *_user_scale; 0467 0468 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES 0469 /// checkxor for the candidates (to avoid having twice the same contents) 0470 std::set<siscone::Creference> cand_refs; 0471 #endif 0472 }; 0473 0474 } 0475 0476 0477 #endif
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|