ktjet_timing.cc

Go to the documentation of this file.
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 
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 
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 
00041   for (int iev = 0; iev < nev; iev++) {
00042   vector<KtJet::KtLorentzVector> jets;
00043   string line;
00044   int  ndone = 0;
00045   while (getline(cin, line)) {
00046       //cout << line<<endl;
00047     istringstream linestream(line);
00048     if (line == "#END") {
00049       ndone += 1;
00050       if (ndone == combine) {break;}
00051     }
00052     if (line.substr(0,1) == "#") {continue;}
00053     valarray<double> fourvec(4);
00054     linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00055     if (massless) {
00056       linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00057       fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00058     else {
00059       linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00060     }
00061     KtJet::KtLorentzVector p(fourvec[0],fourvec[1],fourvec[2],fourvec[3]);
00062     jets.push_back(p);
00063   }
00064   
00065 
00066   // set KtEvent flags
00067   int type  = 4; // PP
00068   int angle = 2; // delta R
00069   int recom = 1; // E
00070   //double rparameter = 1.0;
00071 
00072   for (int i = 0; i < repeat ; i++) {
00073     // Construct the KtEvent object 
00074     KtJet::KtEvent ev(jets,type,angle,recom,ktR);
00075 
00076     if (i!=0) {continue;}
00077     cout << "Number of particles = "<< jets.size() << endl;
00078     cout << "Algorithm: KtJet (long.inv.) with R = " << ktR << endl;
00079 
00080     // Print out the number of final state jets
00081     //std::cout << "Number of final state jets: " << ev.getNJets() << std::endl;
00082     if (write) {
00084       std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00085       
00087       std::vector<KtJet::KtLorentzVector>::const_iterator itr = jets.begin();
00088       for( ; itr != jets.end() ; ++itr) {
00089         std::cout << "Jets Pt2: " << pow2((*itr).perp()) << std::endl; 
00090       }
00091     }
00092 
00093     if (inclkt >= 0.0) {
00094       // Retrieve the final state jets from KtEvent sorted by Pt
00095       std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00096     
00097       // Print out index, rap, phi, pt 
00098       for (size_t j = 0; j < jets.size(); j++) {
00099         if (jets[j].perp() < inclkt) {break;}
00100         double phi = jets[j].phi();
00101         if (phi < 0.0) {phi += fastjet::twopi;}
00102         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rapidity(),phi,jets[j].perp());
00103       }
00104     }
00105 
00106     if (excln > 0) {
00107       ev.findJetsN(excln);
00108       vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00109       cout << "Printing "<<excln<<" exclusive jets\n";
00110       for (size_t j = 0; j < jets.size(); j++) {
00111         printf("%5u %15.8f %15.8f %15.8f\n",j,
00112                jets[j].rapidity(),jets[j].phi(),jets[j].perp());
00113       }
00114     }
00115 
00116     if (excld > 0.0) {
00117       ev.findJetsD(excld);
00118       vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00119       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00120       for (size_t j = 0; j < jets.size(); j++) {
00121         printf("%5u %15.8f %15.8f %15.8f\n",j,
00122                jets[j].rapidity(),jets[j].phi(),jets[j].perp());
00123       }
00124     }
00125 
00126 
00127   }
00128 
00129 }
00130 }

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