43 #include "fastjet/PseudoJet.hh"
44 #include "fastjet/ClusterSequenceArea.hh"
54 const vector<fastjet::PseudoJet> &);
57 int main (
int argc,
char ** argv) {
61 vector<fastjet::PseudoJet> hard_event, full_event;
72 while (getline(cin, line)) {
73 istringstream linestream(line);
76 if (line.substr(0,4) ==
"#END") {
break;}
77 if (line.substr(0,9) ==
"#SUBSTART") {
79 if (nsub == 1) hard_event = full_event;
82 if (line.substr(0,1) ==
"#") {
continue;}
83 valarray<double> fourvec(4);
84 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
86 psjet.set_user_index(0);
88 if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);}
92 if (nsub == 1) hard_event = full_event;
96 cerr <<
"Error: read empty event\n";
110 double ghost_etamax = 6.0;
111 int active_area_repeats = 1;
112 double ghost_area = 0.01;
125 vector<fastjet::PseudoJet> inclusive_jets = clust_seq.
inclusive_jets(ptmin);
129 cout <<
"Hard event only"<<endl;
130 cout <<
"Number of input particles: "<<hard_event.size()<<endl;
132 cout <<
"Printing inclusive jets with pt > "<< ptmin<<
" GeV\n";
133 cout <<
"---------------------------------------\n";
145 inclusive_jets = clust_seq_full.inclusive_jets(ptmin);
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";
165 const vector<fastjet::PseudoJet> & unsorted_jets ) {
168 vector<fastjet::PseudoJet> jets =
sorted_by_pt(unsorted_jets);
171 vector<fastjet::PseudoJet> corrected_jets(jets.size());
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++) {
184 double area = jets[j].area();
187 double pt_corr = jets[j].perp() - area*median_pt_per_area;
192 if (sub_4vect.
perp2() >= jets[j].perp2() ||
193 sub_4vect.E() >= jets[j].E()) {
195 corrected_jets[j] = 0.0 * jets[j];
198 corrected_jets[j] = jets[j] - sub_4vect;
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());
211 cout <<
"median pt_over_area = " << median_pt_per_area << endl;
212 cout <<
"median pt_over_area4vector = " << median_pt_per_area4vector << endl << endl;