Fork me on GitHub

source: git/external/fastjet/plugins/SISCone/split_merge.h@ b906e2b

ImprovedOutputFile Timing dual_readout llp 3.3.3pre12
Last change on this file since b906e2b was 273e668, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

upgrade FastJet to version 3.1.0

  • Property mode set to 100644
File size: 16.7 KB
Line 
1// -*- C++ -*-
2///////////////////////////////////////////////////////////////////////////////
3// File: split_merge.h //
4// Description: header file for splitting/merging (contains the CJet class) //
5// This file is part of the SISCone project. //
6// For more details, see http://projects.hepforge.org/siscone //
7// //
8// Copyright (c) 2006 Gavin Salam and Gregory Soyez //
9// //
10// This program is free software; you can redistribute it and/or modify //
11// it under the terms of the GNU General Public License as published by //
12// the Free Software Foundation; either version 2 of the License, or //
13// (at your option) any later version. //
14// //
15// This program is distributed in the hope that it will be useful, //
16// but WITHOUT ANY WARRANTY; without even the implied warranty of //
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
18// GNU General Public License for more details. //
19// //
20// You should have received a copy of the GNU General Public License //
21// along with this program; if not, write to the Free Software //
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
23// //
24// $Revision:: 367 $//
25// $Date:: 2014-09-04 15:57:37 +0200 (Thu, 04 Sep 2014) $//
26///////////////////////////////////////////////////////////////////////////////
27
28#ifndef __SPLIT_MERGE_H__
29#define __SPLIT_MERGE_H__
30
31#include "defines.h"
32#include "geom_2d.h"
33#include "momentum.h"
34#include <stdio.h>
35#include <vector>
36#include <set>
37#include <memory>
38#include <string>
39
40namespace siscone{
41
42/**
43 * \class Cjet
44 * real Jet information.
45 *
46 * This class contains information for one single jet.
47 * That is, first, its momentum carrying information
48 * about its centre and pT, and second, its particle
49 * contents
50 */
51class Cjet{
52 public:
53 /// default ctor
54 Cjet();
55
56 /// default dtor
57 ~Cjet();
58
59 Cmomentum v; ///< jet momentum
60 double pt_tilde; ///< p-scheme pt
61 int n; ///< number of particles inside
62 std::vector<int> contents; ///< particle contents (list of indices)
63
64 /// ordering variable used for ordering and overlap in the
65 /// split--merge. This variable is automatically set either to
66 /// pt_tilde, or to mt or to pt, depending on the siscone
67 /// parameter. Note that the default behaviour is pt_tilde and that
68 /// other chices may lead to infrared unsafe situations.
69 /// Note: we use the square of the varible rather than the variable itself
70 double sm_var2;
71
72 /// covered range in eta-phi
73 Ceta_phi_range range;
74
75 /// pass at which the jet has been found
76 /// It starts at 0 (first pass), -1 means infinite rapidity
77 int pass;
78};
79
80/// ordering of jets in pt (e.g. used in final jets ordering)
81bool jets_pt_less(const Cjet &j1, const Cjet &j2);
82
83
84/// the choices of scale variable that can be used in the split-merge
85/// step, both for ordering the protojets and for measuing their
86/// overlap; pt, Et and mt=sqrt(pt^2+m^2) are all defined in E-scheme
87/// (4-momentum) recombination; pttilde = \sum_{i\in jet} |p_{t,i}|
88///
89/// NB: if one changes the order here, one _MUST_ also change the order
90/// in the SISCone plugin
91enum Esplit_merge_scale {
92 SM_pt, ///< transverse momentum (E-scheme), IR unsafe
93 SM_Et, ///< transverse energy (E-scheme), not long. boost inv.
94 ///< original run-II choice [may not be implemented]
95 SM_mt, ///< transverse mass (E-scheme), IR safe except
96 ///< in decays of two identical narrow heavy particles
97 SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should
98 ///< be IR safe in all cases
99};
100
101/// return the name of the split-merge scale choice
102std::string split_merge_scale_name(Esplit_merge_scale sms);
103
104/**
105 * \class Csplit_merge_ptcomparison
106 * comparison of jets for split--merge ordering
107 *
108 * a class that allows us to carry out comparisons of pt of jets, using
109 * information from exact particle contents where necessary.
110 */
111class Csplit_merge_ptcomparison{
112public:
113 /// default ctor
114 Csplit_merge_ptcomparison() :
115 particles(0), split_merge_scale(SM_pttilde){};
116
117 /// return the name corresponding to the SM scale variable
118 std::string SM_scale_name() const {
119 return split_merge_scale_name(split_merge_scale);}
120
121 std::vector<Cmomentum> * particles; ///< pointer to the list of particles
122 std::vector<double> * pt; ///< pointer to the pt of the particles
123
124 /// comparison between 2 jets
125 bool operator()(const Cjet &jet1, const Cjet &jet2) const;
126
127 /**
128 * get the difference between 2 jets, calculated such that rounding
129 * errors will not affect the result even if the two jets have
130 * almost the same content (so that the difference is below the
131 * rounding errors)
132 *
133 * \param j1 first jet
134 * \param j2 second jet
135 * \param v jet1-jet2
136 * \param pt_tilde jet1-jet2 pt_tilde
137 */
138 void get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const;
139
140 /// the following parameter controls the variable we're using for
141 /// the split-merge process i.e. the variable we use for
142 /// 1. ordering jet candidates;
143 /// 2. computing the overlap fraction of two candidates.
144 /// The default value uses pttile (p-scheme pt). Other alternatives are
145 /// pt, mt=sqrt(pt^2+m^2)=sqrt(E^2-pz^2) or Et.
146 /// NOTE: Modifying the default choice can have nasty effects:
147 /// - using pt leads to some IR unsafety when we have two jets,
148 /// e.g. back-to-back, with the same pt. In that case, their ordering
149 /// in pt is random and can be affected by the addition of a
150 /// soft particle. Hence, we highly recommand to keep this to
151 /// the default value i.e. to use pt only for the purpose of
152 /// investigating the IR issue
153 /// - using Et is safe but does not respect boost invariance
154 /// - using mt solves the IR unsafety issues with the pt variable
155 /// for QCD jets but the IR unsafety remains for nack-to-back
156 /// jets of unstable narrow-width particles (e.g. Higgs).
157 /// Therefore, keeping the default value is strongly advised.
158 Esplit_merge_scale split_merge_scale;
159};
160
161
162// iterator types
163/// iterator definition for the jet candidates structure
164typedef std::multiset<siscone::Cjet,Csplit_merge_ptcomparison>::iterator cjet_iterator;
165
166/// iterator definition for the jet structure
167typedef std::vector<siscone::Cjet>::iterator jet_iterator;
168
169
170
171/**
172 * \class Csplit_merge
173 * Class used to split and merge jets.
174 */
175class Csplit_merge{
176 public:
177 /// default ctor
178 Csplit_merge();
179
180 /// default dtor
181 ~Csplit_merge();
182
183
184 //////////////////////////////
185 // initialisation functions //
186 //////////////////////////////
187
188 /**
189 * initialisation function
190 * \param _particles list of particles
191 * \param protocones list of protocones (initial jet candidates)
192 * \param R2 cone radius (squared)
193 * \param ptmin minimal pT allowed for jets
194 * \return 0 on success, 1 on error
195 */
196 int init(std::vector<Cmomentum> &_particles, std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
197
198 /**
199 * initialisation function for particle list
200 * \param _particles list of particles
201 * \return 0 on success, 1 on error
202 */
203 int init_particles(std::vector<Cmomentum> &_particles);
204
205 /**
206 * build initial list of left particles
207 */
208 int init_pleft();
209
210 /**
211 * use a pt-dependent boundary for splitting
212 * When called with true, the criterium for splitting two protojets
213 * will be to compare D1^2/kt1^2 vs. D2^2/kt2^2, the (anti-)kt-weighted
214 * distance instead of the plain distance D1^2 vs. D2^2.
215 * This can be set in order to produce more circular hard jets,
216 * with the same underlying philosophy as for the anti-kt algorithm.
217 * We thus expect a behaviour closer to the IterativeCone one.
218 * By default, we use the standard D1^2 vs. D2^2 comparison and this
219 * function is not called.
220 */
221 inline int set_pt_weighted_splitting(bool _use_pt_weighted_splitting){
222 use_pt_weighted_splitting = _use_pt_weighted_splitting;
223 return 0;
224 }
225
226 ////////////////////////
227 // cleaning functions //
228 ////////////////////////
229
230 /// partial clearance
231 int partial_clear();
232
233 /// full clearance
234 int full_clear();
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
282
283 /////////////////////////////////
284 // main parts of the algorithm //
285 /////////////////////////////////
286
287 /**
288 * build the list 'p_uncol_hard' from p_remain by clustering
289 * collinear particles and removing particles softer than
290 * stable_cone_soft_pt2_cutoff
291 * note that thins in only used for stable-cone detection
292 * so the parent_index field is unnecessary
293 */
294 int merge_collinear_and_remove_soft();
295
296 /**
297 * add a list of protocones
298 * \param protocones list of protocones (initial jet candidates)
299 * \param R2 cone radius (squared)
300 * \param ptmin minimal pT allowed for jets
301 * \return 0 on success, 1 on error
302 */
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);
316
317 /**
318 * really do the splitting and merging
319 * At the end, the vector jets is filled with the jets found.
320 * the 'contents' field of each jets contains the indices
321 * of the particles included in that jet.
322 * \param overlap_tshold threshold for splitting/merging transition
323 * \param ptmin minimal pT allowed for jets
324 * \return the number of jets is returned
325 */
326 int perform(double overlap_tshold, double ptmin=0.0);
327
328
329 //////////////////////////////
330 // save and debug functions //
331 //////////////////////////////
332
333 /// save final jets
334 /// \param flux stream to save the jet contentss
335 int save_contents(FILE *flux);
336
337 /// show jets/candidates status
338 int show();
339
340 // particle information
341 int n; ///< number of particles
342 std::vector<Cmomentum> particles; ///< list of particles
343 std::vector<double> pt; ///< list of particles' pt
344 int n_left; ///< numer of particles that does not belong to any jet
345 std::vector<Cmomentum> p_remain; ///< list of particles remaining to deal with
346 std::vector<Cmomentum> p_uncol_hard; ///< list of particles remaining with collinear clustering
347 int n_pass; ///< index of the run
348
349 /// minimal difference in squared distance between a particle and
350 /// two overlapping protojets when doing a split (useful when
351 /// testing approx. collinear safety)
352 double most_ambiguous_split;
353
354 // jets information
355 std::vector<Cjet> jets; ///< list of jets
356
357 // working entries
358 int *indices; ///< maximal size array for indices works
359 int idx_size; ///< number of elements in indices1
360
361 /// The following flag indicates that identical protocones
362 /// are to be merged automatically each time around the split-merge
363 /// loop and before anything else happens.
364 ///
365 /// This flag is only effective if ALLOW_MERGE_IDENTICAL_PROTOCONES
366 /// is set in 'defines.h'
367 /// Note that this lead to infrared-unsafety so it is disabled
368 /// by default
369 bool merge_identical_protocones;
370
371 /// member used for detailed comparisons of pt's
372 Csplit_merge_ptcomparison ptcomparison;
373
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 ///
386 /// This is not collinear-safe so you should not use this
387 /// variable unless you really know what you are doing
388 /// Note that the cut-off is set on the variable squared.
389 double stable_cone_soft_pt2_cutoff;
390
391 private:
392 /**
393 * get the overlap between 2 jets
394 * \param j1 first jet
395 * \param j2 second jet
396 * \param v returned overlap^2 (determined by the choice of SM variable)
397 * \return true if overlapping, false if disjoint
398 */
399 bool get_overlap(const Cjet &j1, const Cjet &j2, double *v);
400
401
402 /**
403 * split the two given jets.
404 * during this procedure, the jets j1 & j2 are replaced
405 * by 2 new jets. Common particles are associted to the
406 * closest initial jet.
407 * \param it_j1 iterator of the first jet in 'candidates'
408 * \param it_j2 iterator of the second jet in 'candidates'
409 * \param j1 first jet (Cjet instance)
410 * \param j2 second jet (Cjet instance)
411 * \return true on success, false on error
412 */
413 bool split(cjet_iterator &it_j1, cjet_iterator &it_j2);
414
415 /**
416 * merge the two given jet.
417 * during this procedure, the jets j1 & j2 are replaced
418 * by 1 single jets containing both of them.
419 * \param it_j1 iterator of the first jet in 'candidates'
420 * \param it_j2 iterator of the second jet in 'candidates'
421 * \return true on success, false on error
422 */
423 bool merge(cjet_iterator &it_j1, cjet_iterator &it_j2);
424
425 /**
426 * Check whether or not a jet has to be inserted in the
427 * list of protojets. If it has, set its sm_variable and
428 * insert it to the list of protojets.
429 * \param jet jet to insert
430 */
431 bool insert(Cjet &jet);
432
433 /**
434 * given a 4-momentum and its associated pT, return the
435 * variable tht has to be used for SM
436 * \param v 4 momentum of the protojet
437 * \param pt_tilde pt_tilde of the protojet
438 */
439 double get_sm_var2(Cmomentum &v, double &pt_tilde);
440
441 // jet information
442 /// list of jet candidates
443 std::auto_ptr<std::multiset<Cjet,Csplit_merge_ptcomparison> > candidates;
444
445 /// minimal pt2
446 double pt_min2;
447
448 /**
449 * do we have or not to use the pt-weighted splitting
450 * (see description for set_pt_weighted_splitting)
451 * This will be false by default
452 */
453 bool use_pt_weighted_splitting;
454
455 /// use a user-defined scale to order the stable cones and jet
456 /// candidates
457 const Cuser_scale_base *_user_scale;
458
459#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
460 /// checkxor for the candidates (to avoid having twice the same contents)
461 std::set<Creference> cand_refs;
462#endif
463};
464
465}
466
467
468#endif
Note: See TracBrowser for help on using the repository browser.