fastjet_areas.cc

Go to the documentation of this file.
00001 
00002 //STARTHEADER
00003 // $Id: fastjet_areas.cc 958 2007-11-28 17:43:48Z cacciari $
00004 //
00005 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00027 //  Foundation, Inc.:
00028 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00029 //----------------------------------------------------------------------
00030 //ENDHEADER
00031 
00032 
00033 
00034 
00035 
00036 
00037 
00038 //----------------------------------------------------------------------
00039 // fastjet areas example program. 
00040 //
00041 // Compile it with: make fastjet_areas
00042 // run it with    : ./fastjet_areas < data/single-event.dat
00043 //
00044 //----------------------------------------------------------------------
00045 #include "fastjet/PseudoJet.hh"
00046 #include "fastjet/ClusterSequenceArea.hh"
00047 #include "fastjet/ClusterSequencePassiveArea.hh"
00048 
00049 // get info on how fastjet was configured
00050 #include "fastjet/config.h"
00051 
00052 #ifdef ENABLE_PLUGIN_SISCONE
00053 #include "fastjet/SISConePlugin.hh"
00054 #endif
00055 
00056 #include<iostream> // needed for io
00057 #include<sstream>  // needed for internal io
00058 #include<vector> 
00059 
00060 using namespace std;
00061 
00062 // a declaration of a function that pretty prints a list of jets
00063 void print_jets (const fastjet::ClusterSequenceAreaBase &, 
00064                  const vector<fastjet::PseudoJet> &);
00065 
00067 int main (int argc, char ** argv) {
00068   
00069   vector<fastjet::PseudoJet> input_particles;
00070   
00071   // read in input particles
00072   double px, py , pz, E;
00073   while (cin >> px >> py >> pz >> E) {
00074     // create a fastjet::PseudoJet with these components and put it onto
00075     // back of the input_particles vector
00076     input_particles.push_back(fastjet::PseudoJet(px,py,pz,E)); 
00077   }
00078 
00079   // create an object that represents your choice of jet algorithm, and 
00080   // the associated parameters
00081   double Rparam = 1.0;
00082   fastjet::Strategy strategy = fastjet::Best;
00083   fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, strategy);
00084   //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, Rparam, strategy);
00085   //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, Rparam, strategy);
00086   //fastjet::JetDefinition jet_def(new fastjet::SISConePlugin(1.0));
00087 
00088   // create an object that specifies how we to define the area
00089   fastjet::AreaDefinition area_def;
00090   bool use_voronoi = false;
00091   if (!use_voronoi) {
00092     double ghost_etamax = 7.0;
00093     int    active_area_repeats = 1;
00094     //double ghost_area    = 0.01;
00095     //int    active_area_repeats = 100;
00096     double ghost_area    = 0.05;
00097     //area_def = fastjet::GhostedAreaSpec(ghost_etamax, active_area_repeats, 
00098     fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, 
00099                                         ghost_area);
00100     area_def = fastjet::AreaDefinition(fastjet::passive_area,ghost_spec);
00101   } else {
00102     double effective_Rfact = 1.0;
00103     area_def = fastjet::VoronoiAreaSpec(effective_Rfact);
00104   }
00105 
00106   // run the jet clustering with the above jet definition
00107   fastjet::ClusterSequenceArea clust_seq(input_particles, 
00108                                              jet_def, area_def);
00109   // run the jet clustering with the above jet definition
00110   //fastjet::ClusterSequencePassiveArea clust_seq(input_particles, jet_def, 
00111   //                                              area_def.active_spec());
00112 
00113   //cout << clust_seq.empty_area(4.0) << endl;
00114   //cout << clust_seq.n_empty_jets(4.0) << endl;
00115 
00116   // tell the user what was done
00117   cout << "Jet definition was: " << jet_def.description() << endl;
00118   cout << "Area definition was: " << area_def.description() << endl;
00119   cout << "Strategy adopted by FastJet was "<<
00120        clust_seq.strategy_string()<<endl<<endl;
00121 
00122   // extract the inclusive jets with pt > 5 GeV, sorted by pt
00123   double ptmin = 5.0;
00124   vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
00125 
00126   // print them out
00127   cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n";
00128   cout << "---------------------------------------\n";
00129   print_jets(clust_seq, inclusive_jets);
00130   cout << endl;
00131 
00132   
00133   cout << "Number of unclustered particles: " 
00134        << clust_seq.unclustered_particles().size() << endl;
00135 
00136 
00137 }
00138 
00139 
00140 //----------------------------------------------------------------------
00142 void print_jets (const fastjet::ClusterSequenceAreaBase & clust_seq, 
00143                  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 (size_t j = 0; j < jets.size(); j++) {
00150 
00151     double area     = clust_seq.area(jets[j]);
00152     double area_error = clust_seq.area_error(jets[j]);
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 

Generated on Tue Dec 18 17:05:02 2007 for fastjet by  doxygen 1.5.2