Page MenuHomeHEPForge
Files F10275477

No OneTemporary
Actions

View Options

diff --git a/main.cc b/main.cc
index 39ed5f7..6925cf1 100644
--- a/main.cc
+++ b/main.cc
@@ -1,125 +1,126 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <algorithm>
#include "TFile.h"
#include "LHEF.h"
#include "src/Analysis.hpp"
#include "Common.h"
#include "src/Event.hpp"
// Common definitions. will be set somewhere else in the end.
// These definitions are passed to the clases when the constructors are
// called
// Define Jet Measure
-static double Rparam(0.4);
+static double Rparam(0.6);
static fastjet::Strategy strategy = fastjet::Best;
static fastjet::RecombinationScheme recomb_scheme = fastjet::E_scheme;
// Define jet for fastjet
static fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, Rparam, recomb_scheme, strategy);
-static double jetptmin = 20.;
+static double jetptmin = 30.;
int main(int argn, char** argv) {
// Check the user input. LHE reader will check later if the file
// actually exists.
if (argn != 2) {
std::cerr << "\n# Usage:\n./HEJ input_file\n\n";
return EXIT_FAILURE;
}
std::string infile = argv[1];
// Time Profiling Init.
clock_t StartTime, EndTime;
StartTime = clock();
// Initialize histogrmas.
RHEJ::histograms hist;
// Initialize process type tester and pdf classes.
LHEF::Reader* reader = new LHEF::Reader(infile);
double E_beam = reader->heprup.EBMUP.first;
if(reader->heprup.EBMUP.second != E_beam){
std::cerr << "asymmetric beam: " << E_beam << " ---> <--- "
<< reader->heprup.EBMUP.second << '\n';
return EXIT_FAILURE;
}
RHEJ::Event rhej{E_beam, jet_def, jetptmin};
+ rhej.pdf("NNPDF23_nlo_as_0119", 0, 1, 1);
// keep track of the event number
int nevent(0),nfkl(0),nfkl2(0);
double sumOfwt(0.),wterror(0.);
// Loop over the events in the inputfile
while(reader->readEvent()){
// if (nevent > 50000)
// break;
nevent++;
if (nevent % 10000 == 0)
std::cout << "event number " << nevent << std::endl;
// Get Event Weight
double weight = reader->hepeup.weight();
// calculate HEJ reweighting factor
double reweight = rhej.reweight(reader->hepeup);
switch(rhej.error()){
case RHEJ::error::not_FKL:
continue;
case RHEJ::error::bad_jets:
nfkl++;
continue;
default:
// count fkl events which pass jet cuts
nfkl++;
nfkl2++;
}
// Add into the sum of weights.
double wt_tmp = reweight*weight;
sumOfwt+= wt_tmp;
wterror+= wt_tmp*wt_tmp;
hist.fill(rhej);
if (sumOfwt<0 || sumOfwt>=10E+30){
std::cerr<<"gwt "<<rhej.weights().gwt<<" s2jwt "<<rhej.weights().s2jwt<<" splitwt "<<rhej.weights().splitwt<<" qwt "<<rhej.weights().qwt<<std::endl;
std::cerr<<"trials "<<rhej.trials()<<std::endl;
std::cerr<<"reweight "<<reweight<<" weight "<<weight<<std::endl;
std::cerr << "Line "<< __LINE__ << " in file "<<__FILE__<<"\n";
}
} // main event loop
TFile out_file("RHEJ.root", "RECREATE");
hist.write();
out_file.Close();
// std::cout<<" pdf calls wroks "<<pdf<<" jet clustering works! " <<jets.size()<<std::endl;
std::cout<<"Sum of weights = "<<sumOfwt<< " +/- "<<sqrt(wterror)<<std::endl;
std::cout<<"Events processed: "<<nevent<<" Number of FKL events: "<<nfkl<< " non-FKL "<<nevent-nfkl<<" passed FKL events "<<nfkl2<<std::endl;
EndTime = clock();
std::cout << "\nTask Runtime: " << static_cast <double> (EndTime - StartTime) / CLOCKS_PER_SEC << " seconds.\n" << std::endl;
// Return success
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Fri, Apr 4, 9:37 PM (15 h, 53 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4737419
Default Alt Text
(3 KB)

Event Timeline