FastJet  3.0.3
 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 }