FastJet 3.0.6
fastjet_timing.cc
00001 //STARTHEADER
00002 // $Id: fastjet_timing.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_timing.cc: Program to help time and test the fastjet package
00032 /// 
00033 /// It reads files containing multiple events in the format 
00034 /// p1x p1y p1z E1
00035 /// p2x p2y p2z E2
00036 /// ...
00037 /// #END
00038 /// 
00039 /// An example input file containing 10 events is included as 
00040 /// data/Pythia-PtMin1000-LHC-10ev.dat
00041 ///
00042 /// Usage:
00043 ///   fastjet_timing [-strategy NUMBER] [-repeat nrepeats] [-massive] \
00044 ///                  [-combine nevents] [-r Rparameter] [-incl ptmin] [...] \
00045 ///                  < data_file
00046 ///
00047 /// where the clustering can be repeated to aid timing and multiple
00048 /// events can be combined to get to larger multiplicities. Some options:
00049 ///
00050 ///   -strategy N   indicate stratgey from the enum FjStrategy (see
00051 ///                 FjClusterSequence.hh).
00052 ///
00053 ///   -combine nev  for combining multiple events from the data file in order
00054 ///                 to get to large multiplicities.
00055 ///
00056 ///   -incl ptmin   output of all inclusive jets with pt > ptmin is obtained
00057 ///                 with the -incl option.
00058 ///
00059 ///   -excld dcut   output of all exclusive jets as obtained in a clustering
00060 ///                 with dcut
00061 ///
00062 ///   -massless     read in only the 3-momenta and deduce energies assuming
00063 ///                 that particles are massless
00064 ///
00065 ///   -write        for writing out detailed clustering sequence (valuable
00066 ///                 for testing purposes)
00067 ///
00068 ///   -unique_write writes out the sequence of dij's according to the
00069 ///                 "unique_history_order" (useful for verifying consistency
00070 ///                 between different clustering strategies).
00071 ///
00072 ///   -cam          switch to the inclusive Cambridge/Aachen algorithm --
00073 ///                 note that the option -excld dcut provides a clustering
00074 ///                 up to the dcut which is the minimum squared
00075 ///                 distance between any pair of jets.
00076 ///
00077 #include "fastjet/PseudoJet.hh"
00078 #include "fastjet/ClusterSequence.hh"
00079 #include<iostream>
00080 #include<sstream>
00081 #include<valarray>
00082 #include<vector>
00083 #include <cstdio>
00084 #include <cstdlib>
00085 #include<cstddef> // for size_t
00086 #include "CmdLine.hh"
00087 
00088 using namespace std;
00089 
00090 // to avoid excessive typing, define an abbreviation for the 
00091 // fastjet namespace
00092 namespace fj = fastjet;
00093 
00094 inline double pow2(const double x) {return x*x;}
00095 
00096 /// a program to test and time the kt algorithm as implemented in fastjet
00097 int main (int argc, char ** argv) {
00098 
00099   CmdLine cmdline(argc,argv);
00100   // allow the use to specify the fj::Strategy either through the
00101   // -clever or the -strategy options (both will take numerical
00102   // values); the latter will override the former.
00103   fj::Strategy  strategy  = fj::Strategy(cmdline.int_val("-strategy",
00104                                         cmdline.int_val("-clever", fj::Best)));
00105   int  repeat  = cmdline.int_val("-repeat",1);
00106   int  combine = cmdline.int_val("-combine",1);
00107   bool write   = cmdline.present("-write");
00108   bool unique_write = cmdline.present("-unique_write");
00109   bool hydjet  = cmdline.present("-hydjet");
00110   double ktR   = cmdline.double_val("-r",1.0);
00111   double inclkt = cmdline.double_val("-incl",-1.0);
00112   int    excln  = cmdline.int_val   ("-excln",-1);
00113   double excld  = cmdline.double_val("-excld",-1.0);
00114   double etamax = cmdline.double_val("-etamax",1.0e305);
00115   bool   show_constituents = cmdline.present("-const");
00116   bool   massless = cmdline.present("-massless");
00117   int  nev     = cmdline.int_val("-nev",1);
00118   bool add_dense_coverage = cmdline.present("-dense");
00119 
00120   // The following option causes the Cambridge algo to be used.
00121   // Note that currently the only output that works sensibly here is
00122   // "-incl 0"
00123   fj::JetAlgorithm jet_algorithm;
00124   if (cmdline.present("-cam")) {
00125     jet_algorithm = fj::cambridge_algorithm;
00126   } else {
00127     jet_algorithm = fj::kt_algorithm;
00128   }
00129 
00130   if (!cmdline.all_options_used()) {cerr << 
00131       "Error: some options were not recognized"<<endl; 
00132     exit(-1);}
00133 
00134 
00135   for (int iev = 0; iev < nev; iev++) {
00136   vector<fj::PseudoJet> jets;
00137   string line;
00138   int  ndone = 0;
00139   while (getline(cin, line)) {
00140       //cout << line<<endl;
00141     istringstream linestream(line);
00142     if (line == "#END") {
00143       ndone += 1;
00144       if (ndone == combine) {break;}
00145     }
00146     if (line.substr(0,1) == "#") {continue;}
00147     valarray<double> fourvec(4);
00148     if (hydjet) {
00149       // special reading from hydjet.txt event record (though actually
00150       // this is supposed to be a standard pythia event record, so
00151       // being able to read from it is perhaps not so bad an idea...)
00152       int ii, istat,id,m1,m2,d1,d2;
00153       double mass;
00154       linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00155                  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00156       // current file contains mass of particle as 4th entry
00157       if (istat == 1) {
00158         fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00159                           +pow2(fourvec[2])+pow2(mass));
00160       }
00161     } else {
00162       if (massless) {
00163         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00164         fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00165       else {
00166         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00167       }
00168     }
00169     fj::PseudoJet psjet(fourvec);
00170     if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00171   }
00172 
00173   // add a fake underlying event which is very soft, uniformly distributed
00174   // in eta,phi so as to allow one to reconstruct the area that is associated
00175   // with each jet.
00176   if (add_dense_coverage) {
00177     srand(2);
00178     int nphi = 60;
00179     int neta = 100;
00180     double kt = 1e-1;
00181     for (int iphi = 0; iphi<nphi; iphi++) {
00182       for (int ieta = -neta; ieta<neta+1; ieta++) {
00183         double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00184         double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00185         kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00186         double pminus = kt*exp(-eta);
00187         double pplus  = kt*exp(+eta);
00188         double px = kt*sin(phi);
00189         double py = kt*cos(phi);
00190         //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00191         fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00192         jets.push_back(mom);
00193       }
00194     }
00195   }
00196   
00197   fj::JetDefinition jet_def(jet_algorithm, ktR, strategy);
00198 
00199   for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00200     fj::ClusterSequence clust_seq(jets,jet_def,write);
00201     if (irepeat != 0) {continue;}
00202     cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00203     cout << "strategy used =  "<< clust_seq.strategy_string()<< endl;
00204 
00205     // now provide some nice output...
00206     if (inclkt >= 0.0) {
00207       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00208       for (size_t j = 0; j < jets.size(); j++) {
00209         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00210         if (show_constituents) {
00211           vector<fj::PseudoJet> const_jets = jets[j].constituents();
00212           for (size_t k = 0; k < const_jets.size(); k++) {
00213             printf("        jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00214                    const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00215           }
00216           cout << "\n\n";
00217         }
00218       }
00219     }
00220 
00221     if (excln > 0) {
00222       vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00223  
00224       cout << "Printing "<<excln<<" exclusive jets\n";
00225       for (size_t j = 0; j < jets.size(); j++) {
00226         printf("%5u %15.8f %15.8f %15.8f\n",
00227                //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00228                j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00229       }
00230     }
00231 
00232     if (excld > 0.0) {
00233       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00234       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00235       for (size_t j = 0; j < jets.size(); j++) {
00236         printf("%5u %15.8f %15.8f %15.8f\n",
00237                j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00238       }
00239     }
00240     
00241     // useful for testing that recombination sequences are unique
00242     if (unique_write) {
00243       vector<int> unique_history = clust_seq.unique_history_order();
00244       // construct the inverse of the above mapping
00245       vector<int> inv_unique_history(clust_seq.history().size());
00246       for (unsigned int i = 0; i < unique_history.size(); i++) {
00247         inv_unique_history[unique_history[i]] = i;}
00248 
00249       for (unsigned int i = 0; i < unique_history.size(); i++) {
00250         fj::ClusterSequence::history_element el = 
00251           clust_seq.history()[unique_history[i]];
00252         int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00253         int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00254         printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00255       }
00256     }
00257 
00258   } // irepeat
00259 
00260   } // iev
00261 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends