FastJet 3.0.6
ktjet_timing.cc
00001 //----------------------------------------------------------------------
00002 // ktjet_timing.cc: Program similar to fastjet_timing.cc, to be used
00003 //                  for timing and result comparisons. For usage, see
00004 //                  explanations given at the start of fastjet_timing.cc
00005 // 
00006 // Note: - some options present in fastjet_timing are missing
00007 //       - one or two behave differently, notably -write
00008 //
00009 #include<iostream>
00010 #include<sstream>
00011 #include<valarray>
00012 #include<vector>
00013 #include<cstddef> // for size_t
00014 #include "CmdLine.hh"
00015 #include "fastjet/internal/numconsts.hh"
00016 
00017 
00018 /** Need to include these KtJet Headers */
00019 #include "KtJet/KtEvent.h"
00020 #include "KtJet/KtLorentzVector.h"
00021 using namespace std;
00022 using namespace KtJet;
00023 
00024 inline double pow2(const double x) {return x*x;}
00025 
00026 /// a program to test and time the kt algorithm as implemented in ktjet
00027 int main (int argc, char ** argv) {
00028 
00029   CmdLine cmdline(argc,argv);
00030   //bool clever = !cmdline.present("-dumb");
00031   int  repeat = cmdline.int_val("-repeat",1);
00032   int  combine = cmdline.int_val("-combine",1);
00033   bool write   = cmdline.present("-write");
00034   double ktR   = cmdline.double_val("-r",1.0);
00035   double inclkt = cmdline.double_val("-incl",-1.0);
00036   int    excln  = cmdline.int_val   ("-excln",-1);
00037   double excld  = cmdline.double_val("-excld",-1.0);
00038   int  nev     = cmdline.int_val("-nev",1);
00039   bool   massless = cmdline.present("-massless");
00040   bool   get_all_dij   = cmdline.present("-get-all-dij");
00041 
00042 
00043   for (int iev = 0; iev < nev; iev++) {
00044   vector<KtJet::KtLorentzVector> jets;
00045   string line;
00046   int  ndone = 0;
00047   while (getline(cin, line)) {
00048       //cout << line<<endl;
00049     istringstream linestream(line);
00050     if (line == "#END") {
00051       ndone += 1;
00052       if (ndone == combine) {break;}
00053     }
00054     if (line.substr(0,1) == "#") {continue;}
00055     valarray<double> fourvec(4);
00056     linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00057     if (massless) {
00058       linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00059       fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00060     else {
00061       linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00062     }
00063     KtJet::KtLorentzVector p(fourvec[0],fourvec[1],fourvec[2],fourvec[3]);
00064     jets.push_back(p);
00065   }
00066   
00067   // set KtEvent flags
00068   int type, angle, recom;
00069   ostringstream info;
00070   if (cmdline.present("-eekt")) {
00071     type  = 1; // e+e-
00072     angle = 1; // angular
00073     recom = 1; // E
00074     info << "Algorithm: KtJet e+e- kt algorithm" ;
00075   } else {
00076     type  = 4; // PP
00077     angle = 2; // delta R
00078     recom = 1; // E
00079     info << "Algorithm: KtJet (long.inv.) with R = " << ktR ;
00080   }
00081   //double rparameter = 1.0;
00082 
00083   for (int i = 0; i < repeat ; i++) {
00084     // Construct the KtEvent object 
00085     KtJet::KtEvent ev(jets,type,angle,recom,ktR);
00086 
00087     if (i!=0) {continue;}
00088     int nparticles = jets.size();
00089     cout << "Number of particles = "<< nparticles << endl;
00090     cout << info.str() << endl;
00091 
00092     // Print out the number of final state jets
00093     //std::cout << "Number of final state jets: " << ev.getNJets() << std::endl;
00094     if (write) {
00095       /** Retrieve the final state jets from KtEvent sorted by Pt*/
00096       std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00097       
00098       /** Print out jets 4-momentum and Pt */
00099       std::vector<KtJet::KtLorentzVector>::const_iterator itr = jets.begin();
00100       for( ; itr != jets.end() ; ++itr) {
00101         std::cout << "Jets Pt2: " << pow2((*itr).perp()) << std::endl; 
00102       }
00103     }
00104 
00105     if (inclkt >= 0.0) {
00106       // Retrieve the final state jets from KtEvent sorted by Pt
00107       std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00108     
00109       // Print out index, rap, phi, pt 
00110       for (size_t j = 0; j < jets.size(); j++) {
00111         if (jets[j].perp() < inclkt) {break;}
00112         double phi = jets[j].phi();
00113         if (phi < 0.0) {phi += fastjet::twopi;}
00114         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rapidity(),phi,jets[j].perp());
00115       }
00116     }
00117 
00118     if (excln > 0) {
00119       ev.findJetsN(excln);
00120       vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00121       cout << "Printing "<<excln<<" exclusive jets\n";
00122       for (size_t j = 0; j < jets.size(); j++) {
00123         double phi = jets[j].phi();
00124         if (phi < 0) phi += fastjet::twopi;
00125         printf("%5u %15.8f %15.8f %15.8f\n",j,
00126                jets[j].rapidity(),phi,jets[j].perp());
00127       }
00128     }
00129 
00130     if (excld > 0.0) {
00131       ev.findJetsD(excld);
00132       vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00133       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00134       for (size_t j = 0; j < jets.size(); j++) {
00135         double phi = jets[j].phi();
00136         if (phi < 0) phi += fastjet::twopi;
00137         printf("%5u %15.8f %15.8f %15.8f\n",j,
00138                jets[j].rapidity(),phi,jets[j].perp());
00139       }
00140     }
00141 
00142     if (get_all_dij) {
00143       for (int i = nparticles-1; i > 0; i--) {
00144         printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, ev.getDMerge(i));
00145       }
00146     }
00147 
00148   }
00149 
00150 }
00151 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends