fastjet 2.4.5
fastjet_subjets.cc
Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: fastjet_subjets.cc 1555 2009-05-08 11:48:38Z cacciari $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 
00032 //----------------------------------------------------------------------
00033 // fastjet example program to show how to access subjets;
00034 // 
00035 // See also fastjet_higgs_decomp.cc to see the use of subjets for
00036 // identifying boosted higgs (and other objects)
00037 //
00038 // Compile it with: make fastjet_subjets
00039 // run it with    : ./fastjet_subjets < data/single-event.dat
00040 //
00041 //----------------------------------------------------------------------
00042 #include "fastjet/PseudoJet.hh"
00043 #include "fastjet/ClusterSequence.hh"
00044 #include<iostream> // needed for io
00045 #include<sstream>  // needed for internal io
00046 #include<vector> 
00047 #include <cstdio>
00048 
00049 using namespace std;
00050 
00051 // to avoid excessive typing, define an abbreviation for the
00052 // fastjet namespace
00053 namespace fj = fastjet;
00054 
00055 
00056 // a declaration of a function that pretty prints a list of jets
00057 void print_jets (const fj::ClusterSequence &, 
00058                  const vector<fj::PseudoJet> &);
00059 
00060 // and this pretty prinst a single jet
00061 void print_jet (const fj::ClusterSequence & clust_seq, 
00062                 const fj::PseudoJet & jet);
00063 
00064 // pretty print the jets and their subjets
00065 void print_jets_and_sub (const fj::ClusterSequence & clust_seq, 
00066                          const vector<fj::PseudoJet> & jets,
00067                          double dcut);
00068 
00070 int main (int argc, char ** argv) {
00071   
00072   vector<fj::PseudoJet> input_particles;
00073   
00074   // read in input particles
00075   double px, py , pz, E;
00076   while (cin >> px >> py >> pz >> E) {
00077     // create a fj::PseudoJet with these components and put it onto
00078     // back of the input_particles vector
00079     input_particles.push_back(fj::PseudoJet(px,py,pz,E)); 
00080   }
00081   
00082   // We'll do two subjet analyses, physically rather different
00083   double R = 1.0;
00084   //fj::JetDefinition kt_def(fj::kt_algorithm, R);
00085   fj::JetDefinition cam_def(fj::cambridge_algorithm, R);
00086 
00087   // run the jet clustering with the above jet definition
00088   //fj::ClusterSequence kt_seq(input_particles, kt_def);
00089   fj::ClusterSequence cam_seq(input_particles, cam_def);
00090 
00091   // tell the user what was done
00092   cout << "Ran " << cam_def.description() << endl;
00093   cout << "Strategy adopted by FastJet was "<<
00094        cam_seq.strategy_string()<<endl<<endl;
00095 
00096   // extract the inclusive jets with pt > 5 GeV
00097   double ptmin = 5.0;
00098   vector<fj::PseudoJet> inclusive_jets = cam_seq.inclusive_jets(ptmin);
00099 
00100   // for the subjets  
00101   double smallR = 0.4;
00102   double dcut_cam = pow(smallR/R,2);
00103 
00104   // print them out
00105   cout << "Printing inclusive jets (R = "<<R<<") with pt > "<< ptmin<<" GeV\n";
00106   cout << "and their subjets with smallR = " << smallR << "\n";
00107   cout << "---------------------------------------\n";
00108   print_jets_and_sub(cam_seq, inclusive_jets, dcut_cam);
00109   cout << endl;
00110 
00111 
00112   // print them out
00113   vector<fj::PseudoJet> exclusive_jets = cam_seq.exclusive_jets(dcut_cam);
00114   cout << "Printing exclusive jets with dcut = "<< dcut_cam<<" \n";
00115   cout << "--------------------------------------------\n";
00116   print_jets(cam_seq, exclusive_jets);
00117 
00118 
00119 }
00120 
00121 
00122 //----------------------------------------------------------------------
00124 void print_jets (const fj::ClusterSequence & clust_seq, 
00125                  const vector<fj::PseudoJet> & jets) {
00126 
00127   // sort jets into increasing pt
00128   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00129 
00130   // label the columns
00131   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", 
00132          "phi", "pt", "n constituents");
00133   
00134   // print out the details for each jet
00135   for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00136     printf("%5u ",i);
00137     print_jet(clust_seq, sorted_jets[i]);
00138   }
00139 
00140 }
00141 
00142 
00143 //----------------------------------------------------------------------
00145 void print_jets_and_sub (const fj::ClusterSequence & clust_seq, 
00146                          const vector<fj::PseudoJet> & jets,
00147                          double dcut) {
00148 
00149   // sort jets into increasing pt
00150   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00151 
00152   // label the columns
00153   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", 
00154          "phi", "pt", "n constituents");
00155   
00156   // print out the details for each jet
00157   for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00158     printf("%5u       ",i);
00159     print_jet(clust_seq, sorted_jets[i]);
00160     vector<fj::PseudoJet> subjets = 
00161       sorted_by_pt(clust_seq.exclusive_subjets(sorted_jets[i], dcut));
00162     for (unsigned int j = 0; j < subjets.size(); j++) {
00163       printf("    -sub-%02u ",j);
00164       print_jet(clust_seq, subjets[j]);
00165     }
00166   }
00167 
00168 }
00169 
00170 
00171 //----------------------------------------------------------------------
00173 void print_jet (const fj::ClusterSequence & clust_seq, 
00174                 const fj::PseudoJet & jet) {
00175   int n_constituents = clust_seq.constituents(jet).size();
00176   printf("%15.8f %15.8f %15.8f %8u\n",
00177          jet.rap(), jet.phi(), jet.perp(), n_constituents);
00178 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines