43 #include "fastjet/PseudoJet.hh"
44 #include "fastjet/ClusterSequenceArea.hh"
49 int main (
int argc,
char ** argv) {
57 vector<fastjet::PseudoJet> hard_event, full_event;
63 double particle_maxrap = 5.0;
67 while (getline(cin, line)) {
68 istringstream linestream(line);
71 if (line.substr(0,4) ==
"#END") {
break;}
72 if (line.substr(0,9) ==
"#SUBSTART") {
74 if (nsub == 1) hard_event = full_event;
77 if (line.substr(0,1) ==
"#") {
continue;}
79 linestream >> px >> py >> pz >> E;
84 if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
88 if (nsub == 1) hard_event = full_event;
92 cerr <<
"Error: read empty event\n";
108 double ghost_maxrap = 6.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));
173 double range_maxrap = 4.5;
176 bool use_4vector_area =
true;
179 clust_seq_bkgd.get_median_rho_and_sigma(range, use_4vector_area, rho, sigma);
186 cout <<
"Main clustering:" << endl;
188 cout <<
" Area: " << area_def.description() << endl;
189 cout <<
" Particles up to |y|=" << particle_maxrap << endl;
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;
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());
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");
227 for (
unsigned int i=0; i<full_jets.size(); i++){
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(),
236 subtracted_jet.
rap(), subtracted_jet.
phi(),
237 subtracted_jet.
perp());