Skip to content
Snippets Groups Projects
Commit bcf2c240 authored by Dmitry Kalinkin's avatar Dmitry Kalinkin
Browse files

rm benchmarks/dis/generator

parent b76f32aa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
benchmarks/dis/generator/gen_central_electrons.cxx deleted 100644 → 0
+ 0
− 83
View file @ b76f32aa
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include "HepMC3/Print.h"
#include <iostream>
#include<random>
#include<cmath>
#include <math.h>
#include <TMath.h>
using namespace HepMC3;
/** Generate electrons in the central region.
* This is for testing detectors in the "barrel" region.
*/
void gen_central_electrons(int n_events = 100,
const char* out_fname = "central_electrons.hepmc")
{
double cos_theta_min = std::cos( 50.0*(M_PI/180.0));
double cos_theta_max = std::cos(130.0*(M_PI/180.0));
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Random number generator
TRandom *r1 = new TRandom();
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
// pdgid 11 - electron
// pdgid 111 - pi0
// pdgid 2212 - proton
GenParticlePtr p1 =
std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
GenParticlePtr p2 = std::make_shared<GenParticle>(
FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
Double_t p = r1->Uniform(1.0, 10.0);
Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
Double_t th = std::acos(costh);
Double_t px = p * std::cos(phi) * std::sin(th);
Double_t py = p * std::sin(phi) * std::sin(th);
Double_t pz = p * std::cos(th);
// Generates random vectors, uniformly distributed over the surface of a
// sphere of given radius, in this case momentum.
// r1->Sphere(px, py, pz, p);
//std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
// type 1 is final state
// pdgid 11 - electron 0.510 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(
FourVector(
px, py, pz,
sqrt(p*p + (0.000511 * 0.000511))),
11, 1);
GenVertexPtr v1 = std::make_shared<GenVertex>();
v1->add_particle_in(p1);
v1->add_particle_in(p2);
v1->add_particle_out(p3);
evt.add_vertex(v1);
if (events_parsed == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
hepmc_output.write_event(evt);
if (events_parsed % 10000 == 0) {
std::cout << "Event: " << events_parsed << std::endl;
}
evt.clear();
}
hepmc_output.close();
std::cout << "Events parsed and written: " << events_parsed << std::endl;
}
benchmarks/dis/generator/pythia_dis.cxx deleted 100644 → 0
+ 0
− 263
View file @ b76f32aa
#include "Pythia8/Pythia.h"
#include "Pythia8Plugins/HepMC3.h"
#include <cassert>
#include <unistd.h>
using namespace Pythia8;
#include "clipp.h"
using namespace clipp;
using std::string;
//______________________________________________________________________________
enum class mode { none, help, list, part };
struct settings {
double E_electron = 18.0; // GeV
double E_ion = 275.0; // GeV
std::string outfile = "dis.hepmc";
int ion_PID = 2212;
int electron_PID = 11;
bool help = false;
bool success = false;
double Q2_min = 4.0;
int N_events = 1000;
mode selected = mode::none;
};
template <typename T>
void print_usage(T cli, const char* argv0)
{
// used default formatting
std::cout << "Usage:\n" << usage_lines(cli, argv0) << "\nOptions:\n" << documentation(cli) << '\n';
}
//______________________________________________________________________________
template <typename T>
void print_man_page(T cli, const char* argv0)
{
// all formatting options (with their default values)
auto fmt = clipp::doc_formatting{}
.start_column(8) // column where usage lines and documentation starts
.doc_column(20) // parameter docstring start col
.indent_size(4) // indent of documentation lines for children of a
// documented group
.line_spacing(0) // number of empty lines after single documentation lines
.paragraph_spacing(1) // number of empty lines before and after paragraphs
.flag_separator(", ") // between flags of the same parameter
.param_separator(" ") // between parameters
.group_separator(" ") // between groups (in usage)
.alternative_param_separator("|") // between alternative flags
.alternative_group_separator(" | ") // between alternative groups
.surround_group("(", ")") // surround groups with these
.surround_alternatives("(", ")") // surround group of alternatives with these
.surround_alternative_flags("", "") // surround alternative flags with these
.surround_joinable("(", ")") // surround group of joinable flags with these
.surround_optional("[", "]") // surround optional parameters with these
.surround_repeat("", "..."); // surround repeatable parameters with these
auto mp = make_man_page(cli, argv0, fmt);
mp.prepend_section("DESCRIPTION", "Simple pythia dis generator.");
mp.append_section("EXAMPLES", " $ pythia_dis [output_file]");
std::cout << mp << "\n";
}
//______________________________________________________________________________
settings cmdline_settings(int argc, char* argv[])
{
settings s;
auto lastOpt = " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
value("file", s.outfile).if_missing([] {
std::cout << "You need to provide an output filename argument!\n";
}) % "output file");
auto cli = (command("help").set(s.selected, mode::help) | lastOpt);
assert(cli.flags_are_prefix_free());
auto res = parse(argc, argv, cli);
// if( res.any_error() ) {
// s.success = false;
// std::cout << make_man_page(cli, argv[0]).prepend_section("error: ",
// " The best
// thing since
// sliced bread.");
// return s;
//}
s.success = true;
if (s.selected == mode::help) {
print_man_page<decltype(cli)>(cli, argv[0]);
};
return s;
}
//______________________________________________________________________________
int main(int argc, char* argv[])
{
settings s = cmdline_settings(argc, argv);
if (!s.success) {
return 1;
}
if (s.selected == mode::help) {
return 0;
}
// Beam energies, minimal Q2, number of events to generate.
double eProton = s.E_ion;
double eElectron = s.E_electron;
double Q2min = s.Q2_min;
int nEvent = s.N_events;
// Generator. Shorthand for event.
Pythia pythia;
Event& event = pythia.event;
// Set up incoming beams, for frame with unequal beam energies.
pythia.readString("Beams:frameType = 2");
// BeamA = proton.
pythia.readString("Beams:idA = " + std::to_string(s.ion_PID));
pythia.settings.parm("Beams:eA", eProton);
// BeamB = electron.
pythia.readString("Beams:idB = 11");
pythia.settings.parm("Beams:eB", eElectron);
// Set up DIS process within some phase space.
// Neutral current (with gamma/Z interference).
pythia.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
// Uncomment to allow charged current.
// pythia.readString("WeakBosonExchange:ff2ff(t:W) = on");
// Phase-space cut: minimal Q2 of process.
pythia.settings.parm("PhaseSpace:Q2Min", Q2min);
// Set dipole recoil on. Necessary for DIS + shower.
pythia.readString("SpaceShower:dipoleRecoil = on");
// Allow emissions up to the kinematical limit,
// since rate known to match well to matrix elements everywhere.
pythia.readString("SpaceShower:pTmaxMatch = 2");
// QED radiation off lepton not handled yet by the new procedure.
pythia.readString("PDF:lepton = off");
pythia.readString("TimeShower:QEDshowerByL = off");
// Initialize.
pythia.init();
// Interface for conversion from Pythia8::Event to HepMC one.
HepMC3::Pythia8ToHepMC3 toHepMC;
// Specify file where HepMC events will be stored.
HepMC3::WriterAscii ascii_io(argv[1]);
cout << endl << endl << endl;
// Histograms.
double Wmax = sqrt(4. * eProton * eElectron);
Hist Qhist("Q [GeV]", 100, 0., 50.);
Hist Whist("W [GeV]", 100, 0., Wmax);
Hist xhist("x", 100, 0., 1.);
Hist yhist("y", 100, 0., 1.);
Hist pTehist("pT of scattered electron [GeV]", 100, 0., 50.);
Hist pTrhist("pT of radiated parton [GeV]", 100, 0., 50.);
Hist pTdhist("ratio pT_parton/pT_electron", 100, 0., 5.);
double sigmaTotal(0.), errorTotal(0.);
bool wroteRunInfo = false;
// Get the inclusive x-section by summing over all process x-sections.
double xs = 0.;
for (int i = 0; i < pythia.info.nProcessesLHEF(); ++i) {
xs += pythia.info.sigmaLHEF(i);
}
// Begin event loop.
for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
if (!pythia.next())
continue;
double sigmaSample = 0., errorSample = 0.;
// Four-momenta of proton, electron, virtual photon/Z^0/W^+-.
Vec4 pProton = event[1].p();
Vec4 peIn = event[4].p();
Vec4 peOut = event[6].p();
Vec4 pPhoton = peIn - peOut;
// Q2, W2, Bjorken x, y.
double Q2 = -pPhoton.m2Calc();
double W2 = (pProton + pPhoton).m2Calc();
double x = Q2 / (2. * pProton * pPhoton);
double y = (pProton * pPhoton) / (pProton * peIn);
// Fill kinematics histograms.
Qhist.fill(sqrt(Q2));
Whist.fill(sqrt(W2));
xhist.fill(x);
yhist.fill(y);
pTehist.fill(event[6].pT());
// pT spectrum of partons being radiated in shower.
for (int i = 0; i < event.size(); ++i)
if (event[i].statusAbs() == 43) {
pTrhist.fill(event[i].pT());
pTdhist.fill(event[i].pT() / event[6].pT());
}
// Get event weight(s).
double evtweight = pythia.info.weight();
// Do not print zero-weight events.
if (evtweight == 0.)
continue;
// Create a GenRunInfo object with the necessary weight names and write
// them to the HepMC3 file only once.
if (!wroteRunInfo) {
shared_ptr<HepMC3::GenRunInfo> genRunInfo;
genRunInfo = make_shared<HepMC3::GenRunInfo>();
vector<string> weight_names = pythia.info.weightNameVector();
genRunInfo->set_weight_names(weight_names);
ascii_io.set_run_info(genRunInfo);
ascii_io.write_run_info();
wroteRunInfo = true;
}
// Construct new empty HepMC event.
HepMC3::GenEvent hepmcevt;
// Work with weighted (LHA strategy=-4) events.
double normhepmc = 1.;
if (abs(pythia.info.lhaStrategy()) == 4)
normhepmc = 1. / double(1e9 * nEvent);
// Work with unweighted events.
else
normhepmc = xs / double(1e9 * nEvent);
// Set event weight
// hepmcevt.weights().push_back(evtweight*normhepmc);
// Fill HepMC event
toHepMC.fill_next_event(pythia, &hepmcevt);
// Add the weight of the current event to the cross section.
sigmaTotal += evtweight * normhepmc;
sigmaSample += evtweight * normhepmc;
errorTotal += pow2(evtweight * normhepmc);
errorSample += pow2(evtweight * normhepmc);
// Report cross section to hepmc
shared_ptr<HepMC3::GenCrossSection> xsec;
xsec = make_shared<HepMC3::GenCrossSection>();
// First add object to event, then set cross section. This order ensures
// that the lengths of the cross section and the weight vector agree.
hepmcevt.set_cross_section(xsec);
xsec->set_cross_section(sigmaTotal * 1e9, pythia.info.sigmaErr() * 1e9);
// Write the HepMC event to file. Done with it.
ascii_io.write_event(hepmcevt);
// End of event loop. Statistics and histograms.
}
pythia.stat();
cout << Qhist << Whist << xhist << yhist << pTehist << pTrhist << pTdhist;
return 0;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment