Fork me on GitHub

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

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 4b9a2dc was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 14.3 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:: 268 $//
25// $Date:: 2009-03-12 21:24:16 +0100 (Thu, 12 Mar 2009) $//
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 te 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 do 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 /////////////////////////////////
238 // main parts of the algorithm //
239 /////////////////////////////////
240
241 /**
242 * build the list 'p_uncol_hard' from p_remain by clustering
243 * collinear particles and removing particles softer than
244 * stable_cone_soft_pt2_cutoff
245 * note that thins in only used for stable-cone detection
246 * so the parent_index field is unnecessary
247 */
248 int merge_collinear_and_remove_soft();
249
250 /**
251 * add a list of protocones
252 * \param protocones list of protocones (initial jet candidates)
253 * \param R2 cone radius (squared)
254 * \param ptmin minimal pT allowed for jets
255 * \return 0 on success, 1 on error
256 */
257 int add_protocones(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
258
259 /**
260 * really do the splitting and merging
261 * At the end, the vector jets is filled with the jets found.
262 * the 'contents' field of each jets contains the indices
263 * of the particles included in that jet.
264 * \param overlap_tshold threshold for splitting/merging transition
265 * \param ptmin minimal pT allowed for jets
266 * \return the number of jets is returned
267 */
268 int perform(double overlap_tshold, double ptmin=0.0);
269
270
271 //////////////////////////////
272 // save and debug functions //
273 //////////////////////////////
274
275 /// save final jets
276 /// \param flux stream to save the jet contentss
277 int save_contents(FILE *flux);
278
279 /// show jets/candidates status
280 int show();
281
282 // particle information
283 int n; ///< number of particles
284 std::vector<Cmomentum> particles; ///< list of particles
285 std::vector<double> pt; ///< list of particles' pt
286 int n_left; ///< numer of particles that does not belong to any jet
287 std::vector<Cmomentum> p_remain; ///< list of particles remaining to deal with
288 std::vector<Cmomentum> p_uncol_hard; ///< list of particles remaining with collinear clustering
289 int n_pass; ///< index of the run
290
291 /// minimal difference in squared distance between a particle and
292 /// two overlapping protojets when doing a split (useful when
293 /// testing approx. collinear safety)
294 double most_ambiguous_split;
295
296 // jets information
297 std::vector<Cjet> jets; ///< list of jets
298
299 // working entries
300 int *indices; ///< maximal size array for indices works
301 int idx_size; ///< number of elements in indices1
302
303 /// The following flag indicates that identical protocones
304 /// are to be merged automatically each time around the split-merge
305 /// loop and before anything else happens.
306 ///
307 /// This flag is only effective if ALLOW_MERGE_IDENTICAL_PROTOCONES
308 /// is set in 'defines.h'
309 /// Note that this lead to infrared-unsafety so it is disabled
310 /// by default
311 bool merge_identical_protocones;
312
313 /// member used for detailed comparisons of pt's
314 Csplit_merge_ptcomparison ptcomparison;
315
316 /// stop split--merge when the SM_var of the hardest protojet
317 /// is below this cut-off.
318 /// This is not collinear-safe so you should not use this
319 /// variable unless you really know what you are doing
320 /// 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_hard
324 /// this is meant to allow removing soft particles in the
325 /// stable-cone search.
326 double stable_cone_soft_pt2_cutoff;
327
328 private:
329 /**
330 * get the overlap between 2 jets
331 * \param j1 first jet
332 * \param j2 second jet
333 * \param v returned overlap^2 (determined by the choice of SM variable)
334 * \return true if overlapping, false if disjoint
335 */
336 bool get_overlap(const Cjet &j1, const Cjet &j2, double *v);
337
338
339 /**
340 * split the two given jets.
341 * during this procedure, the jets j1 & j2 are replaced
342 * by 2 new jets. Common particles are associted to the
343 * closest initial jet.
344 * \param it_j1 iterator of the first jet in 'candidates'
345 * \param it_j2 iterator of the second jet in 'candidates'
346 * \param j1 first jet (Cjet instance)
347 * \param j2 second jet (Cjet instance)
348 * \return true on success, false on error
349 */
350 bool split(cjet_iterator &it_j1, cjet_iterator &it_j2);
351
352 /**
353 * merge the two given jet.
354 * during this procedure, the jets j1 & j2 are replaced
355 * by 1 single jets containing both of them.
356 * \param it_j1 iterator of the first jet in 'candidates'
357 * \param it_j2 iterator of the second jet in 'candidates'
358 * \return true on success, false on error
359 */
360 bool merge(cjet_iterator &it_j1, cjet_iterator &it_j2);
361
362 /**
363 * Check whether or not a jet has to be inserted in the
364 * list of protojets. If it has, set its sm_variable and
365 * insert it to the list of protojets.
366 * \param jet jet to insert
367 */
368 bool insert(Cjet &jet);
369
370 /**
371 * given a 4-momentum and its associated pT, return the
372 * variable tht has to be used for SM
373 * \param v 4 momentum of the protojet
374 * \param pt_tilde pt_tilde of the protojet
375 */
376 double get_sm_var2(Cmomentum &v, double &pt_tilde);
377
378 // jet information
379 /// list of jet candidates
380 std::auto_ptr<std::multiset<Cjet,Csplit_merge_ptcomparison> > candidates;
381
382 /// minimal pt2
383 double pt_min2;
384
385 /**
386 * do we have or not to use the pt-weighted splitting
387 * (see description for set_pt_weighted_splitting)
388 * This will be false by default
389 */
390 bool use_pt_weighted_splitting;
391
392#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
393 /// checkxor for the candidates (to avoid having twice the same contents)
394 std::set<Creference> cand_refs;
395#endif
396};
397
398}
399
400
401#endif
Note: See TracBrowser for help on using the repository browser.