Main Page | Namespace List | Class Hierarchy | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 
00079     // Print out the number of final state jets
00080     //std::cout << "Number of final state jets: " << ev.getNJets() << std::endl;
00081     if (write) {
00083       std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00084       
00086       std::vector<KtJet::KtLorentzVector>::const_iterator itr = jets.begin();
00087       for( ; itr != jets.end() ; ++itr) {
00088         std::cout << "Jets Pt2: " << pow2((*itr).perp()) << std::endl; 
00089       }
00090     }
00091 
00092     if (inclkt >= 0.0) {
00093       // Retrieve the final state jets from KtEvent sorted by Pt
00094       std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00095     
00096       // Print out index, rap, phi, pt 
00097       for (size_t j = 0; j < jets.size(); j++) {
00098         if (jets[j].perp() < inclkt) {break;}
00099         double phi = jets[j].phi();
00100         if (phi < 0.0) {phi += fastjet::twopi;}
00101         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rapidity(),phi,jets[j].perp());
00102       }
00103     }
00104 
00105     if (excln > 0) {
00106       ev.findJetsN(excln);
00107       vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00108       cout << "Printing "<<excln<<" exclusive jets\n";
00109       for (size_t j = 0; j < jets.size(); j++) {
00110         printf("%5u %15.8f %15.8f %15.8f\n",j,
00111                jets[j].rapidity(),jets[j].phi(),jets[j].perp());
00112       }
00113     }
00114 
00115     if (excld > 0.0) {
00116       ev.findJetsD(excld);
00117       vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00118       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00119       for (size_t j = 0; j < jets.size(); j++) {
00120         printf("%5u %15.8f %15.8f %15.8f\n",j,
00121                jets[j].rapidity(),jets[j].phi(),jets[j].perp());
00122       }
00123     }
00124 
00125 
00126   }
00127 
00128 }
00129 }

Generated on Mon Apr 2 20:57:48 2007 for fastjet by  doxygen 1.4.2