FastJet 3.0.2
|
00001 00002 //STARTHEADER 00003 // $Id: fastjet_subtraction.cc 2577 2011-09-13 15:11:38Z salam $ 00004 // 00005 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00006 // 00007 //---------------------------------------------------------------------- 00008 // This file is part of FastJet. 00009 // 00010 // FastJet is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU General Public License as published by 00012 // the Free Software Foundation; either version 2 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // The algorithms that underlie FastJet have required considerable 00016 // development and are described in hep-ph/0512210. If you use 00017 // FastJet as part of work towards a scientific publication, please 00018 // include a citation to the FastJet paper. 00019 // 00020 // FastJet is distributed in the hope that it will be useful, 00021 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00022 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00023 // GNU General Public License for more details. 00024 // 00025 // You should have received a copy of the GNU General Public License 00026 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00027 //---------------------------------------------------------------------- 00028 //ENDHEADER 00029 00030 00031 00032 00033 00034 00035 00036 //---------------------------------------------------------------------- 00037 // fastjet subtraction example program. 00038 // 00039 // Compile it with: make fastjet_subtraction 00040 // run it with : ./fastjet_subtraction < data/Pythia-dijet-ptmin100-lhc-pileup-1ev.dat 00041 // 00042 //---------------------------------------------------------------------- 00043 #include "fastjet/PseudoJet.hh" 00044 #include "fastjet/ClusterSequenceArea.hh" 00045 #include<iostream> // needed for io 00046 #include<sstream> // needed for internal io 00047 #include<vector> 00048 00049 using namespace std; 00050 00051 // A declaration of a function that pretty prints a list of jets 00052 // The subtraction is also performed inside this function 00053 void print_jets (const fastjet::ClusterSequenceAreaBase &, 00054 const vector<fastjet::PseudoJet> &); 00055 00056 /// an example program showing how to use fastjet 00057 int main (int argc, char ** argv) { 00058 00059 // since we use here simulated data we can split the hard event 00060 // from the full (i.e. with pileup added) one 00061 vector<fastjet::PseudoJet> hard_event, full_event; 00062 00063 // maximum rapidity that we accept for the input particles 00064 double etamax = 6.0; 00065 00066 // read in input particles. Keep the hard event generated by PYTHIA 00067 // separated from the full event, so as to be able to gauge the 00068 // "goodness" of the subtraction from the full event which also 00069 // includes pileup 00070 string line; 00071 int nsub = 0; 00072 while (getline(cin, line)) { 00073 istringstream linestream(line); 00074 // take substrings to avoid problems when there are extra "pollution" 00075 // characters (e.g. line-feed). 00076 if (line.substr(0,4) == "#END") {break;} 00077 if (line.substr(0,9) == "#SUBSTART") { 00078 // if more sub events follow, make copy of hard one here 00079 if (nsub == 1) hard_event = full_event; 00080 nsub += 1; 00081 } 00082 if (line.substr(0,1) == "#") {continue;} 00083 valarray<double> fourvec(4); 00084 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3]; 00085 fastjet::PseudoJet psjet(fourvec); 00086 psjet.set_user_index(0); 00087 // push event onto back of full_event vector 00088 if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);} 00089 } 00090 00091 // if we have read in only one event, copy it across here... 00092 if (nsub == 1) hard_event = full_event; 00093 00094 // if there was nothing in the event 00095 if (nsub == 0) { 00096 cerr << "Error: read empty event\n"; 00097 exit(-1); 00098 } 00099 00100 00101 // create an object that represents your choice of jet algorithm and 00102 // the associated parameters 00103 double R = 0.7; 00104 fastjet::Strategy strategy = fastjet::Best; 00105 fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R, strategy); 00106 //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, R, strategy); 00107 //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R, strategy); 00108 00109 // create an object that specifies how we to define the (active) area 00110 double ghost_etamax = 6.0; 00111 int active_area_repeats = 1; 00112 double ghost_area = 0.01; 00113 fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, 00114 ghost_area); 00115 fastjet::AreaDefinition area_def(ghost_spec); 00116 //fastjet::AreaDefinition area_def(fastjet::VoronoiAreaSpec(1.0)); 00117 00118 // run the jet clustering with the above jet definition. hard event first 00119 fastjet::ClusterSequenceArea clust_seq(hard_event, jet_def, area_def); 00120 //fastjet::ClusterSequencePassiveArea clust_seq(hard_event, jet_def); 00121 00122 00123 // extract the inclusive jets with pt > 5 GeV, sorted by pt 00124 double ptmin = 5.0; 00125 vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin); 00126 00127 // print them out 00128 cout << "" <<endl; 00129 cout << "Hard event only"<<endl; 00130 cout << "Number of input particles: "<<hard_event.size()<<endl; 00131 cout << "Strategy used: "<<clust_seq.strategy_string()<<endl; 00132 cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n"; 00133 cout << "---------------------------------------\n"; 00134 print_jets(clust_seq, inclusive_jets); 00135 cout << endl; 00136 00137 00138 // repeat everything on the full event 00139 00140 // run the jet clustering with the above jet definition 00141 fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def); 00142 00143 // extract the inclusive jets with pt > 20 GeV, sorted by pt 00144 ptmin = 20.0; 00145 inclusive_jets = clust_seq_full.inclusive_jets(ptmin); 00146 00147 // print them out 00148 cout << "Full event, with pileup, and its subtraction"<<endl; 00149 cout << "Number of input particles: "<<full_event.size()<<endl; 00150 cout << "Strategy used: "<<clust_seq_full.strategy_string()<<endl; 00151 cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV (before subtraction)\n"; 00152 cout << "---------------------------------------\n"; 00153 print_jets(clust_seq_full, inclusive_jets); 00154 cout << endl; 00155 00156 00157 } 00158 00159 00160 //---------------------------------------------------------------------- 00161 /// a function that pretty prints a list of jets, and performs the subtraction 00162 /// in two different ways, using a generic ClusterSequenceAreaBase 00163 /// type object. 00164 void print_jets (const fastjet::ClusterSequenceAreaBase & clust_seq, 00165 const vector<fastjet::PseudoJet> & unsorted_jets ) { 00166 00167 // sort jets into increasing pt 00168 vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets); 00169 00170 // the corrected jets will go in here 00171 vector<fastjet::PseudoJet> corrected_jets(jets.size()); 00172 00173 // get median pt per unit area -- you must decide the rapidity range 00174 // in which you measure the activity. (If doing a ghost-based active 00175 // area, make sure that your ghosts go up at least to \sim range+R). 00176 double range = 5.0; 00177 double median_pt_per_area = clust_seq.median_pt_per_unit_area(range); 00178 double median_pt_per_area4vector = clust_seq.median_pt_per_unit_area_4vector(range); 00179 00180 printf(" ijet rap phi Pt area Pt corr (rap corr phi corr Pt corr)ext\n"); 00181 for (size_t j = 0; j < jets.size(); j++) { 00182 00183 // get area of each jet 00184 double area = jets[j].area(); 00185 00186 // "standard" correction. Subtract only the Pt 00187 double pt_corr = jets[j].perp() - area*median_pt_per_area; 00188 00189 // "extended" correction 00190 fastjet::PseudoJet sub_4vect = 00191 median_pt_per_area4vector*jets[j].area_4vector(); 00192 if (sub_4vect.perp2() >= jets[j].perp2() || 00193 sub_4vect.E() >= jets[j].E()) { 00194 // if the correction is too large, set the jet to zero 00195 corrected_jets[j] = 0.0 * jets[j]; 00196 } else { 00197 // otherwise do an E-scheme subtraction 00198 corrected_jets[j] = jets[j] - sub_4vect; 00199 } 00200 // NB We could also leave out the above "if": 00201 // corrected_jets[j] = jets[j] - sub_4vect; 00202 // but the result would be different, since we would not avoid 00203 // jets with negative Pt or energy 00204 00205 printf("%5u %7.3f %7.3f %9.3f %7.3f %9.3f %7.3f %7.3f %9.3f\n", 00206 j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, pt_corr, 00207 corrected_jets[j].rap(),corrected_jets[j].phi(), corrected_jets[j].perp()); 00208 } 00209 00210 cout << endl; 00211 cout << "median pt_over_area = " << median_pt_per_area << endl; 00212 cout << "median pt_over_area4vector = " << median_pt_per_area4vector << endl << endl; 00213 00214 00215 } 00216 00217 00218 00219 00220 00221 00222