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