00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "fastjet/PseudoJet.hh"
00046 #include "fastjet/ClusterSequenceActiveArea.hh"
00047 #include<iostream>
00048 #include<sstream>
00049 #include<vector>
00050
00051 using namespace std;
00052
00053
00054
00055 void print_jets (const fastjet::ClusterSequenceActiveArea &,
00056 const vector<fastjet::PseudoJet> &);
00057
00059 int main (int argc, char ** argv) {
00060
00061
00062
00063 vector<fastjet::PseudoJet> hard_event, full_event;
00064
00065
00066 double etamax = 6.0;
00067
00068
00069
00070
00071
00072 string line;
00073 int nsub = 0;
00074 while (getline(cin, line)) {
00075 istringstream linestream(line);
00076
00077
00078 if (line.substr(0,4) == "#END") {break;}
00079 if (line.substr(0,9) == "#SUBSTART") {
00080
00081 if (nsub == 1) hard_event = full_event;
00082 nsub += 1;
00083 }
00084 if (line.substr(0,1) == "#") {continue;}
00085 valarray<double> fourvec(4);
00086 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00087 fastjet::PseudoJet psjet(fourvec);
00088 psjet.set_user_index(0);
00089
00090 if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);}
00091 }
00092
00093
00094 if (nsub == 1) hard_event = full_event;
00095
00096
00097 if (nsub == 0) {
00098 cerr << "Error: read empty event\n";
00099 exit(-1);
00100 }
00101
00102
00103
00104
00105 double R = 0.7;
00106 fastjet::Strategy strategy = fastjet::Best;
00107 fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R, strategy);
00108
00109
00110 double ghost_etamax = 6.0;
00111 int active_area_repeats = 1;
00112 double ghost_area = 0.01;
00113 fastjet::ActiveAreaSpec area_spec(ghost_etamax, active_area_repeats,
00114 ghost_area);
00115
00116
00117 fastjet::ClusterSequenceActiveArea clust_seq(hard_event,
00118 jet_def, area_spec);
00119
00120
00121
00122 double ptmin = 5.0;
00123 vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
00124
00125
00126 cout << "" <<endl;
00127 cout << "Hard event only"<<endl;
00128 cout << "Number of input particles: "<<hard_event.size()<<endl;
00129 cout << "Strategy used: "<<clust_seq.strategy_string()<<endl;
00130 cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n";
00131 cout << "---------------------------------------\n";
00132 print_jets(clust_seq, inclusive_jets);
00133 cout << endl;
00134
00135
00136
00137
00138
00139 fastjet::ClusterSequenceActiveArea clust_seq_full(full_event,
00140 jet_def, area_spec);
00141
00142
00143 ptmin = 20.0;
00144 inclusive_jets = clust_seq_full.inclusive_jets(ptmin);
00145
00146
00147 cout << "Full event, with pileup, and its subtraction"<<endl;
00148 cout << "Number of input particles: "<<full_event.size()<<endl;
00149 cout << "Strategy used: "<<clust_seq_full.strategy_string()<<endl;
00150 cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV (before subtraction)\n";
00151 cout << "---------------------------------------\n";
00152 print_jets(clust_seq_full, inclusive_jets);
00153 cout << endl;
00154
00155
00156 }
00157
00158
00159
00162 void print_jets (const fastjet::ClusterSequenceActiveArea & clust_seq,
00163 const vector<fastjet::PseudoJet> & unsorted_jets ) {
00164
00165
00166 vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets);
00167
00168
00169 vector<fastjet::PseudoJet> corrected_jets(jets.size());
00170
00171
00172
00173
00174
00175 double median_pt_per_area = clust_seq.pt_per_unit_area();
00176
00177 printf(" ijet rap phi Pt area Pt corr (rap corr phi corr Pt corr)ext\n");
00178 for (size_t j = 0; j < jets.size(); j++) {
00179
00180
00181 double area = clust_seq.area(jets[j]);
00182
00183
00184 double pt_corr = jets[j].perp() - area*median_pt_per_area;
00185
00186
00187 fastjet::PseudoJet area_4vect =
00188 median_pt_per_area*clust_seq.area_4vector(jets[j]);
00189 if (area_4vect.perp2() >= jets[j].perp2() ||
00190 area_4vect.E() >= jets[j].E()) {
00191
00192 corrected_jets[j] = 0.0 * jets[j];
00193 } else {
00194
00195 corrected_jets[j] = jets[j] - area_4vect;
00196 }
00197
00198
00199
00200
00201
00202 printf("%5u %7.3f %7.3f %9.3f %7.3f %9.3f %7.3f %7.3f %9.3f\n",
00203 j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, pt_corr,
00204 corrected_jets[j].rap(),corrected_jets[j].phi(), corrected_jets[j].perp());
00205 }
00206
00207 cout << endl;
00208 cout << "median pt_over_area = " << median_pt_per_area << endl << endl;
00209
00210
00211 }
00212
00213
00214
00215
00216
00217
00218