SISCone  2.0.6
area.cpp
1 // -*- C++ -*-
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) $//
27 
28 #include "area.h"
29 #include "momentum.h"
30 #include <stdlib.h>
31 #include <iostream>
32 
33 namespace siscone{
34 using 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 //--------------
52  active_area = passive_area = 0.0;
53 }
54 
55 // jet-initiated ctor
56 //-------------------
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 //--------------
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 //-------------
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 //-------------
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  ********************************************************************************************/
125 int 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  ********************************************************************************************/
217 int 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  ********************************************************************************************/
291 int 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 }
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Mon May 6 2013 11:30:35 for SISCone by  Doxygen 1.8.3.1