modified: dis/src/pythia_dis.cc

dis/src/pythia_dis.cc

-<<<<<<< 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