FastJet 3.4.3
Loading...
Searching...
No Matches
09-user_info.cc
Go to the documentation of this file.
1//----------------------------------------------------------------------
2/// \file
3/// \page Example09 09 - adding user information to a fastjet::PseudoJet
4///
5/// This example illustrates how it is possible to associate
6/// user-defined information to a fastjet::PseudoJet (beyond the
7/// simple user index), using a class derived from
8/// fastjet::UserInfoBase.
9///
10/// Note that in this example we have chosen to use this
11/// user-defined information to obtain properties of the constituents
12/// of the reconstructed jet (e.g. if the event is made of a hard
13/// interaction and pileup, what part of the reconstructed jet comes
14/// from the hard interaction). To do that, we also show how to
15/// introduce a user-defined fastjet::Selector. For some applications,
16/// it might also be useful to define new recombination schemes using
17/// the extra information.
18///
19/// run it with : ./09-user_info < data/Pythia-dijet-ptmin100-lhc-pileup-1ev.dat
20///
21/// (Note that this event consists of many sub-events, the first one
22/// being the "hard" interaction and the following being minbias
23/// events composing the pileup. It has the specificity that it also
24/// contains the PDG id of the particles)
25///
26/// Source code: 09-user_info.cc
27//----------------------------------------------------------------------
28
29//FJSTARTHEADER
30// $Id$
31//
32// Copyright (c) 2005-2024, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
33//
34//----------------------------------------------------------------------
35// This file is part of FastJet.
36//
37// FastJet is free software; you can redistribute it and/or modify
38// it under the terms of the GNU General Public License as published by
39// the Free Software Foundation; either version 2 of the License, or
40// (at your option) any later version.
41//
42// The algorithms that underlie FastJet have required considerable
43// development. They are described in the original FastJet paper,
44// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
45// FastJet as part of work towards a scientific publication, please
46// quote the version you use and include a citation to the manual and
47// optionally also to hep-ph/0512210.
48//
49// FastJet is distributed in the hope that it will be useful,
50// but WITHOUT ANY WARRANTY; without even the implied warranty of
51// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
52// GNU General Public License for more details.
53//
54// You should have received a copy of the GNU General Public License
55// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
56//----------------------------------------------------------------------
57//FJENDHEADER
58
59#include "fastjet/ClusterSequence.hh"
60#include "fastjet/Selector.hh"
61#include <iostream> // needed for io
62#include <sstream> // needed for io
63#include <cstdio> // needed for io
64
65using namespace std;
66using namespace fastjet;
67
68//------------------------------------------------------------------------
69// the user information
70//
71// To associate extra information to a PseudoJet, one first has to
72// create a class, derived from UserInfoBase, that contains
73// that information.
74//
75// In our simple example, we shall use 2 informations
76// - the PDG id associated with the particle
77// - the "vertex number" associated with the particle
78class MyUserInfo : public PseudoJet::UserInfoBase{
79public:
80 // default ctor
81 // - pdg_id the PDG id of the particle
82 // - vertex_number theid of the vertex it originates from
83 MyUserInfo(const int & pdg_id_in, const int & vertex_number_in) :
84 _pdg_id(pdg_id_in), _vertex_number(vertex_number_in){}
85
86 /// access to the PDG id
87 int pdg_id() const { return _pdg_id;}
88
89 /// access to the vertex number
90 int vertex_number() const { return _vertex_number;}
91
92protected:
93 int _pdg_id; // the associated pdg id
94 int _vertex_number; // the associated vertex number
95};
96
97
98//------------------------------------------------------------------------
99// Select pi0 and photons
100//
101// This shows how we can build a Selector that uses the user-defined
102// information to select particles that are either pi0's or photons
103// (we choose this purely for simplicity).
104//
105// To create a user-defined Selector, the first step is to
106// create its associated "worker" class, i.e. to derive a class from
107// SelectorWorker. Then (see below), we just write a function
108// (SelectorIsPi0Gamma()) that creates a Selector with the
109// appropriate worker class.
110class SW_IsPi0Gamma : public SelectorWorker{
111public:
112 // default ctor
113 SW_IsPi0Gamma(){}
114
115 // the selector's description
116 string description() const{
117 return "neutral pions or photons";
118 }
119
120 // keeps the ones that have the pdg id of the pi0
121 bool pass(const PseudoJet &p) const{
122 // This is how we access the extra information associated with the
123 // particles we test.
124 // p.user_info<MyUserInfo>()
125 // returns a reference to the user-defined information (of type
126 // MyUserInfo, as mentioned explicitly). It ensures automatically
127 // that there is an associated user info compatible with the
128 // requested type (and throws an error if it is not the case)
129 //
130 // We can then access the "pdg_id" member of MyUserInfo to
131 // extract the targeted information.
132 const int & pdgid = p.user_info<MyUserInfo>().pdg_id();
133 return (pdgid == 111) || (pdgid == 22);
134 }
135};
136
137// the function that allows to write simply
138// Selector sel = SelectorIsPi0Gamma();
139Selector SelectorIsPi0Gamma(){
140 return Selector(new SW_IsPi0Gamma());
141}
142
143//------------------------------------------------------------------------
144// Select particles from a given vertex number
145//
146// This is the same kind of construct as just above except that we
147// select on particles that are originated from a given vertex. The
148// test event has been structured as a superposition of sub-events
149// (the 0th being the hard interaction) and each particle will be
150// associated a vertex number. This Selector allows to select
151// particles corresponding to a given vertex number.
152//
153// As in the previous case, we start with the worker class and then
154// write a function for the Selector itself.
155class SW_VertexNumber : public SelectorWorker{
156public:
157 // ctor from the vertex we want to keep
158 SW_VertexNumber(const int & vertex_number) : _vertex_number(vertex_number){}
159
160 // the selector's description
161 string description() const{
162 ostringstream oss;
163 oss << "vertex number " << _vertex_number;
164 return oss.str();
165 }
166
167 // keeps the ones that have the correct vertex number
168 bool pass(const PseudoJet &p) const{
169 // This is how we access the extra information associated with the
170 // particles we test.
171 // p.user_info<MyUserInfo>()
172 // returns a reference to the user-defined information (of type
173 // MyUserInfo, as mentioned explicitly). It ensures automatically
174 // that there is an associated user info compatible with the
175 // requested type (and throws an error if it is not the case)
176 //
177 // We can then access the "vertex_number" member of MyUserInfo to
178 // extract the targeted information.
179 return p.user_info<MyUserInfo>().vertex_number()==_vertex_number;
180 }
181
182private:
183 int _vertex_number; // the vertex number we're selecting
184};
185
186// The function that allows to write e.g.
187// Selector sel = !SelectorVertexNumber(0);
188// to select particles from all vertices except the 0th.
189Selector SelectorVertexNumber(const int & vertex_number){
190 return Selector(new SW_VertexNumber(vertex_number));
191}
192
193
194//------------------------------------------------------------------------
195// The example code associating user-info to the particles in the event
196int main(){
197 // read in input particles
198 //----------------------------------------------------------
199 vector<PseudoJet> input_particles;
200
201 double px, py , pz, E;
202 string str;
203 int vertex_number=-1;
204 int pdg_id = 21;
205 while (getline(cin, str)){
206 // if the line does not start with #, it's a particle
207 // read its momentum and pdg id
208 if (str[0] != '#'){
209 istringstream iss(str);
210 if (! (iss >> px >> py >> pz >> E >> pdg_id)){
211 cerr << "Wrong file format: particles must be specified with" << endl;
212 cerr << " px py pz E pdg_id" << endl;
213 }
214
215 // first create a PseudoJet with the correct momentum
216 PseudoJet p(px,py,pz,E);
217
218 // associate to that our user-defined extra information
219 // which is done using
220 // PseudoJet::set_user_info()
221 //
222 // IMPORTANT NOTE: set_user_info(...) takes a pointer as an
223 // argument. It will "own" that pointer i.e. will delete it when
224 // all the PseudoJet's using it will be deleted.
225 //
226 // NB: once you've done p.set_user_info(my_user_info_ptr), you must
227 // not call p2.set_user_info(my_user_info_ptr) with the same pointer
228 // because p and p2 will both attempt to delete it when they go out
229 // of scope causing a double-free corruption error. Instead do
230 // p2.user_info_shared_ptr() = p.user_info_shared_ptr();
231 p.set_user_info(new MyUserInfo(pdg_id, vertex_number));
232 PseudoJet p2; // defined only to make the above documentation consistent!
233
234 input_particles.push_back(p);
235 continue;
236 }
237
238 // check if we start a new sub-event
239 if ((str.length()>=9) && (str.compare(0,9,"#SUBSTART")==0)){
240 vertex_number++;
241 }
242 }
243
244
245 // create a jet definition:
246 // a jet algorithm with a given radius parameter
247 //----------------------------------------------------------
248 double R = 0.6;
250
251
252 // run the jet clustering with the above jet definition
253 //----------------------------------------------------------
254 ClusterSequence clust_seq(input_particles, jet_def);
255
256
257 // get the resulting jets ordered in pt
258 //----------------------------------------------------------
259 double ptmin = 25.0;
260 vector<PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(ptmin));
261
262
263 // tell the user what was done
264 // - the description of the algorithm used
265 // - extract the inclusive jets with pt > 5 GeV
266 // show the output as
267 // {index, rap, phi, pt}
268 //----------------------------------------------------------
269 cout << "Ran " << jet_def.description() << endl;
270
271 // label the columns
272 printf("%5s %15s %15s %15s %15s %15s\n","jet #",
273 "rapidity", "phi", "pt",
274 "pt_hard", "pt_pi0+gamma");
275
276 // a selection on the 1st vertex
277 Selector sel_vtx0 = SelectorVertexNumber(0);
278
279 // a selection on the pi0
280 Selector sel_pi0gamma = SelectorIsPi0Gamma();
281
282 // print out the details for each jet
283 for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
284 const PseudoJet & full = inclusive_jets[i];
285 const vector<PseudoJet> constituents = full.constituents();
286
287 // get the contribution from the 1st vertex
288 PseudoJet hard = join(sel_vtx0(constituents));
289
290 // get the contribution from the pi0's
291 PseudoJet pi0gamma = join(sel_pi0gamma(constituents));
292
293 // print the result
294 printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
295 full.rap(), full.phi(), full.perp(),
296 hard.perp(), pi0gamma.perp());
297 }
298
299 return 0;
300}
int main()
an example program showing how to use fastjet
Definition 01-basic.cc:52
deals with clustering
class that is intended to hold a full definition of the jet clusterer
a base class to hold extra user information in a PseudoJet
Definition PseudoJet.hh:423
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition PseudoJet.cc:681
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition PseudoJet.hh:138
const L & user_info() const
returns a reference to the dynamic cast conversion of user_info to type L.
Definition PseudoJet.hh:478
double phi() const
returns phi (in the range 0..2pi)
Definition PseudoJet.hh:123
double perp() const
returns the scalar transverse momentum
Definition PseudoJet.hh:158
default selector worker is an abstract virtual base class
Definition Selector.hh:59
virtual std::string description() const
returns a description of the worker
Definition Selector.hh:94
virtual bool pass(const PseudoJet &jet) const =0
returns true if a given object passes the selection criterion, and is the main function that needs to...
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition Selector.hh:149
the FastJet namespace
@ antikt_algorithm
like the k_t but with distance measures dij = min(1/kti^2,1/ktj^2) Delta R_{ij}^2 / R^2 diB = 1/kti^2
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition PseudoJet.cc:871