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"
58 vector<PseudoJet> hard_event, full_event;
64 double particle_maxrap = 5.0;
68 while (getline(cin, line)) {
69 istringstream linestream(line);
72 if (line.substr(0,4) ==
"#END") {
break;}
73 if (line.substr(0,9) ==
"#SUBSTART") {
75 if (nsub == 1) hard_event = full_event;
78 if (line.substr(0,1) ==
"#") {
continue;}
80 linestream >> px >> py >> pz >> E;
85 if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
89 if (nsub == 1) hard_event = full_event;
93 cerr <<
"Error: read empty event\n";
109 double ghost_maxrap = 6.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));
174 bkgd_estimator.set_particles(full_event);
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;
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;
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());
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");
216 vector<PseudoJet> subtracted_jets = subtractor(full_jets);
218 for (
unsigned int i=0; i<full_jets.size(); i++){
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(),
224 subtracted_jets[i].rap(), subtracted_jets[i].phi(),
225 subtracted_jets[i].perp());