FastJet 3.0.4
fastjet_subjets.cc
00001 //STARTHEADER
00002 // $Id: fastjet_subjets.cc 2577 2011-09-13 15:11:38Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 //----------------------------------------------------------------------
00031 // fastjet example program to show how to access subjets;
00032 // 
00033 // See also fastjet_higgs_decomp.cc to see the use of subjets for
00034 // identifying boosted higgs (and other objects)
00035 //
00036 // Compile it with: make fastjet_subjets
00037 // run it with    : ./fastjet_subjets < data/single-event.dat
00038 //
00039 //----------------------------------------------------------------------
00040 #include "fastjet/PseudoJet.hh"
00041 #include "fastjet/ClusterSequence.hh"
00042 #include<iostream> // needed for io
00043 #include<sstream>  // needed for internal io
00044 #include<vector> 
00045 #include <cstdio>
00046 
00047 using namespace std;
00048 
00049 // to avoid excessive typing, define an abbreviation for the
00050 // fastjet namespace
00051 namespace fj = fastjet;
00052 
00053 
00054 // a declaration of a function that pretty prints a list of jets
00055 void print_jets (const vector<fj::PseudoJet> &);
00056 
00057 // and this pretty prinst a single jet
00058 void print_jet (const fj::PseudoJet & jet);
00059 
00060 // pretty print the jets and their subjets
00061 void print_jets_and_sub (const vector<fj::PseudoJet> & jets,
00062                          double dcut);
00063 
00064 /// an example program showing how to use fastjet
00065 int main (int argc, char ** argv) {
00066   
00067   vector<fj::PseudoJet> input_particles;
00068   
00069   // read in input particles
00070   double px, py , pz, E;
00071   while (cin >> px >> py >> pz >> E) {
00072     // create a fj::PseudoJet with these components and put it onto
00073     // back of the input_particles vector
00074     input_particles.push_back(fj::PseudoJet(px,py,pz,E)); 
00075   }
00076   
00077   // We'll do two subjet analyses, physically rather different
00078   double R = 1.0;
00079   //fj::JetDefinition kt_def(fj::kt_algorithm, R);
00080   fj::JetDefinition cam_def(fj::cambridge_algorithm, R);
00081 
00082   // run the jet clustering with the above jet definition
00083   //fj::ClusterSequence kt_seq(input_particles, kt_def);
00084   fj::ClusterSequence cam_seq(input_particles, cam_def);
00085 
00086   // tell the user what was done
00087   cout << "Ran " << cam_def.description() << endl;
00088   cout << "Strategy adopted by FastJet was "<<
00089        cam_seq.strategy_string()<<endl<<endl;
00090 
00091   // extract the inclusive jets with pt > 5 GeV
00092   double ptmin = 5.0;
00093   vector<fj::PseudoJet> inclusive_jets = cam_seq.inclusive_jets(ptmin);
00094 
00095   // for the subjets  
00096   double smallR = 0.4;
00097   double dcut_cam = pow(smallR/R,2);
00098 
00099   // print them out
00100   cout << "Printing inclusive jets (R = "<<R<<") with pt > "<< ptmin<<" GeV\n";
00101   cout << "and their subjets with smallR = " << smallR << "\n";
00102   cout << "---------------------------------------\n";
00103   print_jets_and_sub(inclusive_jets, dcut_cam);
00104   cout << endl;
00105 
00106 
00107   // print them out
00108   vector<fj::PseudoJet> exclusive_jets = cam_seq.exclusive_jets(dcut_cam);
00109   cout << "Printing exclusive jets with dcut = "<< dcut_cam<<" \n";
00110   cout << "--------------------------------------------\n";
00111   print_jets(exclusive_jets);
00112 
00113 
00114 }
00115 
00116 
00117 //----------------------------------------------------------------------
00118 /// a function that pretty prints a list of jets
00119 void print_jets (const vector<fj::PseudoJet> & jets) {
00120 
00121   // sort jets into increasing pt
00122   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00123 
00124   // label the columns
00125   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", 
00126          "phi", "pt", "n constituents");
00127   
00128   // print out the details for each jet
00129   for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00130     printf("%5u ",i);
00131     print_jet(sorted_jets[i]);
00132   }
00133 
00134 }
00135 
00136 
00137 //----------------------------------------------------------------------
00138 /// a function that pretty prints a list of jets
00139 void print_jets_and_sub (const vector<fj::PseudoJet> & jets,
00140                          double dcut) {
00141 
00142   // sort jets into increasing pt
00143   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00144 
00145   // label the columns
00146   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", 
00147          "phi", "pt", "n constituents");
00148   
00149   // print out the details for each jet
00150   for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00151     printf("%5u       ",i);
00152     print_jet(sorted_jets[i]);
00153     vector<fj::PseudoJet> subjets = sorted_by_pt(sorted_jets[i].exclusive_subjets(dcut));
00154     for (unsigned int j = 0; j < subjets.size(); j++) {
00155       printf("    -sub-%02u ",j);
00156       print_jet(subjets[j]);
00157     }
00158   }
00159 
00160 }
00161 
00162 
00163 //----------------------------------------------------------------------
00164 /// print a single jet
00165 void print_jet (const fj::PseudoJet & jet) {
00166   int n_constituents = jet.constituents().size();
00167   printf("%15.8f %15.8f %15.8f %8u\n",
00168          jet.rap(), jet.phi(), jet.perp(), n_constituents);
00169 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends