Fork me on GitHub

source: git/external/fastjet/plugins/SISCone/area.cc

Last change on this file 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.2 KB
Line 
1// -*- C++ -*-
2///////////////////////////////////////////////////////////////////////////////
3// File: area.h //
4// Description: header file for the computation of jet area //
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:: 149 $//
25// $Date:: 2007-03-15 00:13:58 +0100 (Thu, 15 Mar 2007) $//
26///////////////////////////////////////////////////////////////////////////////
27
28#include "area.h"
29#include "momentum.h"
30#include <stdlib.h>
31#include <iostream>
32
33namespace siscone{
34using namespace std;
35
36/*******************************************************
37 * Cjet_area implementation *
38 * real Jet information, including its area(s) *
39 * *
40 * This class contains information for one single jet. *
41 * That is, first, its momentum carrying information *
42 * about its centre and pT, and second, its particle *
43 * contents (from CJeT). *
44 * Compared to the Cjet class, it also includes the *
45 * passive and active areas of the jet computed using *
46 * the Carea class. *
47 *******************************************************/
48
49// default ctor
50//--------------
51Cjet_area::Cjet_area(){
52 active_area = passive_area = 0.0;
53}
54
55// jet-initiated ctor
56//-------------------
57Cjet_area::Cjet_area(Cjet &j){
58 v = j.v;
59 n = j.n;
60 contents = j.contents;
61
62 pass = j.pass;
63
64 pt_tilde = j.pt_tilde;
65 sm_var2 = j.sm_var2;
66
67 active_area = passive_area = 0.0;
68}
69
70// default dtor
71//--------------
72Cjet_area::~Cjet_area(){
73
74}
75
76
77/******************************************************************
78 * Csiscone_area implementation *
79 * class for the computation of jet areas. *
80 * *
81 * This is the class user should use whenever you want to compute *
82 * the jet area (passive and active). *
83 * It uses the SISCone algorithm to perform the jet analysis. *
84 ******************************************************************/
85
86// default ctor
87//-------------
88Carea::Carea(){
89 grid_size = 60; // 3600 particles added
90 grid_eta_max = 6.0; // maybe having twice more points in eta than in phi should be nice
91 grid_shift = 0.5; // 50% "shacking"
92
93 pt_soft = 1e-100;
94 pt_shift = 0.05;
95 pt_soft_min = 1e-90;
96}
97
98// default dtor
99//-------------
100Carea::~Carea(){
101
102}
103
104/*
105 * compute the jet areas from a given particle set.
106 * The parameters of this method are the ones which control the jet clustering alghorithm.
107 * Note that the pt_min is not allowed here soince the jet-area determination involves soft
108 * particles/jets and thus is used internally.
109 * - _particles list of particles
110 * - _radius cone radius
111 * - _f shared energy threshold for splitting&merging
112 * - _n_pass_max maximum number of passes (0=full search, the default)
113 * - _split_merge_scale the scale choice for the split-merge procedure
114 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation.
115 * SM_Et is IR safe but not boost invariant and not implemented(!)
116 * SM_mt is IR safe for hadronic events, but not for decays of two
117 * back-to-back particles of identical mass
118 * SM_pttilde
119 * is always IR safe, and also boost invariant (default)
120 * - _hard_only when this is set on, only hard jets are computed
121 * and not the purely ghosted jets (default: false)
122 * return the jets together with their areas
123 * The return value is the number of jets (including pure-ghost ones if they are included)
124 ********************************************************************************************/
125int Carea::compute_areas(std::vector<Cmomentum> &_particles, double _radius, double _f,
126 int _n_pass_max, Esplit_merge_scale _split_merge_scale,
127 bool _hard_only){
128
129 vector<Cmomentum> all_particles;
130
131 // put "hardest cut-off" if needed
132 // this avoids computation of ghosted jets when not required and
133 // significantly shortens the SM.
134 if (_hard_only){
135 SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min;
136 }
137
138 // clear potential previous runs
139 jet_areas.clear();
140
141 // put initial set of particles in the list
142 int n_hard = _particles.size();
143 all_particles = _particles;
144
145 // build the set of ghost particles and add them to the particle list
146 int i,j;
147 double eta_g,phi_g,pt_g;
148
149 for (i=0;i<grid_size;i++){
150 for (j=0;j<grid_size;j++){
151 eta_g = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
152 phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
153 pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))));
154 all_particles.push_back(Cmomentum(pt_g*cos(phi_g),pt_g*sin(phi_g),pt_g*sinh(eta_g),pt_g*cosh(eta_g)));
155 }
156 }
157
158 // run clustering with all particles.
159 // the split-merge here dynamically accounts for the purely soft jets.
160 // we therefore end up with the active area for the jets
161 int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale);
162
163 // save jets in the Cjet_area structure
164 // and determine their size
165 // jet contents is ordered by increasing index of the initial
166 // particles. Hence, we look for the first particle with index >= n_hard
167 // and deduce the number of ghosts in the jet, hence the area.
168 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size);
169
170 for (i=0;i<(int) jets.size();i++){
171 jet_areas.push_back(jets[i]);
172 j=0;
173 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
174 jet_areas[i].active_area = (jets[i].n-j)*area_factor;
175 }
176
177 // determine passive jet area
178 // for that onem we use the pt_min cut in order to remove purely
179 // soft jets from the SM procedure
180 recompute_jets(_f, pt_soft_min);
181
182 // for the area computation, we assume the jete order is the same!
183 for (i=0;i<(int) jets.size();i++){
184 j=0;
185 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
186 jet_areas[i].passive_area = (jets[i].n-j)*area_factor;
187 }
188
189 // Note:
190 // there surely remain purely soft je at the end
191 // their active area is 0 by default (see ctor)
192
193 jets.clear();
194
195 return n_jets;
196}
197
198/*
199 * compute the passive jet areas from a given particle set.
200 * The parameters of this method are the ones which control the jet clustering alghorithm.
201 * Note that the pt_min is not allowed here soince the jet-area determination involves soft
202 * particles/jets and thus is used internally.
203 * - _particles list of particles
204 * - _radius cone radius
205 * - _f shared energy threshold for splitting&merging
206 * - _n_pass_max maximum number of passes (0=full search, the default)
207 * - _split_merge_scale the scale choice for the split-merge procedure
208 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation.
209 * SM_Et is IR safe but not boost invariant and not implemented(!)
210 * SM_mt is IR safe for hadronic events, but not for decays of two
211 * back-to-back particles of identical mass
212 * SM_pttilde
213 * is always IR safe, and also boost invariant (default)
214 * return the jets together with their passive areas
215 * The return value is the number of jets
216 ********************************************************************************************/
217int Carea::compute_passive_areas(std::vector<Cmomentum> &_particles, double _radius, double _f,
218 int _n_pass_max, Esplit_merge_scale _split_merge_scale){
219
220 vector<Cmomentum> all_particles;
221
222 // in the case of passive area, we do not need
223 // to put the ghosts in the stable-cone search
224 // (they do no influence the list of stable cones)
225 // Here's how it goes...
226 stable_cone_soft_pt2_cutoff = pt_soft_min*pt_soft_min;
227
228 // clear potential previous runs
229 jet_areas.clear();
230
231 // put initial set of particles in the list
232 int n_hard = _particles.size();
233 all_particles = _particles;
234
235 // build the set of ghost particles and add them to the particle list
236 int i,j;
237 double eta_g,phi_g,pt_g;
238
239 for (i=0;i<grid_size;i++){
240 for (j=0;j<grid_size;j++){
241 eta_g = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
242 phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
243 pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))));
244 all_particles.push_back(Cmomentum(pt_g*cos(phi_g),pt_g*sin(phi_g),pt_g*sinh(eta_g),pt_g*cosh(eta_g)));
245 }
246 }
247
248 // determine passive jet area
249 // for that onem we use the pt_min cut in order to remove purely
250 // soft jets from the SM procedure
251 int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, pt_soft_min, _split_merge_scale);
252
253 // save jets in the Cjet_area structure
254 // and determine their size
255 // jet contents is ordered by increasing index of the initial
256 // particles. Hence, we look for the first particle with index >= n_hard
257 // and deduce the number of ghosts in the jet, hence the area.
258 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size);
259 for (i=0;i<(int) jets.size();i++){
260 j=0;
261 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
262 jet_areas[i].passive_area = (jets[i].n-j)*area_factor;
263 }
264
265 jets.clear();
266
267 return n_jets;
268}
269
270/*
271 * compute the active jet areas from a given particle set.
272 * The parameters of this method are the ones which control the jet clustering alghorithm.
273 * Note that the pt_min is not allowed here soince the jet-area determination involves soft
274 * particles/jets and thus is used internally.
275 * - _particles list of particles
276 * - _radius cone radius
277 * - _f shared energy threshold for splitting&merging
278 * - _n_pass_max maximum number of passes (0=full search, the default)
279 * - _split_merge_scale the scale choice for the split-merge procedure
280 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation.
281 * SM_Et is IR safe but not boost invariant and not implemented(!)
282 * SM_mt is IR safe for hadronic events, but not for decays of two
283 * back-to-back particles of identical mass
284 * SM_pttilde
285 * is always IR safe, and also boost invariant (default)
286 * - _hard_only when this is set on, only hard jets are computed
287 * and not the purely ghosted jets (default: false)
288 * return the jets together with their active areas
289 * The return value is the number of jets (including pure-ghost ones if they are included)
290 ********************************************************************************************/
291int Carea::compute_active_areas(std::vector<Cmomentum> &_particles, double _radius, double _f,
292 int _n_pass_max, Esplit_merge_scale _split_merge_scale,
293 bool _hard_only){
294
295 vector<Cmomentum> all_particles;
296
297 // put "hardest cut-off" if needed
298 // this avoids computation of ghosted jets when not required and
299 // significantly shortens the SM.
300 if (_hard_only){
301 SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min;
302 }
303
304 // clear potential previous runs
305 jet_areas.clear();
306
307 // put initial set of particles in the list
308 int n_hard = _particles.size();
309 all_particles = _particles;
310
311 // build the set of ghost particles and add them to the particle list
312 int i,j;
313 double eta_g,phi_g,pt_g;
314
315 for (i=0;i<grid_size;i++){
316 for (j=0;j<grid_size;j++){
317 eta_g = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
318 phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
319 pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))));
320 all_particles.push_back(Cmomentum(pt_g*cos(phi_g),pt_g*sin(phi_g),pt_g*sinh(eta_g),pt_g*cosh(eta_g)));
321 }
322 }
323
324 // run clustering with all particles.
325 // the split-merge here dynamically accounts for the purely soft jets.
326 // we therefore end up with the active area for the jets
327 int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale);
328
329 // save jets in the Cjet_area structure
330 // and determine their size
331 // jet contents is ordered by increasing index of the initial
332 // particles. Hence, we look for the first particle with index >= n_hard
333 // and deduce the number of ghosts in the jet, hence the area.
334 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size);
335
336 for (i=0;i<(int) jets.size();i++){
337 jet_areas.push_back(jets[i]);
338 j=0;
339 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
340 jet_areas[i].active_area = (jets[i].n-j)*area_factor;
341 }
342
343 jets.clear();
344
345 return n_jets;
346}
347
348}
Note: See TracBrowser for help on using the repository browser.