Skip to content
Snippets Groups Projects
Commit c5e879e3 authored by lintao@ihep.ac.cn's avatar lintao@ihep.ac.cn
Browse files

migrate the HepEvt reader.

parent 637085ae
No related branches found
No related tags found
No related merge requests found
......@@ -10,8 +10,8 @@ gaudi_add_module(GenAlgo
src/GenReader.cpp
src/StdHepRdr.cpp
src/GenPrinter.cpp
# src/LCAscHepRdr.cc
# src/HepevtRdr.cpp
src/LCAscHepRdr.cc
src/HepevtRdr.cpp
src/SLCIORdr.cpp
src/HepMCRdr.cpp
src/GtGunTool.cpp
......
......@@ -7,12 +7,10 @@
#include "IMPL/MCParticleImpl.h"
#include "plcio/MCParticle.h" //plcio
#include "plcio/MCParticleObj.h"
#include "plcio/MCParticleCollection.h"
#include "plcio/DoubleThree.h"
#include "plcio/FloatThree.h"
#include "plcio/EventHeaderCollection.h"
#include "edm4hep/MCParticle.h" //edm4hep
#include "edm4hep/MCParticleObj.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/EventHeaderCollection.h"
......@@ -23,7 +21,7 @@
using namespace lcio;
using namespace IMPL;
using namespace plcio;
using namespace edm4hep;
using namespace std;
typedef enum HEPFILEFORMATS
......@@ -35,28 +33,57 @@ typedef enum HEPFILEFORMATS
} HEPFILEFORMAT;
HepevtRdr::HepevtRdr(string name){
DECLARE_COMPONENT(HepevtRdr)
m_hepevt_rdr = new UTIL::LCAscHepRdr(name.c_str(), hepevt);
m_processed_event=0;
std::cout<<"initial hepevt_rdr"<<std::endl;
HepevtRdr::~HepevtRdr(){
delete m_hepevt_rdr;
}
HepevtRdr::~HepevtRdr(){
delete m_hepevt_rdr;
StatusCode HepevtRdr::initialize() {
StatusCode sc;
if (not configure_gentool()) {
error() << "failed to initialize." << endmsg;
return StatusCode::FAILURE;
}
return sc;
}
StatusCode HepevtRdr::finalize() {
StatusCode sc;
if (not finish()) {
error() << "Failed to finalize." << endmsg;
return StatusCode::FAILURE;
}
return sc;
}
bool HepevtRdr::configure_gentool(){
int format = hepevt;
if (m_format == "HEPEvt") {
format = HEPEvt;
} else if (m_format == "hepevt") {
format = hepevt;
}
m_hepevt_rdr = new UTIL::LCAscHepRdr(m_filename.value().c_str(), format);
m_processed_event=0;
std::cout<<"initial hepevt_rdr"<<std::endl;
return true;
}
bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
LCCollectionVec* mc_vec = m_hepevt_rdr->readEvent();
if(mc_vec==nullptr) return false;
m_processed_event ++;
int n_mc = mc_vec->getNumberOfElements();
int n_mc = mc_vec->size();
std::cout<<"Read event :"<< m_processed_event <<", mc size :"<< n_mc <<std::endl;
std::map<int, int> pmcid_lmcid;
for (int i=0; i < n_mc; i++){
MCParticleImpl* mc = (MCParticleImpl*) mc_vec->getElementAt(i);
//std::cout<<"At mc :"<< i <<std::endl;
plcio::MCParticle mcp = event.m_mc_vec.create();
// std::cout<<"At mc :"<< i <<std::endl;
auto mcp = event.m_mc_vec.create();
pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i));
//std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl;
......@@ -68,8 +95,8 @@ bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
mcp.setMass (mc->getMass());
mcp.setVertex (mc->getVertex());
mcp.setEndpoint (mc->getEndpoint());
mcp.setMomentum (FloatThree(float(mc->getMomentum()[0]), float(mc->getMomentum()[1]), float(mc->getMomentum()[2]) ));
mcp.setMomentumAtEndpoint (FloatThree(float(mc->getMomentumAtEndpoint()[0]), float(mc->getMomentumAtEndpoint()[1]), float(mc->getMomentumAtEndpoint()[2]) ));
mcp.setMomentum (edm4hep::Vector3f(float(mc->getMomentum()[0]), float(mc->getMomentum()[1]), float(mc->getMomentum()[2]) ));
mcp.setMomentumAtEndpoint (edm4hep::Vector3f(float(mc->getMomentumAtEndpoint()[0]), float(mc->getMomentumAtEndpoint()[1]), float(mc->getMomentumAtEndpoint()[2]) ));
mcp.setSpin (mc->getSpin());
mcp.setColorFlow (mc->getColorFlow());
}
......@@ -79,16 +106,16 @@ bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
MCParticleImpl* mc = (MCParticleImpl*) mc_vec->getElementAt(i);
const MCParticleVec & mc_parents = mc->getParents();
const MCParticleVec & mc_daughters = mc->getDaughters();
plcio::MCParticle pmc = event.m_mc_vec.at(i);
auto pmc = event.m_mc_vec.at(i);
//std::cout<<"mc at "<< i<<", parent size "<<mc_parents.size() <<std::endl;
for(unsigned int j=0; j< mc_parents.size(); j++){int p_id = mc_parents.at(j)->id();
//std::cout<<"parent id "<<p_id<<std::endl;
pmc.addParent( event.m_mc_vec.at( pmcid_lmcid.at(p_id) ) );
pmc.addToParents( event.m_mc_vec.at( pmcid_lmcid.at(p_id) ) );
}
//std::cout<<"mc at "<< i<<", daughter size "<<mc_daughters.size() <<std::endl;
for(unsigned int j=0; j< mc_daughters.size(); j++){int d_id = mc_daughters.at(j)->id();
//std::cout<<"daughter id "<<d_id<<std::endl;
pmc.addDaughter( event.m_mc_vec.at( pmcid_lmcid.at(d_id) ) );
pmc.addToDaughters( event.m_mc_vec.at( pmcid_lmcid.at(d_id) ) );
}
}
......@@ -98,13 +125,9 @@ bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
}
bool HepevtRdr::isEnd(){
return false;
}
bool HepevtRdr::configure(){
return true;
return false;
}
bool HepevtRdr::finish(){
return true;
return true;
}
#ifndef HepevtRdr_h
#define HepevtRdr_h 1
#include "GaudiKernel/AlgTool.h"
#include "GenReader.h"
#include "GenEvent.h"
......@@ -8,19 +10,27 @@
#include "EVENT/LCIO.h"
#include "LCAscHepRdr.h"
class HepevtRdr: public GenReader{
class HepevtRdr: public extends<AlgTool, GenReader> {
public:
HepevtRdr(string name);
using extends::extends;
~HepevtRdr();
bool configure();
StatusCode initialize() override;
StatusCode finalize() override;
bool configure_gentool();
bool mutate(MyHepMC::GenEvent& event);
bool finish();
bool isEnd();
private:
UTIL::LCAscHepRdr* m_hepevt_rdr;
long m_total_event;
long m_processed_event;
UTIL::LCAscHepRdr* m_hepevt_rdr = nullptr;
long m_total_event = -1;
long m_processed_event = -1;
// input file name
Gaudi::Property<std::string> m_filename{this, "Input"};
Gaudi::Property<std::string> m_format{this, "Format"};
};
#endif
......
......@@ -56,11 +56,16 @@ namespace UTIL{
//
// Read the event, check for errors
//
int NHEP; // number of entries
int NOUT ; // number of outgoing particles
int BRE ; // beam remnants
double WEIGHT ; // weight
inputFile >> NHEP >> NOUT >> BRE >> WEIGHT;
int NHEP = -1; // number of entries
int NOUT = -1; // number of outgoing particles
int BRE = -1; // beam remnants
double WEIGHT = -1; // weight
std::string line; // modified by Tao
std::getline(inputFile, line);
std::stringstream ss_(line);
ss_ >> NHEP >> NOUT >> BRE >> WEIGHT;
std::cout << "NHEP: " << NHEP << std::endl;
if( inputFile.eof() )
{
//
......@@ -74,6 +79,7 @@ namespace UTIL{
// Create a Collection Vector
//
mcVec = new IMPL::LCCollectionVec(LCIO::MCPARTICLE);
std::cout << "mc size: " << mcVec->size() << std::endl;
MCParticleImpl* p;
MCParticleImpl* d;
......@@ -101,12 +107,15 @@ namespace UTIL{
for( int IHEP=0; IHEP<NHEP; IHEP++ )
{
//if ( theFileFormat == HEPEvt)
if ( false)
inputFile >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
std::getline(inputFile, line);
std::cout << "LINE: " << line << std::endl;
std::stringstream ss(line);
if ( theFileFormat == HEPEvt)
ss >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
>> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5;
else
inputFile >> ISTHEP >> IDHEP
ss >> ISTHEP >> IDHEP
>> JMOHEP1 >> JMOHEP2
>> JDAHEP1 >> JDAHEP2
>> PHEP1 >> PHEP2 >> PHEP3
......@@ -114,6 +123,8 @@ namespace UTIL{
>> VHEP1 >> VHEP2 >> VHEP3
>> VHEP4;
std::cout << "ISTHEP: " << ISTHEP << std::endl;
if(inputFile.eof())
return nullptr;
//
......@@ -124,6 +135,7 @@ namespace UTIL{
// PDGID
//
mcp->setPDG(IDHEP);
std::cout << "PDG: " << IDHEP << std::endl;
//
// Momentum vector
//
......@@ -221,6 +233,7 @@ namespace UTIL{
//
for( int IHEP=0; IHEP<NHEP; IHEP++ )
{
continue;
//
// Get the MCParticle
//
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment