diff --git a/dis/src/pythia_dis.cc b/dis/src/pythia_dis.cc index f450a6c48b263a06a932f8f1ec58be5cba313816..9fed61ac07591dca7d0c683ba03091282d3efcd7 100644 --- a/dis/src/pythia_dis.cc +++ b/dis/src/pythia_dis.cc @@ -1,44 +1,3 @@ -<<<<<<< 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. -// Please respect the MCnet Guidelines, see GUIDELINES for details. - -// Author: Stefan Prestel <stefan.prestel@thep.lu.se>. - -// Keywords: LHE file; hepmc; - -// This program (main45.cc) illustrates how a file with HepMC3 events can be -// generated by Pythia8. See main44.cc for how to ouput HepMC2 events instead. -// Note: both main44.cc and main45.cc can use the same main44.cmnd input card. - -#include "Pythia8/Pythia.h" -#include "Pythia8Plugins/HepMC3.h" -#include <unistd.h> - -using namespace Pythia8; - -//========================================================================== - -// Example main programm to illustrate merging. - -int main( int argc, char* argv[] ){ - - // Check that correct number of command-line arguments - if (argc != 2) { - cerr << " Unexpected number of command-line arguments ("<<argc<<"). \n" - << " You are expected to provide the arguments" << endl - << " 1. Output file for HepMC events" << endl - << " Program stopped. " << endl; - return 1; - } - - // Beam energies, minimal Q2, number of events to generate. - double eProton = 250.; - double eElectron = 10.0; - double Q2min = 5.; - int nEvent = 10000; -======= #include "Pythia8/Pythia.h" #include "Pythia8Plugins/HepMC3.h" #include <unistd.h> @@ -147,13 +106,11 @@ int main(int argc, char* argv[]) { return 0; } - - // Beam energies, minimal Q2, number of events to generate. + // 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; @@ -162,11 +119,7 @@ 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.readString("Beams:idA = " + std::to_string(s.ion_PID)); pythia.settings.parm("Beams:eA", eProton); // BeamB = electron. pythia.readString("Beams:idB = 11"); @@ -176,11 +129,7 @@ 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); @@ -205,23 +154,6 @@ 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); - 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); -======= double Wmax = sqrt(4. * eProton * eElectron); Hist Qhist("Q [GeV]", 100, 0., 50.); Hist Whist("W [GeV]", 100, 0., Wmax); @@ -235,13 +167,14 @@ int main(int argc, char* argv[]) { 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) + 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) { - if (!pythia.next()) continue; + if (!pythia.next()) + continue; double sigmaSample = 0., errorSample = 0.; @@ -252,35 +185,37 @@ int main(int argc, char* argv[]) { 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); + 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() ); + 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() ); - } + 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(); + double evtweight = pythia.info.weight(); // Do not print zero-weight events. - if ( evtweight == 0. ) continue; + 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>(); + genRunInfo = make_shared<HepMC3::GenRunInfo>(); vector<string> weight_names = pythia.info.weightNameVector(); genRunInfo->set_weight_names(weight_names); ascii_io.set_run_info(genRunInfo); @@ -288,52 +223,40 @@ int main(int argc, char* argv[]) { 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); + normhepmc = 1. / double(1e9 * nEvent); // Work with unweighted events. else - normhepmc = xs / double(1e9*nEvent); + normhepmc = xs / double(1e9 * nEvent); // Set event weight - //hepmcevt.weights().push_back(evtweight*normhepmc); + // hepmcevt.weights().push_back(evtweight*normhepmc); // Fill HepMC event - toHepMC.fill_next_event( pythia, &hepmcevt ); + 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); + 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 ); + 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); -<<<<<<< HEAD - - -======= ->>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a // End of event loop. Statistics and histograms. } pythia.stat(); cout << Qhist << Whist << xhist << yhist << pTehist << pTrhist << pTdhist; - return 0; -<<<<<<< HEAD - } -======= - } ->>>>>>> 5c96bf0c8bedda31e7d1c8a8f1d24ceb5722db3a