FastJet 3.0.2
|
00001 00002 //STARTHEADER 00003 // $Id: fastjet_areas.cc 2704 2011-11-16 11:11:10Z soyez $ 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 areas example program. 00038 // 00039 // Compile it with: make fastjet_areas 00040 // run it with : ./fastjet_areas < data/single-event.dat 00041 // 00042 //---------------------------------------------------------------------- 00043 #include "fastjet/PseudoJet.hh" 00044 #include "fastjet/ClusterSequenceArea.hh" 00045 #include "fastjet/ClusterSequencePassiveArea.hh" 00046 00047 // get info on how fastjet was configured 00048 #include "fastjet/config.h" 00049 00050 #ifdef ENABLE_PLUGIN_SISCONE 00051 #include "fastjet/SISConePlugin.hh" 00052 #endif 00053 00054 #include<iostream> // needed for io 00055 #include<sstream> // needed for internal io 00056 #include<vector> 00057 00058 using namespace std; 00059 00060 // a declaration of a function that pretty prints a list of jets 00061 void print_jets (const vector<fastjet::PseudoJet> &); 00062 00063 /// an example program showing how to use fastjet 00064 int main () { 00065 00066 vector<fastjet::PseudoJet> input_particles; 00067 00068 // read in input particles 00069 double px, py , pz, E; 00070 while (cin >> px >> py >> pz >> E) { 00071 // create a fastjet::PseudoJet with these components and put it onto 00072 // back of the input_particles vector 00073 input_particles.push_back(fastjet::PseudoJet(px,py,pz,E)); 00074 } 00075 00076 // create an object that represents your choice of jet algorithm, and 00077 // the associated parameters 00078 double Rparam = 1.0; 00079 fastjet::Strategy strategy = fastjet::Best; 00080 fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, strategy); 00081 //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, Rparam, strategy); 00082 //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, Rparam, strategy); 00083 //fastjet::JetDefinition jet_def(new fastjet::SISConePlugin(Rparam,0.75)); 00084 00085 // create an object that specifies how we to define the area 00086 fastjet::AreaDefinition area_def; 00087 bool use_voronoi = false; 00088 if (!use_voronoi) { 00089 double ghost_etamax = 6.0; 00090 double ghost_area = 0.01; 00091 int active_area_repeats = 1; 00092 00093 // now create the object that holds info about ghosts, and from that 00094 // get an area definition 00095 fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, 00096 ghost_area); 00097 area_def = fastjet::AreaDefinition(fastjet::active_area,ghost_spec); 00098 //area_def = fastjet::AreaDefinition(fastjet::passive_area,ghost_spec); 00099 } else { 00100 double effective_Rfact = 1.0; 00101 area_def = fastjet::VoronoiAreaSpec(effective_Rfact); 00102 } 00103 00104 // run the jet clustering with the above jet definition 00105 fastjet::ClusterSequenceArea clust_seq(input_particles, 00106 jet_def, area_def); 00107 // you can also run the individual area classes directly 00108 //fastjet::ClusterSequencePassiveArea clust_seq(input_particles, jet_def, 00109 // area_def.ghost_spec()); 00110 00111 // you may want to find out how much area in a given range (|y|<range) 00112 // is empty of real jets (or corresponds to pure "ghost" jets). 00113 //double range = 4.0; 00114 //cout << clust_seq.empty_area(range) << endl; 00115 //cout << clust_seq.n_empty_jets(range) << endl; 00116 00117 // tell the user what was done 00118 cout << "Jet definition was: " << jet_def.description() << endl; 00119 cout << "Area definition was: " << area_def.description() << endl; 00120 cout << "Strategy adopted by FastJet was "<< 00121 clust_seq.strategy_string()<<endl<<endl; 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 << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n"; 00129 cout << "---------------------------------------\n"; 00130 print_jets(inclusive_jets); 00131 cout << endl; 00132 00133 00134 cout << "Number of unclustered particles: " 00135 << clust_seq.unclustered_particles().size() << endl; 00136 00137 00138 } 00139 00140 00141 //---------------------------------------------------------------------- 00142 /// a function that pretty prints a list of jets 00143 void print_jets (const vector<fastjet::PseudoJet> & unsorted_jets) { 00144 00145 // sort jets into increasing pt 00146 vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets); 00147 00148 printf(" ijet rap phi Pt area +- err\n"); 00149 for (unsigned int j = 0; j < jets.size(); j++) { 00150 00151 double area = jets[j].area(); 00152 double area_error = jets[j].area_error(); 00153 00154 printf("%5u %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(), 00155 jets[j].phi(),jets[j].perp(), area, area_error); 00156 } 00157 00158 00159 } 00160 00161 00162 00163 00164 00165 00166