FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
07-subtraction-old.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example07old 07 - subtracting jet background contamination (old version)
4 ///
5 /// fastjet subtraction example program.
6 ///
7 /// Note that this example is deprecated --- see 07-subtraction.cc
8 /// for the newest version --- so it is not built by default
9 ///
10 /// run it with : ./07-subtraction-old < data/Pythia-Zp2jets-lhc-pileup-1ev.dat
11 ///
12 /// Source code: 07-subtraction-old.cc
13 //----------------------------------------------------------------------
14 
15 //STARTHEADER
16 // $Id: 07-subtraction-old.cc 2577 2011-09-13 15:11:38Z salam $
17 //
18 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
19 //
20 //----------------------------------------------------------------------
21 // This file is part of FastJet.
22 //
23 // FastJet is free software; you can redistribute it and/or modify
24 // it under the terms of the GNU General Public License as published by
25 // the Free Software Foundation; either version 2 of the License, or
26 // (at your option) any later version.
27 //
28 // The algorithms that underlie FastJet have required considerable
29 // development and are described in hep-ph/0512210. If you use
30 // FastJet as part of work towards a scientific publication, please
31 // include a citation to the FastJet paper.
32 //
33 // FastJet is distributed in the hope that it will be useful,
34 // but WITHOUT ANY WARRANTY; without even the implied warranty of
35 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 // GNU General Public License for more details.
37 //
38 // You should have received a copy of the GNU General Public License
39 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
40 //----------------------------------------------------------------------
41 //ENDHEADER
42 
43 #include "fastjet/PseudoJet.hh"
44 #include "fastjet/ClusterSequenceArea.hh"
45 #include <iostream> // needed for io
46 
47 using namespace std;
48 
49 int main (int argc, char ** argv) {
50 
51  // read in input particles
52  //
53  // since we use here simulated data we can split the hard event
54  // from the full (i.e. with pileup added) one
55  //----------------------------------------------------------
56 
57  vector<fastjet::PseudoJet> hard_event, full_event;
58 
59  // read in input particles. Keep the hard event generated by PYTHIA
60  // separated from the full event, so as to be able to gauge the
61  // "goodness" of the subtraction from the full event, which also
62  // includes pileup
63  double particle_maxrap = 5.0;
64 
65  string line;
66  int nsub = 0; // counter to keep track of which sub-event we're reading
67  while (getline(cin, line)) {
68  istringstream linestream(line);
69  // take substrings to avoid problems when there are extra "pollution"
70  // characters (e.g. line-feed).
71  if (line.substr(0,4) == "#END") {break;}
72  if (line.substr(0,9) == "#SUBSTART") {
73  // if more sub events follow, make copy of first one (the hard one) here
74  if (nsub == 1) hard_event = full_event;
75  nsub += 1;
76  }
77  if (line.substr(0,1) == "#") {continue;}
78  double px,py,pz,E;
79  linestream >> px >> py >> pz >> E;
80  // you can construct
81  fastjet::PseudoJet particle(px,py,pz,E);
82 
83  // push event onto back of full_event vector
84  if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
85  }
86 
87  // if we have read in only one event, copy it across here...
88  if (nsub == 1) hard_event = full_event;
89 
90  // if there was nothing in the event
91  if (nsub == 0) {
92  cerr << "Error: read empty event\n";
93  exit(-1);
94  }
95 
96 
97  // create a jet definition for the clustering
98  // We use the anti-kt algorithm with a radius of 0.5
99  //----------------------------------------------------------
100  double R = 0.5;
102 
103  // create an area definition for the clustering
104  //----------------------------------------------------------
105  // ghosts should go up to the acceptance of the detector or
106  // (with infinite acceptance) at least 2R beyond the region
107  // where you plan to investigate jets.
108  double ghost_maxrap = 6.0;
109  fastjet::GhostedAreaSpec area_spec(ghost_maxrap);
110  fastjet::AreaDefinition area_def(fastjet::active_area, area_spec);
111 
112  // run the jet clustering with the above jet and area definitions
113  // for both the hard and full event
114  //
115  // We retrieve the jets above 7 GeV in both case (note that the
116  // 7-GeV cut we be applied again later on after we subtract the jets
117  // from the full event)
118  // ----------------------------------------------------------
119  fastjet::ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def);
120  fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
121 
122  double ptmin = 7.0;
123  vector<fastjet::PseudoJet> hard_jets = sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin));
124  vector<fastjet::PseudoJet> full_jets = sorted_by_pt(clust_seq_full.inclusive_jets(ptmin));
125 
126  // Now turn to the estimation of the background (for the full event)
127  //
128  // This also requires a ClusterSequenceArea.
129  // In general, this ClusterSequenceArea does not need to be the same
130  // as the one used (above) to cluster and extract the jets from the
131  // event:
132  // - We strongly recommend using the kt or Cambridge/Aachen algorithm
133  // (a warning will be issued otherwise)
134  // - The choice of the radius is a bit more subtle. R=0.4 has been
135  // chosen to limit the impact of hard jets; in samples of
136  // dominantly sparse events it may cause the UE/pileup to be
137  // underestimated a little, a slightly larger value (0.5 or 0.6)
138  // may be better.
139  // - For the area definition, we recommend the use of explicit
140  // ghosts (i.e. active_area_explicit_ghosts)
141  // As mentionned in the area example (06-area.cc), ghosts should
142  // extend sufficiently far in rapidity to cover the jets used in
143  // the computation of the background (see also the comment below)
144  //
145  // ----------------------------------------------------------
147  fastjet::GhostedAreaSpec area_spec_bkgd(ghost_maxrap);
148  fastjet::AreaDefinition area_def_bkgd(fastjet::active_area_explicit_ghosts, area_spec_bkgd);
149  fastjet::ClusterSequenceArea clust_seq_bkgd(full_event, jet_def_bkgd, area_def_bkgd);
150 
151  // Once you have the ClusterSequenceArea, you can compute the
152  // background. This is estimated over a given range
153  // (RangeDefinition) i.e. only jets within that range will be used
154  // to estimate the background. You shold thus make sure the ghosts
155  // extend far enough in rapidity to cover the range, a warning will
156  // be issued otherwise.
157  //
158  // The simplest way to define a RangeDefinition is through its
159  // maximal |y| extent but other options are possible e.g. through a
160  // minimal and maximal rapidity and minimal and maximal azimuthal
161  // angle. If needed, you can even define your own ranges (a few are
162  // provided with FastJet)
163  //
164  // Finally, the estimation of the background properties rho (the
165  // average density per unit area) and sigma (the average
166  // fluctuations per unit area) is done using
167  // ClusterSequenceArea::get_median_rho_and_sigma(). This takes
168  // 2 main parameters: the range discussed above and a boolean
169  // controlling the use of 4-vector or scalar areas (we suggest using
170  // 4-vector areas)
171  //
172  // ----------------------------------------------------------
173  double range_maxrap = 4.5; // we have a ghost_maxrap of 6.0, particles up to 5
174  fastjet::RangeDefinition range(range_maxrap);
175 
176  bool use_4vector_area = true;
177 
178  double rho, sigma;
179  clust_seq_bkgd.get_median_rho_and_sigma(range, use_4vector_area, rho, sigma);
180 
181  // show a summary of what was done so far
182  // - the description of the algorithms, areas and ranges used
183  // - the background properties
184  // - the jets in the hard event
185  //----------------------------------------------------------
186  cout << "Main clustering:" << endl;
187  cout << " Ran: " << jet_def.description() << endl;
188  cout << " Area: " << area_def.description() << endl;
189  cout << " Particles up to |y|=" << particle_maxrap << endl;
190  cout << endl;
191 
192  cout << "Background estimation:" << endl;
193  cout << " Ran " << jet_def_bkgd.description() << endl;
194  cout << " Area: " << area_def_bkgd.description() << endl;
195  cout << " Range: " << range.description() << endl;
196  cout << " Giving, for the full event" << endl;
197  cout << " rho = " << rho << endl;
198  cout << " sigma = " << sigma << endl;
199  cout << endl;
200 
201  cout << "Jets above " << ptmin << " GeV in the hard event (" << hard_event.size() << " particles)" << endl;
202  cout << "---------------------------------------\n";
203  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area");
204  for (unsigned int i = 0; i < hard_jets.size(); i++) {
205  printf("%5u %15.8f %15.8f %15.8f %15.8f\n", i,
206  hard_jets[i].rap(), hard_jets[i].phi(), hard_jets[i].perp(),
207  hard_jets[i].area());
208  }
209  cout << endl;
210 
211  // Once the background properties have been computed, subtraction
212  // can be applied on the jets
213  //
214  // This uses ClusterSequenceArea::subtracted_jet(jet, rho), with the
215  // ClusterSequence used to cluster the jet and the background
216  // density we have just computed
217  //
218  // (Note that when using scalar areas, subtracted_pt should be used
219  // instead of subtracted_jet)
220  //
221  // We output the jets before and after subtraction
222  // ----------------------------------------------------------
223  cout << "Jets above " << ptmin << " GeV in the full event (" << full_event.size() << " particles)" << endl;
224  cout << "---------------------------------------\n";
225  printf("%5s %15s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area", "rap_sub", "phi_sub", "pt_sub");
226  unsigned int idx=0;
227  for (unsigned int i=0; i<full_jets.size(); i++){
228  // get the subtracted jet
229  fastjet::PseudoJet subtracted_jet = clust_seq_full.subtracted_jet(full_jets[i], rho);
230 
231  // re-apply the pt cut
232  if (subtracted_jet.perp2() >= ptmin*ptmin){
233  printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", idx,
234  full_jets[i].rap(), full_jets[i].phi(), full_jets[i].perp(),
235  full_jets[i].area(),
236  subtracted_jet.rap(), subtracted_jet.phi(),
237  subtracted_jet.perp());
238  idx++;
239  }
240  }
241 
242  return 0;
243 }
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:123
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:799
General class for user to obtain ClusterSequence with additional area information.
the longitudinally invariant kt algorithm
class that holds a generic area definition
class for holding a range definition specification, given by limits on rapidity and azimuth...
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:108
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...
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:143
Parameters to configure the computation of jet areas using ghosts.
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:141
class that is intended to hold a full definition of the jet clusterer