FastJet  3.3.0
07-subtraction.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example07 07 - subtracting jet background contamination
4 ///
5 /// fastjet subtraction example program.
6 ///
7 /// run it with : ./07-subtraction < data/Pythia-Zp2jets-lhc-pileup-1ev.dat
8 ///
9 /// Source code: 07-subtraction.cc
10 //----------------------------------------------------------------------
11 
12 //STARTHEADER
13 // $Id: 07-subtraction.cc 3670 2014-09-08 14:17:59Z soyez $
14 //
15 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
16 //
17 //----------------------------------------------------------------------
18 // This file is part of FastJet.
19 //
20 // FastJet is free software; you can redistribute it and/or modify
21 // it under the terms of the GNU General Public License as published by
22 // the Free Software Foundation; either version 2 of the License, or
23 // (at your option) any later version.
24 //
25 // The algorithms that underlie FastJet have required considerable
26 // development and are described in hep-ph/0512210. If you use
27 // FastJet as part of work towards a scientific publication, please
28 // include a citation to the FastJet paper.
29 //
30 // FastJet is distributed in the hope that it will be useful,
31 // but WITHOUT ANY WARRANTY; without even the implied warranty of
32 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 // GNU General Public License for more details.
34 //
35 // You should have received a copy of the GNU General Public License
36 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
37 //----------------------------------------------------------------------
38 //ENDHEADER
39 
40 #include "fastjet/PseudoJet.hh"
41 #include "fastjet/ClusterSequenceArea.hh"
42 #include "fastjet/Selector.hh"
43 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
44 #include "fastjet/tools/Subtractor.hh"
45 #include <iostream> // needed for io
46 #include <fastjet/config.h> // for the FASTJET_VERSION_NUMBER preprocessor symbol
47 
48 using namespace std;
49 using namespace fastjet;
50 
51 int main(){
52 
53  // read in input particles
54  //
55  // since we use here simulated data we can split the hard event
56  // from the full (i.e. with pileup added) one
57  //----------------------------------------------------------
58 
59  vector<PseudoJet> hard_event, full_event;
60 
61  // read in input particles. Keep the hard event generated by PYTHIA
62  // separated from the full event, so as to be able to gauge the
63  // "goodness" of the subtraction from the full event, which also
64  // includes pileup
65  double particle_maxrap = 5.0;
66 
67  string line;
68  int nsub = 0; // counter to keep track of which sub-event we're reading
69  while (getline(cin, line)) {
70  istringstream linestream(line);
71  // take substrings to avoid problems when there are extra "pollution"
72  // characters (e.g. line-feed).
73  if (line.substr(0,4) == "#END") {break;}
74  if (line.substr(0,9) == "#SUBSTART") {
75  // if more sub events follow, make copy of first one (the hard one) here
76  if (nsub == 1) hard_event = full_event;
77  nsub += 1;
78  }
79  if (line.substr(0,1) == "#") {continue;}
80  double px,py,pz,E;
81  linestream >> px >> py >> pz >> E;
82  // you can construct
83  PseudoJet particle(px,py,pz,E);
84 
85  // push event onto back of full_event vector
86  if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
87  }
88 
89  // if we have read in only one event, copy it across here...
90  if (nsub == 1) hard_event = full_event;
91 
92  // if there was nothing in the event
93  if (nsub == 0) {
94  cerr << "Error: read empty event\n";
95  exit(-1);
96  }
97 
98 
99  // create a jet definition for the clustering
100  // We use the anti-kt algorithm with a radius of 0.5
101  //----------------------------------------------------------
102  double R = 0.5;
103  JetDefinition jet_def(antikt_algorithm, R);
104 
105  // create an area definition for the clustering
106  //----------------------------------------------------------
107  // ghosts should go up to the acceptance of the detector or
108  // (with infinite acceptance) at least 2R beyond the region
109  // where you plan to investigate jets.
110  double ghost_maxrap = 6.0;
111  GhostedAreaSpec area_spec(ghost_maxrap);
112  AreaDefinition area_def(active_area, area_spec);
113 
114  // run the jet clustering with the above jet and area definitions
115  // for both the hard and full event
116  //
117  // We retrieve the jets above 7 GeV in both case (note that the
118  // 7-GeV cut we be applied again later on after we subtract the jets
119  // from the full event)
120  // ----------------------------------------------------------
121  ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def);
122  ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
123 
124  double ptmin = 7.0;
125  vector<PseudoJet> hard_jets = sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin));
126  vector<PseudoJet> full_jets = sorted_by_pt(clust_seq_full.inclusive_jets(ptmin));
127 
128  // Now turn to the estimation of the background (for the full event)
129  //
130  // There are different ways to do that. In general, this also
131  // requires clustering the particles that will be handled internally
132  // in FastJet.
133  //
134  // The suggested way to proceed is to use a BackgroundEstimator
135  // constructed from the following 3 arguments:
136  // - a jet definition used to cluster the particles.
137  // . We strongly recommend using the kt or Cambridge/Aachen
138  // algorithm (a warning will be issued otherwise)
139  // . The choice of the radius is a bit more subtle. R=0.4 has
140  // been chosen to limit the impact of hard jets; in samples of
141  // dominantly sparse events it may cause the UE/pileup to be
142  // underestimated a little, a slightly larger value (0.5 or
143  // 0.6) may be better.
144  // - An area definition for which we recommend the use of explicit
145  // ghosts (i.e. active_area_explicit_ghosts)
146  // As mentionned in the area example (06-area.cc), ghosts should
147  // extend sufficiently far in rapidity to cover the jets used in
148  // the computation of the background (see also the comment below)
149  // - A Selector specifying the range over which we will keep the
150  // jets entering the estimation of the background (you should
151  // thus make sure the ghosts extend far enough in rapidity to
152  // cover the range, a warning will be issued otherwise).
153  // In this particular example, the two hardest jets in the event
154  // are removed from the background estimation
155  // ----------------------------------------------------------
156  JetDefinition jet_def_bkgd(kt_algorithm, 0.4);
157  AreaDefinition area_def_bkgd(active_area_explicit_ghosts,
158  GhostedAreaSpec(ghost_maxrap));
159  Selector selector = SelectorAbsRapMax(4.5) * (!SelectorNHardest(2));
160  JetMedianBackgroundEstimator bkgd_estimator(selector, jet_def_bkgd, area_def_bkgd);
161 
162  // To help manipulate the background estimator, we also provide a
163  // transformer that allows to apply directly the background
164  // subtraction on the jets. This will use the background estimator
165  // to compute rho for the jets to be subtracted.
166  // ----------------------------------------------------------
167  Subtractor subtractor(&bkgd_estimator);
168 
169  // since FastJet 3.1.0, rho_m is supported natively in background
170  // estimation (both JetMedianBackgroundEstimator and
171  // GridMedianBackgroundEstimator).
172  //
173  // For backward-compatibility reasons its use is by default switched off
174  // (as is the enforcement of m>0 for the subtracted jets). The
175  // following 2 lines of code switch these on. They are strongly
176  // recommended and should become the default in future versions of
177  // FastJet.
178  //
179  // Note that we also illustrate the use of the
180  // FASTJET_VERSION_NUMBER macro
181 #if FASTJET_VERSION_NUMBER >= 30100
182  subtractor.set_use_rho_m(true);
183  subtractor.set_safe_mass(true);
184 #endif
185 
186  // Finally, once we have an event, we can just tell the background
187  // estimator to use that list of particles
188  // This could be done directly when declaring the background
189  // estimator but the usage below can more easily be accomodated to a
190  // loop over a set of events.
191  // ----------------------------------------------------------
192  bkgd_estimator.set_particles(full_event);
193 
194  // show a summary of what was done so far
195  // - the description of the algorithms, areas and ranges used
196  // - the background properties
197  // - the jets in the hard event
198  //----------------------------------------------------------
199  cout << "Main clustering:" << endl;
200  cout << " Ran: " << jet_def.description() << endl;
201  cout << " Area: " << area_def.description() << endl;
202  cout << " Particles up to |y|=" << particle_maxrap << endl;
203  cout << endl;
204 
205  cout << "Background estimation:" << endl;
206  cout << " " << bkgd_estimator.description() << endl << endl;;
207  cout << " Giving, for the full event" << endl;
208  cout << " rho = " << bkgd_estimator.rho() << endl;
209  cout << " sigma = " << bkgd_estimator.sigma() << endl;
210 #if FASTJET_VERSION_NUMBER >= 30100
211  cout << " rho_m = " << bkgd_estimator.rho_m() << endl;
212  cout << " sigma_m = " << bkgd_estimator.sigma_m() << endl;
213 #endif
214  cout << endl;
215 
216  cout << "Jets above " << ptmin << " GeV in the hard event (" << hard_event.size() << " particles)" << endl;
217  cout << "---------------------------------------\n";
218  printf("%5s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "m", "area");
219  for (unsigned int i = 0; i < hard_jets.size(); i++) {
220  printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
221  hard_jets[i].rap(), hard_jets[i].phi(), hard_jets[i].pt(), hard_jets[i].m(),
222  hard_jets[i].area());
223  }
224  cout << endl;
225 
226  // Once the background properties have been computed, subtraction
227  // can be applied on the jets. Subtraction is performed on the
228  // full 4-vector
229  //
230  // We output the jets before and after subtraction
231  // ----------------------------------------------------------
232  cout << "Jets above " << ptmin << " GeV in the full event (" << full_event.size() << " particles)" << endl;
233  cout << "---------------------------------------\n";
234  printf("%5s %15s %15s %15s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "m", "area", "rap_sub", "phi_sub", "pt_sub", "m_sub");
235  unsigned int idx=0;
236 
237  // get the subtracted jets
238  vector<PseudoJet> subtracted_jets = subtractor(full_jets);
239 
240  for (unsigned int i=0; i<full_jets.size(); i++){
241  // re-apply the pt cut
242  if (subtracted_jets[i].pt2() >= ptmin*ptmin){
243  printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", idx,
244  full_jets[i].rap(), full_jets[i].phi(), full_jets[i].pt(), full_jets[i].m(),
245  full_jets[i].area(),
246  subtracted_jets[i].rap(), subtracted_jets[i].phi(),
247  subtracted_jets[i].pt(),
248  subtracted_jets[i].m());
249  idx++;
250  }
251  }
252 
253  return 0;
254 }
Class that helps perform jet background subtraction.
Definition: Subtractor.hh:62
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:774
Class to estimate the pt density of the background per unit area, using the median of the distributio...
General class for user to obtain ClusterSequence with additional area information.
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Definition: Selector.cc:1074
Selector SelectorAbsRapMax(double absrapmax)
select objects with |rap| <= absrapmax
Definition: Selector.cc:880
class that holds a generic area definition
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
the FastJet namespace
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
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
class that is intended to hold a full definition of the jet clusterer