#include "HepMCRdr.h" #include "GenEvent.h" #include "HepMC/IO_GenEvent.h"//HepMC #include "HepMC/GenEvent.h" #include "HepMC/GenVertex.h" #include "HepMC/Polarization.h" #include "edm4hep/MCParticle.h" //edm4hep #include "edm4hep/MCParticleObj.h" #include "edm4hep/MCParticleCollection.h" #include "edm4hep/EventHeaderCollection.h" #include <iostream> #include <vector> #include <fstream> using namespace edm4hep; using namespace std; DECLARE_COMPONENT(HepMCRdr) HepMCRdr::~HepMCRdr(){ delete ascii_in; } bool HepMCRdr::mutate(MyHepMC::GenEvent& event){ HepMC::GenEvent* evt = ascii_in->read_next_event(); if(!evt) return false; m_processed_event ++; int n_mc = evt->particles_size(); //std::cout<<"Read event :"<< m_processed_event <<", mc size :"<< n_mc <<std::endl; std::map<int, int> pmcid_lmcid; int index = 0 ; for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) { //std::cout<<"start mc "<<index<<std::endl; auto mcp = event.m_mc_vec.create(); pmcid_lmcid.insert(std::pair<int, int>((*p)->barcode(),index)); index++; //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl; mcp.setPDG ((*p)->pdg_id()); mcp.setGeneratorStatus ((*p)->status()); mcp.setSimulatorStatus (999); mcp.setCharge (999); mcp.setTime (999); mcp.setMass ((*p)->generated_mass()); if ( (*p)->production_vertex() ){ HepMC::GenVertex* vertex_pro = (*p)->production_vertex(); double three[3] = {vertex_pro->point3d().x(), vertex_pro->point3d().y(), vertex_pro->point3d().z()}; mcp.setVertex (edm4hep::Vector3d (three)); } else mcp.setVertex(edm4hep::Vector3d()); if ( (*p)->end_vertex() ){ HepMC::GenVertex* vertex_end = (*p)->end_vertex(); double three[3] = {vertex_end->point3d().x(), vertex_end->point3d().y(), vertex_end->point3d().z()}; mcp.setEndpoint (edm4hep::Vector3d (three)); } else mcp.setEndpoint (edm4hep::Vector3d()); mcp.setMomentum (edm4hep::Vector3f(float((*p)->momentum().px()), float((*p)->momentum().py()), float((*p)->momentum().pz()) )); mcp.setMomentumAtEndpoint (edm4hep::Vector3f(float((*p)->momentum().px()), float((*p)->momentum().py()), float((*p)->momentum().pz()) )); const HepMC::Polarization & polar = (*p)->polarization(); mcp.setSpin (edm4hep::Vector3f(polar.normal3d().x(), polar.normal3d().y(), polar.normal3d().z()) ); int two[2] = {1, (*p)->flow(1)}; mcp.setColorFlow (edm4hep::Vector2i (two) ); } // second loop for setting parents and daughters index = 0 ; for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) { auto pmc = event.m_mc_vec.at(index); index++; if ( (*p)->production_vertex() ) { for ( HepMC::GenVertex::particle_iterator mother = (*p)->production_vertex()-> particles_begin(HepMC::parents); mother != (*p)->production_vertex()-> particles_end(HepMC::parents); ++mother ) { pmc.addToParents( event.m_mc_vec.at( pmcid_lmcid.at((*mother)->barcode()) ) ); } } if ( (*p)->end_vertex() ) { for ( HepMC::GenVertex::particle_iterator des =(*p)->end_vertex()-> particles_begin(HepMC::descendants); des != (*p)->end_vertex()-> particles_end(HepMC::descendants); ++des ) { pmc.addToDaughters( event.m_mc_vec.at( pmcid_lmcid.at((*des)->barcode()) ) ); } } } event.SetEventHeader( m_processed_event, -99, 9999, "Generator"); //std::cout<<"end event :"<< m_processed_event <<std::endl; delete evt; return true; } bool HepMCRdr::isEnd(){ return false; } bool HepMCRdr::configure_gentool(){ ascii_in = new HepMC::IO_GenEvent(m_filename.value().c_str(),std::ios::in); m_processed_event=0; return true; } bool HepMCRdr::finish(){ return true; } StatusCode HepMCRdr::initialize() { StatusCode sc; if (not configure_gentool()) { error() << "failed to initialize." << endmsg; return StatusCode::FAILURE; } return sc; } StatusCode HepMCRdr::finalize() { StatusCode sc; if (not finish()) { error() << "Failed to finalize." << endmsg; return StatusCode::FAILURE; } return sc; }