FastJet  3.0.6
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 2684 2011-11-14 07:41:44Z 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 
47 using namespace std;
48 using namespace fastjet;
49 
50 int main(){
51 
52  // read in input particles
53  //
54  // since we use here simulated data we can split the hard event
55  // from the full (i.e. with pileup added) one
56  //----------------------------------------------------------
57 
58  vector<PseudoJet> hard_event, full_event;
59 
60  // read in input particles. Keep the hard event generated by PYTHIA
61  // separated from the full event, so as to be able to gauge the
62  // "goodness" of the subtraction from the full event, which also
63  // includes pileup
64  double particle_maxrap = 5.0;
65 
66  string line;
67  int nsub = 0; // counter to keep track of which sub-event we're reading
68  while (getline(cin, line)) {
69  istringstream linestream(line);
70  // take substrings to avoid problems when there are extra "pollution"
71  // characters (e.g. line-feed).
72  if (line.substr(0,4) == "#END") {break;}
73  if (line.substr(0,9) == "#SUBSTART") {
74  // if more sub events follow, make copy of first one (the hard one) here
75  if (nsub == 1) hard_event = full_event;
76  nsub += 1;
77  }
78  if (line.substr(0,1) == "#") {continue;}
79  double px,py,pz,E;
80  linestream >> px >> py >> pz >> E;
81  // you can construct
82  PseudoJet particle(px,py,pz,E);
83 
84  // push event onto back of full_event vector
85  if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
86  }
87 
88  // if we have read in only one event, copy it across here...
89  if (nsub == 1) hard_event = full_event;
90 
91  // if there was nothing in the event
92  if (nsub == 0) {
93  cerr << "Error: read empty event\n";
94  exit(-1);
95  }
96 
97 
98  // create a jet definition for the clustering
99  // We use the anti-kt algorithm with a radius of 0.5
100  //----------------------------------------------------------
101  double R = 0.5;
102  JetDefinition jet_def(antikt_algorithm, R);
103 
104  // create an area definition for the clustering
105  //----------------------------------------------------------
106  // ghosts should go up to the acceptance of the detector or
107  // (with infinite acceptance) at least 2R beyond the region
108  // where you plan to investigate jets.
109  double ghost_maxrap = 6.0;
110  GhostedAreaSpec area_spec(ghost_maxrap);
111  AreaDefinition area_def(active_area, area_spec);
112 
113  // run the jet clustering with the above jet and area definitions
114  // for both the hard and full event
115  //
116  // We retrieve the jets above 7 GeV in both case (note that the
117  // 7-GeV cut we be applied again later on after we subtract the jets
118  // from the full event)
119  // ----------------------------------------------------------
120  ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def);
121  ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
122 
123  double ptmin = 7.0;
124  vector<PseudoJet> hard_jets = sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin));
125  vector<PseudoJet> full_jets = sorted_by_pt(clust_seq_full.inclusive_jets(ptmin));
126 
127  // Now turn to the estimation of the background (for the full event)
128  //
129  // There are different ways to do that. In general, this also
130  // requires clustering the particles that will be handled internally
131  // in FastJet.
132  //
133  // The suggested way to proceed is to use a BackgroundEstimator
134  // constructed from the following 3 arguments:
135  // - a jet definition used to cluster the particles.
136  // . We strongly recommend using the kt or Cambridge/Aachen
137  // algorithm (a warning will be issued otherwise)
138  // . The choice of the radius is a bit more subtle. R=0.4 has
139  // been chosen to limit the impact of hard jets; in samples of
140  // dominantly sparse events it may cause the UE/pileup to be
141  // underestimated a little, a slightly larger value (0.5 or
142  // 0.6) may be better.
143  // - An area definition for which we recommend the use of explicit
144  // ghosts (i.e. active_area_explicit_ghosts)
145  // As mentionned in the area example (06-area.cc), ghosts should
146  // extend sufficiently far in rapidity to cover the jets used in
147  // the computation of the background (see also the comment below)
148  // - A Selector specifying the range over which we will keep the
149  // jets entering the estimation of the background (you should
150  // thus make sure the ghosts extend far enough in rapidity to
151  // cover the range, a warning will be issued otherwise).
152  // In this particular example, the two hardest jets in the event
153  // are removed from the background estimation
154  // ----------------------------------------------------------
155  JetDefinition jet_def_bkgd(kt_algorithm, 0.4);
156  AreaDefinition area_def_bkgd(active_area_explicit_ghosts,
157  GhostedAreaSpec(ghost_maxrap));
158  Selector selector = SelectorAbsRapMax(4.5) * (!SelectorNHardest(2));
159  JetMedianBackgroundEstimator bkgd_estimator(selector, jet_def_bkgd, area_def_bkgd);
160 
161  // To help manipulate the background estimator, we also provide a
162  // transformer that allows to apply directly the background
163  // subtraction on the jets. This will use the background estimator
164  // to compute rho for the jets to be subtracted.
165  // ----------------------------------------------------------
166  Subtractor subtractor(&bkgd_estimator);
167 
168  // Finally, once we have an event, we can just tell the background
169  // estimator to use that list of particles
170  // This could be done directly when declaring the background
171  // estimator but the usage below can more easily be accomodated to a
172  // loop over a set of events.
173  // ----------------------------------------------------------
174  bkgd_estimator.set_particles(full_event);
175 
176  // show a summary of what was done so far
177  // - the description of the algorithms, areas and ranges used
178  // - the background properties
179  // - the jets in the hard event
180  //----------------------------------------------------------
181  cout << "Main clustering:" << endl;
182  cout << " Ran: " << jet_def.description() << endl;
183  cout << " Area: " << area_def.description() << endl;
184  cout << " Particles up to |y|=" << particle_maxrap << endl;
185  cout << endl;
186 
187  cout << "Background estimation:" << endl;
188  cout << " " << bkgd_estimator.description() << endl << endl;;
189  cout << " Giving, for the full event" << endl;
190  cout << " rho = " << bkgd_estimator.rho() << endl;
191  cout << " sigma = " << bkgd_estimator.sigma() << endl;
192  cout << endl;
193 
194  cout << "Jets above " << ptmin << " GeV in the hard event (" << hard_event.size() << " particles)" << endl;
195  cout << "---------------------------------------\n";
196  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area");
197  for (unsigned int i = 0; i < hard_jets.size(); i++) {
198  printf("%5u %15.8f %15.8f %15.8f %15.8f\n", i,
199  hard_jets[i].rap(), hard_jets[i].phi(), hard_jets[i].perp(),
200  hard_jets[i].area());
201  }
202  cout << endl;
203 
204  // Once the background properties have been computed, subtraction
205  // can be applied on the jets. Subtraction is performed on the
206  // full 4-vector
207  //
208  // We output the jets before and after subtraction
209  // ----------------------------------------------------------
210  cout << "Jets above " << ptmin << " GeV in the full event (" << full_event.size() << " particles)" << endl;
211  cout << "---------------------------------------\n";
212  printf("%5s %15s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area", "rap_sub", "phi_sub", "pt_sub");
213  unsigned int idx=0;
214 
215  // get the subtracted jets
216  vector<PseudoJet> subtracted_jets = subtractor(full_jets);
217 
218  for (unsigned int i=0; i<full_jets.size(); i++){
219  // re-apply the pt cut
220  if (subtracted_jets[i].perp2() >= ptmin*ptmin){
221  printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", idx,
222  full_jets[i].rap(), full_jets[i].phi(), full_jets[i].perp(),
223  full_jets[i].area(),
224  subtracted_jets[i].rap(), subtracted_jets[i].phi(),
225  subtracted_jets[i].perp());
226  idx++;
227  }
228  }
229 
230  return 0;
231 }
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:785
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:1022
Class that helps perform jet background subtraction.
Definition: Subtractor.hh:59
Selector SelectorAbsRapMax(double absrapmax)
select objects with |rap| &lt;= absrapmax
Definition: Selector.cc:828
class that holds a generic area definition
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
std::string description() const
return a textual description of the current jet definition
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
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:65
class that is intended to hold a full definition of the jet clusterer