FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
fastjet_subtraction.cc
1 
2 //STARTHEADER
3 // $Id: fastjet_subtraction.cc 2577 2011-09-13 15:11:38Z salam $
4 //
5 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
6 //
7 //----------------------------------------------------------------------
8 // This file is part of FastJet.
9 //
10 // FastJet is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU General Public License as published by
12 // the Free Software Foundation; either version 2 of the License, or
13 // (at your option) any later version.
14 //
15 // The algorithms that underlie FastJet have required considerable
16 // development and are described in hep-ph/0512210. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // include a citation to the FastJet paper.
19 //
20 // FastJet is distributed in the hope that it will be useful,
21 // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 // GNU General Public License for more details.
24 //
25 // You should have received a copy of the GNU General Public License
26 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
27 //----------------------------------------------------------------------
28 //ENDHEADER
29 
30 
31 
32 
33 
34 
35 
36 //----------------------------------------------------------------------
37 // fastjet subtraction example program.
38 //
39 // Compile it with: make fastjet_subtraction
40 // run it with : ./fastjet_subtraction < data/Pythia-dijet-ptmin100-lhc-pileup-1ev.dat
41 //
42 //----------------------------------------------------------------------
43 #include "fastjet/PseudoJet.hh"
44 #include "fastjet/ClusterSequenceArea.hh"
45 #include<iostream> // needed for io
46 #include<sstream> // needed for internal io
47 #include<vector>
48 
49 using namespace std;
50 
51 // A declaration of a function that pretty prints a list of jets
52 // The subtraction is also performed inside this function
54  const vector<fastjet::PseudoJet> &);
55 
56 /// an example program showing how to use fastjet
57 int main (int argc, char ** argv) {
58 
59  // since we use here simulated data we can split the hard event
60  // from the full (i.e. with pileup added) one
61  vector<fastjet::PseudoJet> hard_event, full_event;
62 
63  // maximum rapidity that we accept for the input particles
64  double etamax = 6.0;
65 
66  // read in input particles. Keep the hard event generated by PYTHIA
67  // separated from the full event, so as to be able to gauge the
68  // "goodness" of the subtraction from the full event which also
69  // includes pileup
70  string line;
71  int nsub = 0;
72  while (getline(cin, line)) {
73  istringstream linestream(line);
74  // take substrings to avoid problems when there are extra "pollution"
75  // characters (e.g. line-feed).
76  if (line.substr(0,4) == "#END") {break;}
77  if (line.substr(0,9) == "#SUBSTART") {
78  // if more sub events follow, make copy of hard one here
79  if (nsub == 1) hard_event = full_event;
80  nsub += 1;
81  }
82  if (line.substr(0,1) == "#") {continue;}
83  valarray<double> fourvec(4);
84  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
85  fastjet::PseudoJet psjet(fourvec);
86  psjet.set_user_index(0);
87  // push event onto back of full_event vector
88  if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);}
89  }
90 
91  // if we have read in only one event, copy it across here...
92  if (nsub == 1) hard_event = full_event;
93 
94  // if there was nothing in the event
95  if (nsub == 0) {
96  cerr << "Error: read empty event\n";
97  exit(-1);
98  }
99 
100 
101  // create an object that represents your choice of jet algorithm and
102  // the associated parameters
103  double R = 0.7;
104  fastjet::Strategy strategy = fastjet::Best;
105  fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R, strategy);
106  //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, R, strategy);
107  //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R, strategy);
108 
109  // create an object that specifies how we to define the (active) area
110  double ghost_etamax = 6.0;
111  int active_area_repeats = 1;
112  double ghost_area = 0.01;
113  fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats,
114  ghost_area);
115  fastjet::AreaDefinition area_def(ghost_spec);
116  //fastjet::AreaDefinition area_def(fastjet::VoronoiAreaSpec(1.0));
117 
118  // run the jet clustering with the above jet definition. hard event first
119  fastjet::ClusterSequenceArea clust_seq(hard_event, jet_def, area_def);
120  //fastjet::ClusterSequencePassiveArea clust_seq(hard_event, jet_def);
121 
122 
123  // extract the inclusive jets with pt > 5 GeV, sorted by pt
124  double ptmin = 5.0;
125  vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
126 
127  // print them out
128  cout << "" <<endl;
129  cout << "Hard event only"<<endl;
130  cout << "Number of input particles: "<<hard_event.size()<<endl;
131  cout << "Strategy used: "<<clust_seq.strategy_string()<<endl;
132  cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n";
133  cout << "---------------------------------------\n";
134  print_jets(clust_seq, inclusive_jets);
135  cout << endl;
136 
137 
138  // repeat everything on the full event
139 
140  // run the jet clustering with the above jet definition
141  fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
142 
143  // extract the inclusive jets with pt > 20 GeV, sorted by pt
144  ptmin = 20.0;
145  inclusive_jets = clust_seq_full.inclusive_jets(ptmin);
146 
147  // print them out
148  cout << "Full event, with pileup, and its subtraction"<<endl;
149  cout << "Number of input particles: "<<full_event.size()<<endl;
150  cout << "Strategy used: "<<clust_seq_full.strategy_string()<<endl;
151  cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV (before subtraction)\n";
152  cout << "---------------------------------------\n";
153  print_jets(clust_seq_full, inclusive_jets);
154  cout << endl;
155 
156 
157 }
158 
159 
160 //----------------------------------------------------------------------
161 /// a function that pretty prints a list of jets, and performs the subtraction
162 /// in two different ways, using a generic ClusterSequenceAreaBase
163 /// type object.
164 void print_jets (const fastjet::ClusterSequenceAreaBase & clust_seq,
165  const vector<fastjet::PseudoJet> & unsorted_jets ) {
166 
167  // sort jets into increasing pt
168  vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets);
169 
170  // the corrected jets will go in here
171  vector<fastjet::PseudoJet> corrected_jets(jets.size());
172 
173  // get median pt per unit area -- you must decide the rapidity range
174  // in which you measure the activity. (If doing a ghost-based active
175  // area, make sure that your ghosts go up at least to \sim range+R).
176  double range = 5.0;
177  double median_pt_per_area = clust_seq.median_pt_per_unit_area(range);
178  double median_pt_per_area4vector = clust_seq.median_pt_per_unit_area_4vector(range);
179 
180  printf(" ijet rap phi Pt area Pt corr (rap corr phi corr Pt corr)ext\n");
181  for (size_t j = 0; j < jets.size(); j++) {
182 
183  // get area of each jet
184  double area = jets[j].area();
185 
186  // "standard" correction. Subtract only the Pt
187  double pt_corr = jets[j].perp() - area*median_pt_per_area;
188 
189  // "extended" correction
190  fastjet::PseudoJet sub_4vect =
191  median_pt_per_area4vector*jets[j].area_4vector();
192  if (sub_4vect.perp2() >= jets[j].perp2() ||
193  sub_4vect.E() >= jets[j].E()) {
194  // if the correction is too large, set the jet to zero
195  corrected_jets[j] = 0.0 * jets[j];
196  } else {
197  // otherwise do an E-scheme subtraction
198  corrected_jets[j] = jets[j] - sub_4vect;
199  }
200  // NB We could also leave out the above "if":
201  // corrected_jets[j] = jets[j] - sub_4vect;
202  // but the result would be different, since we would not avoid
203  // jets with negative Pt or energy
204 
205  printf("%5u %7.3f %7.3f %9.3f %7.3f %9.3f %7.3f %7.3f %9.3f\n",
206  j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, pt_corr,
207  corrected_jets[j].rap(),corrected_jets[j].phi(), corrected_jets[j].perp());
208  }
209 
210  cout << endl;
211  cout << "median pt_over_area = " << median_pt_per_area << endl;
212  cout << "median pt_over_area4vector = " << median_pt_per_area4vector << endl << endl;
213 
214 
215 }
216 
217 
218 
219 
220 
221 
222 
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.
double median_pt_per_unit_area(const Selector &selector) const
the median of (pt/area) for jets contained within the selector range, making use also of the info on ...
the longitudinally invariant kt algorithm
class that holds a generic area definition
double median_pt_per_unit_area_4vector(const Selector &selector) const
the median of (pt/area_4vector) for jets contained within the selector range, making use also of the ...
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
base class that sets interface for extensions of ClusterSequence that provide information about the a...
void print_jets(const vector< fastjet::PseudoJet > &)
a function that pretty prints a list of jets
automatic selection of the best (based on N), including the LazyTiled strategies that are new to FJ3...
Parameters to configure the computation of jet areas using ghosts.
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
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
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:732