fastjet_timing_plugins.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: fastjet_timing.cc 293 2006-08-17 19:38:38Z salam $
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 //----------------------------------------------------------------------
00089 #include "fastjet/PseudoJet.hh"
00090 #include "fastjet/ClusterSequence.hh"
00091 #include<iostream>
00092 #include<sstream>
00093 #include<fstream>
00094 #include<valarray>
00095 #include<vector>
00096 #include <cstdlib>
00097 #include<cstddef> // for size_t
00098 #include "CmdLine.hh"
00099 
00100 // get info on how fastjet was configured
00101 #include "fastjet/config.h"
00102 
00103 #ifdef ENABLE_PLUGIN_SISCONE
00104 #include "fastjet/SISConePlugin.hh"
00105 #endif
00106 #ifdef ENABLE_PLUGIN_CDFCONES
00107 #include "fastjet/CDFMidPointPlugin.hh"
00108 #include "fastjet/CDFJetCluPlugin.hh"
00109 #endif
00110 #ifdef ENABLE_PLUGIN_PXCONE
00111 #include "fastjet/PxConePlugin.hh"
00112 #endif
00113 
00114 using namespace std;
00115 
00116 // to avoid excessive typing, define an abbreviation for the 
00117 // fastjet namespace
00118 namespace fj = fastjet;
00119 
00120 inline double pow2(const double x) {return x*x;}
00121 
00123 int main (int argc, char ** argv) {
00124 
00125   CmdLine cmdline(argc,argv);
00126   // allow the use to specify the fj::Strategy either through the
00127   // -clever or the -strategy options (both will take numerical
00128   // values); the latter will override the former.
00129   fj::Strategy  strategy  = fj::Strategy(cmdline.int_val("-strategy",
00130                                         cmdline.int_val("-clever", fj::Best)));
00131   int  repeat  = cmdline.int_val("-repeat",1);
00132   int  combine = cmdline.int_val("-combine",1);
00133   bool write   = cmdline.present("-write");
00134   bool unique_write = cmdline.present("-unique_write");
00135   bool hydjet  = cmdline.present("-hydjet");
00136   double ktR   = cmdline.double_val("-r",1.0);
00137   ktR   = cmdline.double_val("-R",ktR); // allow -r and -R
00138   double inclkt = cmdline.double_val("-incl",-1.0);
00139   int    excln  = cmdline.int_val   ("-excln",-1);
00140   double excld  = cmdline.double_val("-excld",-1.0);
00141   double etamax = cmdline.double_val("-etamax",1.0e305);
00142   bool   show_constituents = cmdline.present("-const");
00143   bool   massless = cmdline.present("-massless");
00144   int  nev     = cmdline.int_val("-nev",1);
00145   bool add_dense_coverage = cmdline.present("-dense");
00146 
00147   bool show_cones = cmdline.present("-cones"); // only works for siscone
00148 
00149   // for cone algorithms
00150   // allow -f and -overlap
00151   double overlap_threshold = cmdline.double_val("-overlap",0.5);
00152   overlap_threshold = cmdline.double_val("-f",overlap_threshold); 
00153   double seed_threshold = cmdline.double_val("-seed",1.0);
00154 
00155   // for printing jets to a file for reading by root
00156   string rootfile = cmdline.value<string>("-root","");
00157 
00158   // The following option causes the Cambridge algo to be used.
00159   // Note that currently the only output that works sensibly here is
00160   // "-incl 0"
00161   fj::JetDefinition jet_def;
00162   if (cmdline.present("-cam")) {
00163     jet_def = fj::JetDefinition(fj::cambridge_algorithm, ktR, strategy);
00164   } else if (cmdline.present("-antikt")) {
00165     jet_def = fj::JetDefinition(fj::antikt_algorithm, ktR, strategy);
00166   } else if (cmdline.present("-midpoint")) {
00167 #ifdef ENABLE_PLUGIN_CDFCONES
00168     typedef fj::CDFMidPointPlugin MPPlug; // for brevity
00169     double cone_area_fraction = 1.0;
00170     int    max_pair_size = 2;
00171     int    max_iterations = 100;
00172     MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
00173     if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
00174     if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
00175     if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
00176     if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00177     jet_def = fj::JetDefinition( new fj::CDFMidPointPlugin (
00178                                       seed_threshold, ktR, 
00179                                       cone_area_fraction, max_pair_size,
00180                                       max_iterations, overlap_threshold,
00181                                       sm_scale));
00182 #else  // ENABLE_PLUGIN_CDFCONES
00183     cerr << "midpoint requested, but not available for this compilation" << endl;
00184 #endif // ENABLE_PLUGIN_CDFCONES
00185   } else if (cmdline.present("-pxcone")) {
00186 #ifdef ENABLE_PLUGIN_PXCONE
00187     double min_jet_energy = 5.0;
00188     jet_def = fj::JetDefinition( new fj::PxConePlugin (
00189                                       ktR, min_jet_energy,
00190                                       overlap_threshold));
00191 #else  // ENABLE_PLUGIN_PXCONE
00192     cerr << "pxcone requested, but not available for this compilation" << endl;
00193     exit(-1);
00194 #endif // ENABLE_PLUGIN_PXCONE
00195   } else if (cmdline.present("-jetclu")) {
00196 #ifdef ENABLE_PLUGIN_CDFCONES
00197     jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin (
00198                                       ktR, overlap_threshold, seed_threshold));
00199 #else  // ENABLE_PLUGIN_CDFCONES
00200     cerr << "jetclu requested, but not available for this compilation" << endl;
00201 #endif // ENABLE_PLUGIN_CDFCONES
00202   } else if (cmdline.present("-siscone")) {
00203 #ifdef ENABLE_PLUGIN_SISCONE
00204     typedef fj::SISConePlugin SISPlug; // for brevity
00205     int npass = cmdline.value("-npass",0);
00206     double sisptmin = cmdline.value("-sisptmin",0.0);
00207     SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
00208     if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
00209     if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
00210     if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
00211     if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
00212     jet_def = fj::JetDefinition(plugin);
00213 #else  // ENABLE_PLUGIN_SISCONE
00214     cerr << "jetclu requested, but not available for this compilation" << endl;
00215 #endif // ENABLE_PLUGIN_SISCONE
00216   } else {
00217     jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy);
00218   }
00219 
00220 
00221 
00222   if (!cmdline.all_options_used()) {cerr << 
00223       "Error: some options were not recognized"<<endl; 
00224     exit(-1);}
00225 
00226 
00227   for (int iev = 0; iev < nev; iev++) {
00228   vector<fj::PseudoJet> jets;
00229   string line;
00230   int  ndone = 0;
00231   while (getline(cin, line)) {
00232       //cout << line<<endl;
00233     istringstream linestream(line);
00234     if (line == "#END") {
00235       ndone += 1;
00236       if (ndone == combine) {break;}
00237     }
00238     if (line.substr(0,1) == "#") {continue;}
00239     valarray<double> fourvec(4);
00240     if (hydjet) {
00241       // special reading from hydjet.txt event record (though actually
00242       // this is supposed to be a standard pythia event record, so
00243       // being able to read from it is perhaps not so bad an idea...)
00244       int ii, istat,id,m1,m2,d1,d2;
00245       double mass;
00246       linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00247                  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00248       // current file contains mass of particle as 4th entry
00249       if (istat == 1) {
00250         fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00251                           +pow2(fourvec[2])+pow2(mass));
00252       }
00253     } else {
00254       if (massless) {
00255         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00256         fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00257       else {
00258         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00259       }
00260     }
00261     fj::PseudoJet psjet(fourvec);
00262     if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00263   }
00264 
00265   // add a fake underlying event which is very soft, uniformly distributed
00266   // in eta,phi so as to allow one to reconstruct the area that is associated
00267   // with each jet.
00268   if (add_dense_coverage) {
00269     srand(2);
00270     int nphi = 60;
00271     int neta = 100;
00272     double kt = 1e-1;
00273     for (int iphi = 0; iphi<nphi; iphi++) {
00274       for (int ieta = -neta; ieta<neta+1; ieta++) {
00275         double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00276         double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00277         kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00278         double pminus = kt*exp(-eta);
00279         double pplus  = kt*exp(+eta);
00280         double px = kt*sin(phi);
00281         double py = kt*cos(phi);
00282         //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00283         fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00284         jets.push_back(mom);
00285       }
00286     }
00287   }
00288   
00289   for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00290     fj::ClusterSequence clust_seq(jets,jet_def,write);
00291     if (irepeat != 0) {continue;}
00292     cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00293     cout << "strategy used =  "<< clust_seq.strategy_string()<< endl;
00294     cout << "Algorithm: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl;
00295 
00296     // now provide some nice output...
00297     if (inclkt >= 0.0) {
00298       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00299       for (size_t j = 0; j < jets.size(); j++) {
00300         //printf("%5u %15.8f %15.8f %15.8e\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00301         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00302         if (show_constituents) {
00303           vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00304           for (size_t k = 0; k < const_jets.size(); k++) {
00305             printf("        jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00306                    const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00307           }
00308           cout << "\n\n";
00309         }
00310         // family tests
00311         //fj::PseudoJet parent1, parent2;
00312         //bool has_parents = clust_seq.has_parents(jets[j],parent1,parent2);
00313         //if (has_parents) {
00314         //  cout << "parent pt's: " << parent1.perp() << " " << parent2.perp() <<endl;
00315         //  fj::PseudoJet child_from1, child_from2;
00316         //  clust_seq.has_child(parent1, child_from1);
00317         //  clust_seq.has_child(parent2, child_from2);
00318         //  cout << "parents' children pt: " << child_from1.perp() << " " << child_from2.perp() << endl;
00319         //  fj::PseudoJet partner_from1, partner_from2;
00320         //  clust_seq.has_partner(parent1, partner_from1);
00321         //  clust_seq.has_partner(parent2, partner_from2);
00322         //  cout << "parents' partners' pt: " << partner_from1.perp() << " " << partner_from2.perp() << endl;
00323         //} else {
00324         //  cout << "has no parents" << endl;
00325         //}
00326       }
00327       if (rootfile != "") {
00328         ofstream ostr(rootfile.c_str());
00329         clust_seq.print_jets_for_root(jets,ostr);
00330       }
00331 
00332     }
00333 
00334     if (excln > 0) {
00335       vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00336  
00337       cout << "Printing "<<excln<<" exclusive jets\n";
00338       for (size_t j = 0; j < jets.size(); j++) {
00339         printf("%5u %15.8f %15.8f %15.8f\n",
00340                //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00341                j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00342       }
00343     }
00344 
00345     if (excld > 0.0) {
00346       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00347       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00348       for (size_t j = 0; j < jets.size(); j++) {
00349         printf("%5u %15.8f %15.8f %15.8f\n",
00350                j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00351       }
00352     }
00353 
00354     
00355     // useful for testing that recombination sequences are unique
00356     if (unique_write) {
00357       vector<int> unique_history = clust_seq.unique_history_order();
00358       // construct the inverse of the above mapping
00359       vector<int> inv_unique_history(clust_seq.history().size());
00360       for (unsigned int i = 0; i < unique_history.size(); i++) {
00361         inv_unique_history[unique_history[i]] = i;}
00362 
00363       for (unsigned int i = 0; i < unique_history.size(); i++) {
00364         fj::ClusterSequence::history_element el = 
00365           clust_seq.history()[unique_history[i]];
00366         int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00367         int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00368         printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00369       }
00370     }
00371 
00372 
00373     // provide some complementary information for SISCone 
00374     if (show_cones) {
00375       const fj::SISConeExtras * extras = 
00376         dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00377       cout << "most ambiguous split (difference in squared dist) = "
00378            << extras->most_ambiguous_split() << endl;
00379       vector<fastjet::PseudoJet> stable_cones(extras->stable_cones()); 
00380       stable_cones = sorted_by_rapidity(stable_cones);
00381       for (unsigned int i = 0; i < stable_cones.size(); i++) {
00382       //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
00383         printf("%5u %15.8f %15.8f %15.8f\n",
00384                i,stable_cones[i].rap(),stable_cones[i].phi(),
00385                stable_cones[i].perp() );
00386       //}
00387       }
00388     }
00389   } // irepeat
00390 
00391   } // iev
00392 }

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