1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
//----------------------------------------------------------------------
/// \file
/// \page Example14 14 - unified use of transformers
///
/// fastjet example program to illustrate the use of the
/// fastjet::Filter and fastjet::Pruner classes in a unified way
/// through their derivation from fastjet::Transformer
///
/// The two hardest jets in a boosted top event, clustered with an
/// (abnormally) large R, are then groomed using different tools. One
/// notes the reduction in the mass of the jets after grooming.
///
/// run it with : ./14-groomers < data/boosted_top_event.dat
///
/// Source code: 14-groomers.cc
//----------------------------------------------------------------------
//STARTHEADER
// $Id: 14-groomers.cc 2684 2011-11-14 07:41:44Z soyez $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development and are described in hep-ph/0512210. If you use
// FastJet as part of work towards a scientific publication, please
// include a citation to the FastJet paper.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
#include <fastjet/PseudoJet.hh>
#include <fastjet/ClusterSequence.hh>
#include <fastjet/Selector.hh>
#include <iostream>
#include "fastjet/tools/Filter.hh"
#include "fastjet/tools/Pruner.hh"
#include <cstdio> // needed for io
using namespace fastjet;
using namespace std;
/// an example program showing how to use Filter and Pruner in FastJet
int main(){
// read in input particles
//----------------------------------------------------------
vector<PseudoJet> input_particles;
double px, py , pz, E;
while (cin >> px >> py >> pz >> E) {
// create a PseudoJet with these components and put it onto
// back of the input_particles vector
input_particles.push_back(PseudoJet(px,py,pz,E));
}
// get the resulting jets ordered in pt
//----------------------------------------------------------
JetDefinition jet_def(cambridge_algorithm, 1.5);
ClusterSequence clust_seq(input_particles, jet_def);
vector<PseudoJet> inclusive_jets =
sorted_by_pt(clust_seq.inclusive_jets(5.0));
// label the columns
printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "mass");
// print out the details for each jet
for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
i, inclusive_jets[i].rap(), inclusive_jets[i].phi(),
inclusive_jets[i].perp(),inclusive_jets[i].m());
}
// simple test to avoid that the example below crashes:
// make sure there is at least 3 jets above our 5 GeV
if (inclusive_jets.size()<2){
cout << "Please provide an event with at least 2 jets above 5 GeV" << endl;
return 1;
}
// We will groom the two hardest jets of the event
//----------------------------------------------------------
vector<PseudoJet> candidates = SelectorNHardest(2)(inclusive_jets);
// create 3 groomers
//----------------------------------------------------------
vector<Transformer *> groomers;
// 1.
// the Cambridge/Aachen filter with Rfilt=0.3
// (simplified version of arXiv:0802.2470)
double Rfilt = 0.3;
unsigned int nfilt = 3;
groomers.push_back(new Filter(JetDefinition(cambridge_algorithm, Rfilt),
SelectorNHardest(nfilt) ) );
// 2.
// Filtering with a pt cut as for trimming (arXiv:0912.1342)
double Rtrim = 0.2;
double ptfrac = 0.03;
groomers.push_back(new Filter(JetDefinition(kt_algorithm, Rtrim),
SelectorPtFractionMin(ptfrac) ) );
// 3.
// Pruning (arXiv:0903.5081)
double zcut = 0.1;
double rcut_factor = 0.5;
groomers.push_back(new Pruner(cambridge_algorithm, zcut, rcut_factor));
// apply the various groomers to the test PseudoJet's
// and show the result
//----------------------------------------------------------
// print out original jet candidates
cout << "\nOriginal jets that will be grooomed: " << endl;
for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
const PseudoJet & c = *jit;
cout << " rap = " << c.rap() << ", phi = " << c.phi() << ", pt = " << c.perp()
<< ", mass = " << c.m()
<< " [" << c.description() << "]" << endl;
}
// loop on groomers
for (unsigned int i=0; i < groomers.size(); i++){
const Transformer & f = *groomers[i];
cout << "\nUsing groomer: " << f.description() << endl;
// loop on jet candidates
for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
const PseudoJet & c = *jit;
// apply groomer f to jet c
PseudoJet j = f(c);
// access properties specific to the given transformer
//
// We first make sure that the jet indeed has a structure
// compatible with the result of a Filter or Pruner (using
// has_structure_of()), and then retrieve the pieces rejected by the
// groomer (using structure_of())
int n_rejected;
if (j.has_structure_of<Filter>()) {
const Filter::StructureType & fj_struct = j.structure_of<Filter>();
n_rejected = fj_struct.rejected().size();
}else {
assert(j.has_structure_of<Pruner>()); // make sure
const Pruner::StructureType & fj_struct = j.structure_of<Pruner>();
n_rejected = fj_struct.rejected().size();
}
// write out result
cout << " rap = " << j.rap() << ", phi = " << j.phi() << ", pt = " << j.perp()
<< " mass = " << j.m()
<< " [kept: " << j.pieces().size() << ", rejected: "
<< n_rejected << "]" << endl;
}
}
// a bit of memory cleaning
for (unsigned int i=0; i < groomers.size(); i++) delete groomers[i];
return 0;
}
|
