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