FastJet 3.0.2
|
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 }