Linear Collider Forum



Home » Simulation » Full Simulations » Converting HepMC to Mokka compatible format
Re: Converting HepMC to Mokka compatible format [message #1955 is a reply to message #1950] Mon, 26 April 2010 00:47 Go to previous messageGo to previous message
weuste
Messages: 2
Registered: April 2010
I had a similar problem. Based on Frank's excerpt from Mokka I wrote a hepevt output formatter.
I seems to work for the long version of hepevt, but the short version (i.e. .HEPEvt) seems to crash Mokka for some reason. Maybe there is still some glitch in here, but it might still be a helpful start for you.

#include "Pythia.h"

#include <iostream>
#include <iomanip>

using namespace Pythia8;
using namespace std;

int main()
{
	Pythia pythia;

	pythia.readString("SLHA:file = sps1a.spc");
	pythia.readString("SUSY:all = on");
	pythia.init("ss_out.lhef");

	pythia.settings.listChanged();
	pythia.particleData.listChanged();

	int fails = 0;
	int nrAllowedFails = 10;
	int evtNr = 0;
	bool longHepEvt = true;

	ofstream out("out.hepevt");

	while (true)
	{
		if (!pythia.next())
		{
			if (pythia.info.atEndOfFile())
			{
				cout << "End of input file reached. Exiting..." << endl;
				break;
			}
			if (++fails > nrAllowedFails)
			{
				cerr << "WARNING: Number allowed fails of Pythia::next() is reached! Will exit program!" << endl;
				break;
			}
			cout << "WARNING: There was an error when calling Pythia::next()! Will try again." << " (" << setw(3) << fails << "/" << nrAllowedFails << " fails)" << endl;
			continue;
		}

		if (++evtNr % 100 == 0)
			cout << "Status: " << setw(6) << evtNr << " (" << setw(3) << fails << "/" << nrAllowedFails << " fails)" << endl;

		out << pythia.event.size() << endl;		// NHEP, nr of particles in this event

		for (int i = 0; i < pythia.event.size(); i++)
		{

			// try new format, i.e. long hepevt format
			// see http://forum.linearcollider.org/index.php?t=msg&th=706&rid=0&S=f812b9be4e054ce14daec834e2b9f2c6
			// for more infos

			// the pythia status code is different from the HepEVT status
			// HepEvt: 0: ignore - 1: final state particle - 2: decayed particle, just included to preserve history - 4+ : reserved/builder specific/user, treat as 0
			out << setw(4) << (pythia.event[i].isFinal() ? 1 : 2) << ' ';	// ISTHEP, status code
			out << setw(8) << pythia.event[i].id() << ' ';					// IDHEP, PDG code
			if (longHepEvt) {
				out << setw(4) << pythia.event[i].mother1()+1 << ' ';			// JMOHEP1, first mother id
				out << setw(4) << pythia.event[i].mother2()+1 << ' ';			// JMOHEP2, last mother id
			}
			out << setw(4) << pythia.event[i].daughter1()+1 << ' ';			// JDAHEP1, first daughter id
			out << setw(4) << pythia.event[i].daughter2()+1 << ' ';			// JDAHEP2, last daughter id

			// now start the info in double, so switch the output to that format
			out.setf(ios_base::scientific);
			out.precision(8);

			out << setw(15) << pythia.event[i].px() << ' ';		// PHEP1, px in GeV/c
			out << setw(15) << pythia.event[i].py() << ' ';		// PHEP2, py in GeV/c
			out << setw(15) << pythia.event[i].pz() << ' ';		// PHEP3, pz in GeV/c
			if (longHepEvt)
				out << setw(15) << pythia.event[i].e() << ' ';		// PHEP4, energy in GeV
			out << setw(15) << pythia.event[i].m() << ' ';		// PHEP5, mass in GeV/cc

			if (longHepEvt) {
				out << setw(15) << pythia.event[i].xProd() << ' ';	// VHEP1, x vertex pos in mm
				out << setw(15) << pythia.event[i].yProd() << ' ';		// VHEP2, y vertex pos in mm
				out << setw(15) << pythia.event[i].zProd() << ' ';		// VHEP3, z vertex pos in mm
				out << setw(15) << pythia.event[i].tProd() << ' ';		// VHEP4, production time in mm/c
			}

			out << endl;

		}
	}
	out.close();

	pythia.statistics();



	return 0;
}
 
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Previous Topic:Converter from Les Houches event xml files to stdhep format?
Next Topic:How to count events in .stdhep file
Goto Forum:
  


Current Time: Sun Dec 16 09:06:31 Pacific Standard Time 2018
.:: Contact :: Home ::.

Powered by: FUDforum 3.0.1.
Copyright ©2001-2010 FUDforum Bulletin Board Software