fastjet_timing.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: fastjet_timing.cc 958 2007-11-28 17:43:48Z 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 //----------------------------------------------------------------------
00079 #include "fastjet/PseudoJet.hh"
00080 #include "fastjet/ClusterSequence.hh"
00081 #include<iostream>
00082 #include<sstream>
00083 #include<valarray>
00084 #include<vector>
00085 #include <cstdlib>
00086 #include<cstddef> // for size_t
00087 #include "CmdLine.hh"
00088 
00089 using namespace std;
00090 
00091 // to avoid excessive typing, define an abbreviation for the 
00092 // fastjet namespace
00093 namespace fj = fastjet;
00094 
00095 inline double pow2(const double x) {return x*x;}
00096 
00098 int main (int argc, char ** argv) {
00099 
00100   CmdLine cmdline(argc,argv);
00101   // allow the use to specify the fj::Strategy either through the
00102   // -clever or the -strategy options (both will take numerical
00103   // values); the latter will override the former.
00104   fj::Strategy  strategy  = fj::Strategy(cmdline.int_val("-strategy",
00105                                         cmdline.int_val("-clever", fj::Best)));
00106   int  repeat  = cmdline.int_val("-repeat",1);
00107   int  combine = cmdline.int_val("-combine",1);
00108   bool write   = cmdline.present("-write");
00109   bool unique_write = cmdline.present("-unique_write");
00110   bool hydjet  = cmdline.present("-hydjet");
00111   double ktR   = cmdline.double_val("-r",1.0);
00112   double inclkt = cmdline.double_val("-incl",-1.0);
00113   int    excln  = cmdline.int_val   ("-excln",-1);
00114   double excld  = cmdline.double_val("-excld",-1.0);
00115   double etamax = cmdline.double_val("-etamax",1.0e305);
00116   bool   show_constituents = cmdline.present("-const");
00117   bool   massless = cmdline.present("-massless");
00118   int  nev     = cmdline.int_val("-nev",1);
00119   bool add_dense_coverage = cmdline.present("-dense");
00120 
00121   // The following option causes the Cambridge algo to be used.
00122   // Note that currently the only output that works sensibly here is
00123   // "-incl 0"
00124   fj::JetAlgorithm jet_algorithm;
00125   if (cmdline.present("-cam")) {
00126     jet_algorithm = fj::cambridge_algorithm;
00127   } else {
00128     jet_algorithm = fj::kt_algorithm;
00129   }
00130 
00131   if (!cmdline.all_options_used()) {cerr << 
00132       "Error: some options were not recognized"<<endl; 
00133     exit(-1);}
00134 
00135 
00136   for (int iev = 0; iev < nev; iev++) {
00137   vector<fj::PseudoJet> jets;
00138   string line;
00139   int  ndone = 0;
00140   while (getline(cin, line)) {
00141       //cout << line<<endl;
00142     istringstream linestream(line);
00143     if (line == "#END") {
00144       ndone += 1;
00145       if (ndone == combine) {break;}
00146     }
00147     if (line.substr(0,1) == "#") {continue;}
00148     valarray<double> fourvec(4);
00149     if (hydjet) {
00150       // special reading from hydjet.txt event record (though actually
00151       // this is supposed to be a standard pythia event record, so
00152       // being able to read from it is perhaps not so bad an idea...)
00153       int ii, istat,id,m1,m2,d1,d2;
00154       double mass;
00155       linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00156                  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00157       // current file contains mass of particle as 4th entry
00158       if (istat == 1) {
00159         fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00160                           +pow2(fourvec[2])+pow2(mass));
00161       }
00162     } else {
00163       if (massless) {
00164         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00165         fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00166       else {
00167         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00168       }
00169     }
00170     fj::PseudoJet psjet(fourvec);
00171     if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00172   }
00173 
00174   // add a fake underlying event which is very soft, uniformly distributed
00175   // in eta,phi so as to allow one to reconstruct the area that is associated
00176   // with each jet.
00177   if (add_dense_coverage) {
00178     srand(2);
00179     int nphi = 60;
00180     int neta = 100;
00181     double kt = 1e-1;
00182     for (int iphi = 0; iphi<nphi; iphi++) {
00183       for (int ieta = -neta; ieta<neta+1; ieta++) {
00184         double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00185         double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00186         kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00187         double pminus = kt*exp(-eta);
00188         double pplus  = kt*exp(+eta);
00189         double px = kt*sin(phi);
00190         double py = kt*cos(phi);
00191         //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00192         fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00193         jets.push_back(mom);
00194       }
00195     }
00196   }
00197   
00198   fj::JetDefinition jet_def(jet_algorithm, ktR, strategy);
00199 
00200   for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00201     fj::ClusterSequence clust_seq(jets,jet_def,write);
00202     if (irepeat != 0) {continue;}
00203     cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00204     cout << "strategy used =  "<< clust_seq.strategy_string()<< endl;
00205 
00206     // now provide some nice output...
00207     if (inclkt >= 0.0) {
00208       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00209       for (size_t j = 0; j < jets.size(); j++) {
00210         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00211         if (show_constituents) {
00212           vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00213           for (size_t k = 0; k < const_jets.size(); k++) {
00214             printf("        jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00215                    const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00216           }
00217           cout << "\n\n";
00218         }
00219       }
00220     }
00221 
00222     if (excln > 0) {
00223       vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00224  
00225       cout << "Printing "<<excln<<" exclusive jets\n";
00226       for (size_t j = 0; j < jets.size(); j++) {
00227         printf("%5u %15.8f %15.8f %15.8f\n",
00228                //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00229                j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00230       }
00231     }
00232 
00233     if (excld > 0.0) {
00234       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00235       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00236       for (size_t j = 0; j < jets.size(); j++) {
00237         printf("%5u %15.8f %15.8f %15.8f\n",
00238                j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00239       }
00240     }
00241     
00242     // useful for testing that recombination sequences are unique
00243     if (unique_write) {
00244       vector<int> unique_history = clust_seq.unique_history_order();
00245       // construct the inverse of the above mapping
00246       vector<int> inv_unique_history(clust_seq.history().size());
00247       for (unsigned int i = 0; i < unique_history.size(); i++) {
00248         inv_unique_history[unique_history[i]] = i;}
00249 
00250       for (unsigned int i = 0; i < unique_history.size(); i++) {
00251         fj::ClusterSequence::history_element el = 
00252           clust_seq.history()[unique_history[i]];
00253         int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00254         int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00255         printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00256       }
00257     }
00258 
00259   } // irepeat
00260 
00261   } // iev
00262 }

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