Skip to content
Snippets Groups Projects
Commit 58fb453c authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Merge branch 'dis_pythia' of eicweb.phy.anl.gov:EIC/benchmarks/physics_benchmarks into dis_pythia

parents 4ad052e0 5c96bf0c
No related branches found
No related tags found
No related merge requests found
---
BasedOnStyle: LLVM
BreakConstructorInitializersBeforeComma: true
ConstructorInitializerAllOnOneLineOrOnePerLine: true
Language: Cpp
BasedOnStyle: Chromium
AccessModifierOffset: -2
AlignAfterOpenBracket: Align
AlignConsecutiveAssignments: true
AlignConsecutiveDeclarations: true
AlignEscapedNewlines: Right
AlignOperands: true
AlignTrailingComments: true
AllowAllParametersOfDeclarationOnNextLine: true
AllowShortBlocksOnASingleLine: false
AllowShortCaseLabelsOnASingleLine: false
AllowShortFunctionsOnASingleLine: All
AllowShortIfStatementsOnASingleLine: false
AllowShortLoopsOnASingleLine: false
AlwaysBreakAfterDefinitionReturnType: None
AlwaysBreakAfterReturnType: None
AlwaysBreakBeforeMultilineStrings: false
AlwaysBreakTemplateDeclarations: true
BinPackArguments: true
BinPackParameters: true
BraceWrapping:
AfterClass: false
AfterControlStatement: false
AfterEnum: false
AfterFunction: true
AfterNamespace: false
AfterObjCDeclaration: false
AfterStruct: false
AfterUnion: false
AfterExternBlock: false
BeforeCatch: false
BeforeElse: false
IndentBraces: false
SplitEmptyFunction: true
SplitEmptyRecord: true
SplitEmptyNamespace: true
BreakBeforeBinaryOperators: None
BreakBeforeBraces: Custom
BreakBeforeInheritanceComma: false
BreakBeforeTernaryOperators: true
BreakConstructorInitializersBeforeComma: false
BreakConstructorInitializers: BeforeColon
BreakAfterJavaFieldAnnotations: false
BreakStringLiterals: true
ColumnLimit: 120
CommentPragmas: '^ IWYU pragma:'
CompactNamespaces: false
ConstructorInitializerAllOnOneLineOrOnePerLine: false
ConstructorInitializerIndentWidth: 4
ContinuationIndentWidth: 4
Cpp11BracedListStyle: true
Standard: Cpp11
#SpaceBeforeParens: ControlStatements
SpaceAfterControlStatementKeyword: true
PointerBindsToType: true
DerivePointerAlignment: false
DisableFormat: false
ExperimentalAutoDetectBinPacking: false
FixNamespaceComments: true
ForEachMacros:
- foreach
- Q_FOREACH
- BOOST_FOREACH
IncludeBlocks: Preserve
IncludeCategories:
- Regex: '^"(llvm|llvm-c|clang|clang-c)/'
Priority: 2
- Regex: '^(<|"(gtest|gmock|isl|json)/)'
Priority: 3
- Regex: '.*'
Priority: 1
IncludeIsMainRegex: '(Test)?$'
IndentCaseLabels: false
IndentPPDirectives: None
IndentWidth: 2
IndentWrappedFunctionNames: false
JavaScriptQuotes: Leave
JavaScriptWrapImports: true
KeepEmptyLinesAtTheStartOfBlocks: true
MacroBlockBegin: ''
MacroBlockEnd: ''
MaxEmptyLinesToKeep: 1
NamespaceIndentation: All
ObjCBlockIndentWidth: 2
ObjCSpaceAfterProperty: false
ObjCSpaceBeforeProtocolList: true
PenaltyBreakAssignment: 2
PenaltyBreakBeforeFirstCallParameter: 19
PenaltyBreakComment: 300
PenaltyBreakFirstLessLess: 120
PenaltyBreakString: 1000
PenaltyExcessCharacter: 1000000
PenaltyReturnTypeOnItsOwnLine: 60
PointerAlignment: Left
ReflowComments: true
SortIncludes: true
SortUsingDeclarations: true
SpaceAfterCStyleCast: false
SpaceAfterTemplateKeyword: true
SpaceBeforeAssignmentOperators: true
SpaceBeforeParens: ControlStatements
#SpaceBeforeRangeBasedForLoopColon: true
SpaceInEmptyParentheses: false
SpacesBeforeTrailingComments: 1
SpacesInAngles: false
SpacesInContainerLiterals: true
SpacesInCStyleCastParentheses: false
SpacesInParentheses: false
SpacesInSquareBrackets: false
Standard: Cpp11
TabWidth: 8
UseTab: Never
...
......@@ -44,6 +44,6 @@ summary:
- ./util/collect_benchmarks.py
artifacts:
paths:
- results/summary.json
- results/*
reports:
junit: ["results/dvcs/report2.xml"]
dis:generate:
stage: initialize
image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
needs: []
timeout: 1 hours
artifacts:
......
......@@ -35,9 +35,9 @@ JUGGLER_FILE_NAME_TAG="dis"
## =============================================================================
## Step 1: Dummy event generator
## TODO better file name that encodes the actual configuration we're running
echo "Compiling dis/src/pythia_dis.cc ..."
g++ dis/src/pythia_dis.cc -o pythia_dis \
-I/usr/local/include \
-I/usr/local/include -Iinclude \
-O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC \
-L/usr/local/lib -Wl,-rpath,/usr/local/lib -lpythia8 -ldl \
-L/usr/local/lib -Wl,-rpath,/usr/local/lib -lHepMC3
......@@ -45,10 +45,12 @@ if [[ "$?" -ne "0" ]] ; then
echo "ERROR compiling pythia"
exit 1
fi
echo "done"
pwd
./pythia_dis
./pythia_dis dis.hepmc
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running script"
echo "ERROR running pythia"
exit 1
fi
......
<<<<<<< HEAD
// main45.cc is a part of the PYTHIA event generator.
// Copyright (C) 2020 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
......@@ -37,6 +38,122 @@ int main( int argc, char* argv[] ){
double eElectron = 10.0;
double Q2min = 5.;
int nEvent = 10000;
=======
#include "Pythia8/Pythia.h"
#include "Pythia8Plugins/HepMC3.h"
#include <unistd.h>
#include <cassert>
using namespace Pythia8;
#include "clipp.h"
using namespace clipp;
using std::string;
//______________________________________________________________________________
enum class mode { none, help, list, part };
struct settings {
double E_electron = 10.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;
>>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a
// Generator. Shorthand for event.
Pythia pythia;
......@@ -45,7 +162,11 @@ int main( int argc, char* argv[] ){
// Set up incoming beams, for frame with unequal beam energies.
pythia.readString("Beams:frameType = 2");
// BeamA = proton.
<<<<<<< HEAD
pythia.readString("Beams:idA = 2212");
=======
pythia.readString("Beams:idA = " +std::to_string(s.ion_PID));
>>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a
pythia.settings.parm("Beams:eA", eProton);
// BeamB = electron.
pythia.readString("Beams:idB = 11");
......@@ -55,7 +176,11 @@ int main( int argc, char* argv[] ){
// Neutral current (with gamma/Z interference).
pythia.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
// Uncomment to allow charged current.
<<<<<<< HEAD
//pythia.readString("WeakBosonExchange:ff2ff(t:W) = on");
=======
// pythia.readString("WeakBosonExchange:ff2ff(t:W) = on");
>>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a
// Phase-space cut: minimal Q2 of process.
pythia.settings.parm("PhaseSpace:Q2Min", Q2min);
......@@ -80,6 +205,7 @@ int main( int argc, char* argv[] ){
cout << endl << endl << endl;
// Histograms.
<<<<<<< HEAD
double Wmax = sqrt(4.* eProton * eElectron);
Hist Qhist("Q [GeV]", 100, 0., 50.);
Hist Whist("W [GeV]", 100, 0., Wmax);
......@@ -95,6 +221,23 @@ int main( int argc, char* argv[] ){
double xs = 0.;
for (int i=0; i < pythia.info.nProcessesLHEF(); ++i)
xs += pythia.info.sigmaLHEF(i);
=======
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);
>>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a
// Begin event loop.
for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
......@@ -176,8 +319,11 @@ int main( int argc, char* argv[] ){
// Write the HepMC event to file. Done with it.
ascii_io.write_event(hepmcevt);
<<<<<<< HEAD
=======
>>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a
// End of event loop. Statistics and histograms.
}
pythia.stat();
......@@ -185,5 +331,9 @@ int main( int argc, char* argv[] ){
return 0;
<<<<<<< HEAD
}
=======
}
>>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a
......@@ -116,6 +116,7 @@ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts,
}
return {parts[first], parts[second]};
}
// Calculate the invariant mass of a given pair of particles
inline double
get_im(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
......@@ -123,6 +124,38 @@ get_im(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
return (particle_pair.first + particle_pair.second).mass();
}
// Calculate the transverse momentum of a given pair of particles
inline double
get_pt(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
double px_pair = (particle_pair.first + particle_pair.second).px();
double py_pair = (particle_pair.first + particle_pair.second).py();
return sqrt(px_pair*px_pair + py_pair*py_pair);
}
// Calculate the azimuthal angle of a given pair of particles
inline double
get_phi(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
double px_pair = (particle_pair.first + particle_pair.second).px();
double py_pair = (particle_pair.first + particle_pair.second).py();
double phi_pair = std::atan2(py_pair,px_pair);
//if(py_pair <= 0.) phi_pair = - phi_pair;
return phi_pair;
}
// Calculate the rapidity of a given pair of particles
inline double
get_y(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
double px_pair = (particle_pair.first + particle_pair.second).px();
double py_pair = (particle_pair.first + particle_pair.second).py();
double pz_pair = (particle_pair.first + particle_pair.second).pz();
double mass_pair = (particle_pair.first + particle_pair.second).mass();
double energy_pair = sqrt(mass_pair*mass_pair + px_pair*px_pair + py_pair*py_pair + pz_pair*pz_pair);
return 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair));
}
} // namespace util
#endif
......@@ -86,48 +86,164 @@ int vm_mass(const std::string& config_name) {
.Define("decay_pair_rec", find_decay_pair, {"p_rec"})
.Define("decay_pair_sim", find_decay_pair, {"p_sim"})
.Define("mass_rec", util::get_im, {"decay_pair_rec"})
.Define("mass_sim", util::get_im, {"decay_pair_sim"});
.Define("mass_sim", util::get_im, {"decay_pair_sim"})
.Define("pt_rec", util::get_pt, {"decay_pair_rec"})
.Define("pt_sim", util::get_pt, {"decay_pair_sim"})
.Define("phi_rec" , util::get_phi, {"decay_pair_rec"})
.Define("phi_sim" , util::get_phi, {"decay_pair_sim"})
.Define("rapidity_rec" , util::get_y, {"decay_pair_rec"})
.Define("rapidity_sim" , util::get_y, {"decay_pair_sim"});
// Define output histograms
auto h_im_rec = d_im.Histo1D(
{"h_im_rec", ";m_{ll'} (GeV);#", 100, -1.1, vm_mass + 5}, "mass_rec");
{"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec");
auto h_im_sim = d_im.Histo1D(
{"h_im_sim", ";m_{ll'} (GeV);#", 100, -1.1, vm_mass + 5}, "mass_sim");
{"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
auto h_pt_rec = d_im.Histo1D(
{"h_pt_rec", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_rec");
auto h_pt_sim = d_im.Histo1D(
{"h_pt_sim", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_sim");
auto h_phi_rec = d_im.Histo1D(
{"h_phi_rec", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_rec");
auto h_phi_sim = d_im.Histo1D(
{"h_phi_sim", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_sim");
auto h_y_rec = d_im.Histo1D(
{"h_y_rec", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_rec");
auto h_y_sim = d_im.Histo1D(
{"h_y_sim", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_sim");
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
// factorize out the plotting code moving forward.
{
TCanvas c{"canvas", "canvas", 800, 800};
gPad->SetLogx(false);
gPad->SetLogy(false);
auto& h0 = *h_im_sim;
auto& h1 = *h_im_rec;
TCanvas c{"canvas", "canvas", 1200, 1200};
c.Divide(2, 2, 0.0001, 0.0001);
//pad 1 mass
c.cd(1);
//gPad->SetLogx(false);
//gPad->SetLogy(false);
auto& h11 = *h_im_sim;
auto& h12 = *h_im_rec;
// histogram style
h0.SetLineColor(plot::kMpBlue);
h0.SetLineWidth(2);
h1.SetLineColor(plot::kMpOrange);
h1.SetLineWidth(2);
h11.SetLineColor(plot::kMpBlue);
h11.SetLineWidth(2);
h12.SetLineColor(plot::kMpOrange);
h12.SetLineWidth(2);
// axes
h0.GetXaxis()->CenterTitle();
h0.GetYaxis()->CenterTitle();
h11.GetXaxis()->CenterTitle();
h11.GetYaxis()->CenterTitle();
// draw everything
h0.DrawClone("hist");
h1.DrawClone("hist same");
h11.DrawClone("hist");
h12.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Invariant mass");
TText* tptr;
auto t = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t->SetFillColorAlpha(kWhite, 0);
t->SetTextFont(43);
t->SetTextSize(25);
tptr = t->AddText("simulated");
tptr->SetTextColor(plot::kMpBlue);
tptr = t->AddText("reconstructed");
tptr->SetTextColor(plot::kMpOrange);
t->Draw();
TText* tptr1;
auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t1->SetFillColorAlpha(kWhite, 0);
t1->SetTextFont(43);
t1->SetTextSize(25);
tptr1 = t1->AddText("simulated");
tptr1->SetTextColor(plot::kMpBlue);
tptr1 = t1->AddText("reconstructed");
tptr1->SetTextColor(plot::kMpOrange);
t1->Draw();
//pad 2 pt
c.cd(2);
//gPad->SetLogx(false);
//gPad->SetLogy(false);
auto& h21 = *h_pt_sim;
auto& h22 = *h_pt_rec;
// histogram style
h21.SetLineColor(plot::kMpBlue);
h21.SetLineWidth(2);
h22.SetLineColor(plot::kMpOrange);
h22.SetLineWidth(2);
// axes
h21.GetXaxis()->CenterTitle();
h21.GetYaxis()->CenterTitle();
// draw everything
h21.DrawClone("hist");
h22.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Transverse Momentum");
TText* tptr2;
auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t2->SetFillColorAlpha(kWhite, 0);
t2->SetTextFont(43);
t2->SetTextSize(25);
tptr2 = t2->AddText("simulated");
tptr2->SetTextColor(plot::kMpBlue);
tptr2 = t2->AddText("reconstructed");
tptr2->SetTextColor(plot::kMpOrange);
t2->Draw();
//pad 3 phi
c.cd(3);
//gPad->SetLogx(false);
//gPad->SetLogy(false);
auto& h31 = *h_phi_sim;
auto& h32 = *h_phi_rec;
// histogram style
h31.SetLineColor(plot::kMpBlue);
h31.SetLineWidth(2);
h32.SetLineColor(plot::kMpOrange);
h32.SetLineWidth(2);
// axes
h31.GetXaxis()->CenterTitle();
h31.GetYaxis()->CenterTitle();
// draw everything
h31.DrawClone("hist");
h32.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "#phi");
TText* tptr3;
auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t3->SetFillColorAlpha(kWhite, 0);
t3->SetTextFont(43);
t3->SetTextSize(25);
tptr3 = t3->AddText("simulated");
tptr3->SetTextColor(plot::kMpBlue);
tptr3 = t3->AddText("reconstructed");
tptr3->SetTextColor(plot::kMpOrange);
t3->Draw();
//pad 4 rapidity
c.cd(4);
//gPad->SetLogx(false);
//gPad->SetLogy(false);
auto& h41 = *h_y_sim;
auto& h42 = *h_y_rec;
// histogram style
h41.SetLineColor(plot::kMpBlue);
h41.SetLineWidth(2);
h42.SetLineColor(plot::kMpOrange);
h42.SetLineWidth(2);
// axes
h41.GetXaxis()->CenterTitle();
h41.GetYaxis()->CenterTitle();
// draw everything
h41.DrawClone("hist");
h42.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "Rapidity");
TText* tptr4;
auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t4->SetFillColorAlpha(kWhite, 0);
t4->SetTextFont(43);
t4->SetTextSize(25);
tptr4 = t4->AddText("simulated");
tptr4->SetTextColor(plot::kMpBlue);
tptr4 = t4->AddText("reconstructed");
tptr4->SetTextColor(plot::kMpOrange);
t4->Draw();
// Print canvas to output file
c.Print(fmt::format("{}vm_mass.png", output_prefix).c_str());
c.Print(fmt::format("{}vm_mass_pt_phi_rapidity.png", output_prefix).c_str());
}
// TODO we're not actually doing an IM fit yet, so for now just return an
......@@ -136,7 +252,7 @@ int vm_mass(const std::string& config_name) {
// write out our test data
eic::util::write_test(vm_mass_resolution_test,
fmt::format("{}vm_mass.json", output_prefix));
fmt::format("{}vm_mass.json", output_prefix));
// That's all!
return 0;
......
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment