source: trunk/SISCone/area.cc@ 20

Last change on this file since 20 was 20, checked in by Pavel Demin, 16 years ago

add SISCone library

File size: 13.5 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: 1.1 $//
25// $Date: 2008-10-02 15:20:24 $//
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 ********************************************************************************************/
124int Carea::compute_areas(std::vector<Cmomentum> &_particles, double _radius, double _f,
125 int _n_pass_max, Esplit_merge_scale _split_merge_scale,
126 bool _hard_only){
127
128 vector<Cmomentum> all_particles;
129
130 // put "hardest cut-off" if needed
131 // this avoids computation of ghosted jets when not required and
132 // significantly shortens the SM.
133 if (_hard_only){
134 SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min;
135 }
136
137 // clear potential previous runs
138 jet_areas.clear();
139
140 // put initial set of particles in the list
141 int n_hard = _particles.size();
142 all_particles = _particles;
143
144 // build the set of ghost particles and add them to the particle list
145 int i,j;
146 double eta,phi,pt;
147
148 for (i=0;i<grid_size;i++){
149 for (j=0;j<grid_size;j++){
150 eta = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
151 phi = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
152 pt = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))));
153 all_particles.push_back(Cmomentum(pt*cos(phi),pt*sin(phi),pt*sinh(eta),pt*cosh(eta)));
154 }
155 }
156
157 // run clustering with all particles.
158 // the split-merge here dynamically accounts for the purely soft jets.
159 // we therefore end up with the active area for the jets
160 compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale);
161
162 // save jets in the Cjet_area structure
163 // and determine their size
164 // jet contents is ordered by increasing index of the initial
165 // particles. Hence, we look for the first particle with index >= n_hard
166 // and deduce the number of ghosts in the jet, hence the area.
167 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size);
168
169 for (i=0;i<(int) jets.size();i++){
170 jet_areas.push_back(jets[i]);
171 j=0;
172 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
173 jet_areas[i].active_area = (jets[i].n-j)*area_factor;
174 }
175
176 // determine passive jet area
177 // for that onem we use the pt_min cut in order to remove purely
178 // soft jets from the SM procedure
179 recompute_jets(_f, pt_soft_min);
180
181 // for the area computation, we assume the jete order is the same!
182 for (i=0;i<(int) jets.size();i++){
183 j=0;
184 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
185 jet_areas[i].passive_area = (jets[i].n-j)*area_factor;
186 }
187
188 // Note:
189 // there surely remain purely soft je at the end
190 // their active area is 0 by default (see ctor)
191
192 jets.clear();
193
194 return 0;
195}
196
197/*
198 * compute the passive jet areas from a given particle set.
199 * The parameters of this method are the ones which control the jet clustering alghorithm.
200 * Note that the pt_min is not allowed here soince the jet-area determination involves soft
201 * particles/jets and thus is used internally.
202 * - _particles list of particles
203 * - _radius cone radius
204 * - _f shared energy threshold for splitting&merging
205 * - _n_pass_max maximum number of passes (0=full search, the default)
206 * - _split_merge_scale the scale choice for the split-merge procedure
207 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation.
208 * SM_Et is IR safe but not boost invariant and not implemented(!)
209 * SM_mt is IR safe for hadronic events, but not for decays of two
210 * back-to-back particles of identical mass
211 * SM_pttilde
212 * is always IR safe, and also boost invariant (default)
213 * return the jets together with their passive areas
214 ********************************************************************************************/
215int Carea::compute_passive_areas(std::vector<Cmomentum> &_particles, double _radius, double _f,
216 int _n_pass_max, Esplit_merge_scale _split_merge_scale){
217
218 vector<Cmomentum> all_particles;
219
220 // clear potential previous runs
221 jet_areas.clear();
222
223 // put initial set of particles in the list
224 int n_hard = _particles.size();
225 all_particles = _particles;
226
227 // build the set of ghost particles and add them to the particle list
228 int i,j;
229 double eta,phi,pt;
230
231 for (i=0;i<grid_size;i++){
232 for (j=0;j<grid_size;j++){
233 eta = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
234 phi = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
235 pt = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))));
236 all_particles.push_back(Cmomentum(pt*cos(phi),pt*sin(phi),pt*sinh(eta),pt*cosh(eta)));
237 }
238 }
239
240 // determine passive jet area
241 // for that onem we use the pt_min cut in order to remove purely
242 // soft jets from the SM procedure
243 compute_jets(all_particles, _radius, _f, _n_pass_max, pt_soft_min, _split_merge_scale);
244
245 // save jets in the Cjet_area structure
246 // and determine their size
247 // jet contents is ordered by increasing index of the initial
248 // particles. Hence, we look for the first particle with index >= n_hard
249 // and deduce the number of ghosts in the jet, hence the area.
250 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size);
251 for (i=0;i<(int) jets.size();i++){
252 j=0;
253 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
254 jet_areas[i].passive_area = (jets[i].n-j)*area_factor;
255 }
256
257 jets.clear();
258
259 return 0;
260}
261
262/*
263 * compute the active jet areas from a given particle set.
264 * The parameters of this method are the ones which control the jet clustering alghorithm.
265 * Note that the pt_min is not allowed here soince the jet-area determination involves soft
266 * particles/jets and thus is used internally.
267 * - _particles list of particles
268 * - _radius cone radius
269 * - _f shared energy threshold for splitting&merging
270 * - _n_pass_max maximum number of passes (0=full search, the default)
271 * - _split_merge_scale the scale choice for the split-merge procedure
272 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation.
273 * SM_Et is IR safe but not boost invariant and not implemented(!)
274 * SM_mt is IR safe for hadronic events, but not for decays of two
275 * back-to-back particles of identical mass
276 * SM_pttilde
277 * is always IR safe, and also boost invariant (default)
278 * - _hard_only when this is set on, only hard jets are computed
279 * and not the purely ghosted jets (default: false)
280 * return the jets together with their active areas
281 ********************************************************************************************/
282int Carea::compute_active_areas(std::vector<Cmomentum> &_particles, double _radius, double _f,
283 int _n_pass_max, Esplit_merge_scale _split_merge_scale,
284 bool _hard_only){
285
286 vector<Cmomentum> all_particles;
287
288 // put "hardest cut-off" if needed
289 // this avoids computation of ghosted jets when not required and
290 // significantly shortens the SM.
291 if (_hard_only){
292 SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min;
293 }
294
295 // clear potential previous runs
296 jet_areas.clear();
297
298 // put initial set of particles in the list
299 int n_hard = _particles.size();
300 all_particles = _particles;
301
302 // build the set of ghost particles and add them to the particle list
303 int i,j;
304 double eta,phi,pt;
305
306 for (i=0;i<grid_size;i++){
307 for (j=0;j<grid_size;j++){
308 eta = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
309 phi = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size);
310 pt = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))));
311 all_particles.push_back(Cmomentum(pt*cos(phi),pt*sin(phi),pt*sinh(eta),pt*cosh(eta)));
312 }
313 }
314
315 // run clustering with all particles.
316 // the split-merge here dynamically accounts for the purely soft jets.
317 // we therefore end up with the active area for the jets
318 compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale);
319
320 // save jets in the Cjet_area structure
321 // and determine their size
322 // jet contents is ordered by increasing index of the initial
323 // particles. Hence, we look for the first particle with index >= n_hard
324 // and deduce the number of ghosts in the jet, hence the area.
325 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size);
326
327 for (i=0;i<(int) jets.size();i++){
328 jet_areas.push_back(jets[i]);
329 j=0;
330 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++;
331 jet_areas[i].active_area = (jets[i].n-j)*area_factor;
332 }
333
334 jets.clear();
335
336 return 0;
337}
338
339}
Note: See TracBrowser for help on using the repository browser.