pythia_dis.cc

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

  // 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 = 2212");
  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;

}