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