From 603368018da7da74b59c63d11bf89ee2a98d0163 Mon Sep 17 00:00:00 2001
From: Wenxing Fang <fangwx@ihep.ac.cn>
Date: Wed, 28 Aug 2019 15:11:36 +0800
Subject: [PATCH] V3
---
.gitignore | 2 +
Examples/options/gen_write.py | 33 ++++
Generator/CMakeLists.txt | 42 ++++-
Generator/src/GenAlgo.cpp | 92 ++++++++++
Generator/src/GenAlgo.h | 49 +++++
Generator/src/GenEvent.cpp | 50 +++++
Generator/src/GenEvent.h | 31 ++++
Generator/src/GenPrinter.cpp | 42 +++++
Generator/src/GenPrinter.h | 19 ++
Generator/src/GenReader.cpp | 5 +
Generator/src/GenReader.h | 19 ++
Generator/src/GenWriter.cpp | 44 +++++
Generator/src/GenWriter.h | 35 ++++
Generator/src/HepMCRdr.cpp | 113 ++++++++++++
Generator/src/HepMCRdr.h | 27 +++
Generator/src/HepevtRdr.cpp | 110 +++++++++++
Generator/src/HepevtRdr.h | 27 +++
Generator/src/IGenTool.cpp | 5 +
Generator/src/IGenTool.h | 15 ++
Generator/src/LCAscHepRdr.cc | 333 ++++++++++++++++++++++++++++++++++
Generator/src/LCAscHepRdr.h | 40 ++++
Generator/src/SLCIORdr.cpp | 169 +++++++++++++++++
Generator/src/SLCIORdr.h | 32 ++++
Generator/src/StdHepRdr.cpp | 104 +++++++++++
Generator/src/StdHepRdr.h | 28 +++
25 files changed, 1465 insertions(+), 1 deletion(-)
create mode 100644 .gitignore
create mode 100644 Examples/options/gen_write.py
create mode 100644 Generator/src/GenAlgo.cpp
create mode 100644 Generator/src/GenAlgo.h
create mode 100644 Generator/src/GenEvent.cpp
create mode 100644 Generator/src/GenEvent.h
create mode 100644 Generator/src/GenPrinter.cpp
create mode 100644 Generator/src/GenPrinter.h
create mode 100644 Generator/src/GenReader.cpp
create mode 100644 Generator/src/GenReader.h
create mode 100644 Generator/src/GenWriter.cpp
create mode 100644 Generator/src/GenWriter.h
create mode 100644 Generator/src/HepMCRdr.cpp
create mode 100644 Generator/src/HepMCRdr.h
create mode 100644 Generator/src/HepevtRdr.cpp
create mode 100644 Generator/src/HepevtRdr.h
create mode 100644 Generator/src/IGenTool.cpp
create mode 100644 Generator/src/IGenTool.h
create mode 100644 Generator/src/LCAscHepRdr.cc
create mode 100644 Generator/src/LCAscHepRdr.h
create mode 100644 Generator/src/SLCIORdr.cpp
create mode 100644 Generator/src/SLCIORdr.h
create mode 100644 Generator/src/StdHepRdr.cpp
create mode 100644 Generator/src/StdHepRdr.h
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..80e594d9
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+./Generator/output/
+./Generator/options/
diff --git a/Examples/options/gen_write.py b/Examples/options/gen_write.py
new file mode 100644
index 00000000..001dd825
--- /dev/null
+++ b/Examples/options/gen_write.py
@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+
+from Gaudi.Configuration import *
+
+from Configurables import CEPCDataSvc
+dsvc = CEPCDataSvc("EventDataSvc")
+
+from Configurables import GenAlgo
+
+read = GenAlgo("read")
+#######################################
+#support format: stdhep, slcio, hepmc #
+#######################################
+#read.Input = "/junofs/users/wxfang/CEPC/CEPCOFF/TestExample/stdhep/nnh_e2e2.e0.p0.00001.stdhep"
+#read.FileFormat = "stdhep"
+read.Input = "/junofs/users/wxfang/CEPC/whizard_apply/ee/ee.slcio"
+read.FileFormat = "slcio"
+read.PrintEvent = True # true for printing mc info
+read.WriteFile = True # true for writting info to root
+
+from Configurables import PodioOutput
+out = PodioOutput("out")
+out.filename = "test.root" #name of output root file
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [read, out],
+ EvtSel = 'NONE',
+ EvtMax = 11,
+ ExtSvc=[dsvc],
+ OutputLevel=DEBUG
+)
diff --git a/Generator/CMakeLists.txt b/Generator/CMakeLists.txt
index ef5f3728..a9a5cd9f 100644
--- a/Generator/CMakeLists.txt
+++ b/Generator/CMakeLists.txt
@@ -1,2 +1,42 @@
+########################################
+gaudi_subdir(Generator v0r0)
+set(GenAlgo_srcs src/GenAlgo.cpp src/GenEvent.cpp src/IGenTool.cpp src/StdHepRdr.cpp src/GenReader.cpp src/GenPrinter.cpp src/LCAscHepRdr.cc src/HepevtRdr.cpp src/SLCIORdr.cpp src/HepMCRdr.cpp)
+set(GenAlgo_incs src)
-gaudi_subdir(Generator v0r0)
\ No newline at end of file
+find_package(ROOT COMPONENTS RIO Tree TreePlayer MathCore Net Graf3d Graf Gpad REQUIRED)
+find_package(LCIO)
+find_package(podio)
+find_package(plcio)
+find_package(HepMC)
+if(ROOT_FOUND)
+ message("found ROOT: ${ROOT_INCLUDE_DIRS} ${ROOT_LIBRARIES}")
+endif(ROOT_FOUND)
+if(LCIO_FOUND)
+ message("found LCIO: ${LCIO_INCLUDE_DIRS} ${LCIO_LIBRARIES}")
+endif(LCIO_FOUND)
+if(podio_FOUND)
+ message("found podio: ${podio_INCLUDE_DIRS} ${podio_LIBRARIES}")
+endif(podio_FOUND)
+if(plcio_FOUND)
+ message("found plcio: ${plcio_INCLUDE_DIRS} ${plcio_LIBRARY_DIR}")
+endif(plcio_FOUND)
+if(HepMC_FOUND)
+ message("found HepMC: ${HepMC_INCLUDE_DIRS} ${HepMC_LIBRARY_DIR}")
+endif(HepMC_FOUND)
+############## for producing plcio library #############
+INCLUDE_DIRECTORIES(${GenAlgo_incs})
+
+gaudi_add_module(GenAlgo ${GenAlgo_srcs}
+INCLUDE_DIRS
+GaudiKernel
+FWCore
+${LCIO_INCLUDE_DIRS}
+${podio_INCLUDE_DIRS}
+${plcio_INCLUDE_DIRS}
+${ROOT_INCLUDE_DIRS}
+HepMC
+LINK_LIBRARIES GaudiKernel ${LCIO_LIBRARIES} ${podio_LIBRARIES} ROOT ${plcio_LIBRARY_DIR}/libplcio.so ${plcio_LIBRARY_DIR}/libplcioDict.so FWCore HepMC
+)
+#gaudi_add_test(Reader FRAMEWORK options/read.py)
+
+###########################
diff --git a/Generator/src/GenAlgo.cpp b/Generator/src/GenAlgo.cpp
new file mode 100644
index 00000000..0d8c4698
--- /dev/null
+++ b/Generator/src/GenAlgo.cpp
@@ -0,0 +1,92 @@
+#include "GenAlgo.h"
+
+#include "GaudiKernel/IEventProcessor.h"
+#include "GaudiKernel/IAppMgrUI.h"
+#include "GaudiKernel/GaudiException.h"
+
+
+#include "plcio/MCParticleCollection.h"//plico
+
+#include <iostream>
+#include <vector>
+#include <fstream>
+
+#include "IGenTool.h"
+#include "GenEvent.h"
+#include "StdHepRdr.h"
+#include "HepevtRdr.h"// not correct still
+#include "SLCIORdr.h"
+#include "HepMCRdr.h"
+#include "GenPrinter.h"
+#include "GenWriter.h"
+
+using namespace std;
+
+DECLARE_COMPONENT(GenAlgo)
+
+GenAlgo::GenAlgo(const std::string& name, ISvcLocator* pSvcLocator): GaudiAlgorithm(name, pSvcLocator) {
+ declareProperty("MCParticleCol", m_hdl, "MCParticle collection (output)");
+ m_evtid = 0;
+
+}
+
+StatusCode
+GenAlgo::initialize() {
+
+ cout << "initialize start" << endl;
+ string generatorName = m_input_file.value();
+ string outputName = m_output_file.value();
+ string format = m_input_format.value();
+ IGenTool* gen_reader;
+ if(format=="stdhep") gen_reader = new StdHepRdr(generatorName);
+ else if(format=="slcio") gen_reader = new SLCIORdr(generatorName);
+ else if(format=="hepmc") gen_reader = new HepMCRdr(generatorName);
+ else{cout << "Error : unsupport format for generator input file" << endl; return StatusCode::FAILURE; }
+ //IGenTool* gen_reader = new HepevtRdr(generatorName);
+ m_genTools.push_back(gen_reader);
+ if(m_print.value()) {
+ IGenTool* gen_printer = new GenPrinter(generatorName);
+ m_genTools.push_back(gen_printer);
+ }
+ //IGenTool* gen_writer = new GenWriter (outputName);
+ //m_genTools.push_back(gen_writer);
+
+ cout << "initialize done" << endl;
+ return StatusCode::SUCCESS;
+
+}
+
+StatusCode
+GenAlgo::execute() {
+ m_evtid++;
+ auto mcCol = m_hdl.createAndPut();
+ MyHepMC::GenEvent m_event(*mcCol);
+
+ for(std::vector<IGenTool*>::iterator it=m_genTools.begin(); it != m_genTools.end(); ++it) {
+ if ((*it)->mutate(m_event)) {}
+ else {
+ cout << "Have read all events, stop now." << endl;
+ auto ep = serviceLocator()->as<IEventProcessor>();
+ if ( !ep ) {
+ error() << "Cannot get IEventProcessor" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ ep->stopRun();
+ return StatusCode::SUCCESS;
+
+ }
+ }
+
+ return StatusCode::SUCCESS;
+
+}
+
+StatusCode
+GenAlgo::finalize() {
+ cout << "finalize" << endl;
+ for(std::vector<IGenTool*>::iterator it=m_genTools.begin(); it != m_genTools.end(); ++it) {
+ if ((*it)->finish()) {}
+ else {cout << "finish Failed" << endl; return StatusCode::FAILURE; }
+ }
+ return StatusCode::SUCCESS;
+}
diff --git a/Generator/src/GenAlgo.h b/Generator/src/GenAlgo.h
new file mode 100644
index 00000000..43b4339b
--- /dev/null
+++ b/Generator/src/GenAlgo.h
@@ -0,0 +1,49 @@
+#ifndef GenAlgo_h
+#define GenAlgo_h
+
+
+#include <GaudiKernel/Algorithm.h>
+#include "GaudiKernel/Property.h"
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "FWCore/DataHandle.h"
+
+#include "GenEvent.h"
+
+class IGenTool;
+namespace plcio {
+ class MCParticleCollection;
+}
+
+
+using namespace std;
+
+class GenAlgo: public GaudiAlgorithm {
+
+ friend class AlgFactory<GenAlgo>;
+public:
+ GenAlgo(const std::string& name, ISvcLocator* pSvcLocator);
+
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+private:
+ Gaudi::Property<std::string> m_input_file{this, "Input", "NULL"};
+ Gaudi::Property<std::string> m_input_format{this, "FileFormat", "NULL"};
+ Gaudi::Property<std::string> m_output_file{this, "OutputRootFile", "NULL"};
+ Gaudi::Property<bool> m_print{this, "PrintEvent", "NULL"};
+ Gaudi::Property<bool> m_do_write{this, "WriteFile", "NULL"};
+
+ std::vector<std::string> m_genToolNames;
+ std::vector<IGenTool*> m_genTools;
+ int m_evtid;
+ int m_evtMax;
+ //MyHepMC::GenEvent m_event;
+ DataHandle<plcio::MCParticleCollection> m_hdl{"MCParticleCol", Gaudi::DataHandle::Writer, this};
+
+
+};
+
+
+#endif
diff --git a/Generator/src/GenEvent.cpp b/Generator/src/GenEvent.cpp
new file mode 100644
index 00000000..1cd5ccb3
--- /dev/null
+++ b/Generator/src/GenEvent.cpp
@@ -0,0 +1,50 @@
+#include "GenEvent.h"
+#include "plcio/MCParticleCollection.h"//plico
+
+using namespace std;
+
+namespace MyHepMC{
+
+//GenEvent::GenEvent(){
+GenEvent::GenEvent(plcio::MCParticleCollection& mcCol)
+ : m_mc_vec(mcCol){
+
+ m_event_id=-1;
+ m_run_id=-1;
+ m_time=-1;
+ m_det_name="";
+
+}
+GenEvent::~GenEvent(){}
+
+void GenEvent::SetEventHeader(long event_id_, long run_id_, float time_, string det_name_){
+ m_event_id = event_id_;
+ m_run_id = run_id_;
+ m_time = time_;
+ m_det_name = det_name_;
+}
+/*
+void GenEvent::SetMCCollection(plcio::MCParticleCollection vec_){
+m_mc_vec = vec_;
+}
+*/
+plcio::MCParticleCollection GenEvent::getMCVec(){
+return m_mc_vec;
+}
+
+long GenEvent::getID(){
+return m_event_id;
+}
+long GenEvent::getRun() {return m_run_id;}
+long GenEvent::getTime() {return m_time;}
+std::string GenEvent::getName() {return m_det_name;}
+
+void GenEvent::ReSet(){
+
+ m_event_id=-1;
+ m_run_id=-1;
+ m_time=-1;
+ m_det_name="";
+ m_mc_vec.clear();
+}
+}
diff --git a/Generator/src/GenEvent.h b/Generator/src/GenEvent.h
new file mode 100644
index 00000000..dae1b377
--- /dev/null
+++ b/Generator/src/GenEvent.h
@@ -0,0 +1,31 @@
+#ifndef GenEvent_h
+#define GenEvent_h 1
+
+#include "plcio/MCParticleCollection.h"//plico
+
+namespace MyHepMC {
+
+class GenEvent{
+ public:
+ //GenEvent();
+ GenEvent(plcio::MCParticleCollection& mcCol);
+ ~GenEvent();
+ void SetEventHeader(long event_id_, long run_id_, float time_, std::string det_name_);
+ //void SetMCCollection(plcio::MCParticleCollection vec_);
+ long getID();
+ long getRun();
+ long getTime();
+ void ReSet();
+ std::string getName();
+ plcio::MCParticleCollection getMCVec();
+ plcio::MCParticleCollection& m_mc_vec;
+ //plcio::MCParticleCollection m_mc_vec;
+ private:
+ long m_event_id;
+ long m_run_id;
+ float m_time;
+ std::string m_det_name;
+};
+
+}
+#endif
diff --git a/Generator/src/GenPrinter.cpp b/Generator/src/GenPrinter.cpp
new file mode 100644
index 00000000..827f5b40
--- /dev/null
+++ b/Generator/src/GenPrinter.cpp
@@ -0,0 +1,42 @@
+#include "GenPrinter.h"
+#include "GenEvent.h"
+
+GenPrinter::GenPrinter(string name){}
+
+
+GenPrinter::~GenPrinter(){
+}
+
+bool GenPrinter::mutate(MyHepMC::GenEvent& event){
+ std::cout << "print mc info for event "<< event.getID() << ", mc size ="<< event.m_mc_vec.size() << std::endl;
+ for ( int i =0; i < event.m_mc_vec.size(); i++ ) {
+ auto p = event.m_mc_vec.at(i);
+ std::cout<< "PDG :"<< p.getPDG ()<<std::endl
+ << "id :"<< p.id ()<<std::endl
+ << "ID :"<< p.getObjectID().index <<std::endl
+ << "GeneratorStatus :"<< p.getGeneratorStatus ()<<std::endl
+ << "SimulatorStatus :"<< p.getSimulatorStatus ()<<std::endl
+ << "Charge :"<< p.getCharge ()<<std::endl
+ << "Time :"<< p.getTime ()<<std::endl
+ << "Mass :"<< p.getMass ()<<std::endl
+ << "Vertex :"<< p.getVertex ()[0]<<std::endl
+ << "Endpoint :"<< p.getEndpoint ()[1]<<std::endl
+ << "Momentum :"<< p.getMomentum ()[2]<<std::endl
+ << "MomentumAtEndpoint:"<< p.getMomentumAtEndpoint()[0]<<std::endl
+ << "Spin :"<< p.getSpin ()[1]<<std::endl
+ << "ColorFlow :"<< p.getColorFlow ()[1]<<std::endl
+ << "Parent size :"<< p.parents_size ()<<std::endl
+ << "Daughter size :"<< p.daughters_size ()<<std::endl;
+ //for(unsigned int j=0; j<p.parents_size(); j++) std::cout << " for parent: "<< j << ",PDG="<< p.getParents(j).getPDG() << ",id=:"<< p.getParents(j).id()<<std::endl;
+ for (auto it = p.parents_begin(), end = p.parents_end(); it != end ; ++it ) std::cout << " for parent ,PDG="<< it->getPDG() << ",id=:"<< it->getObjectID().index<<std::endl;
+ }
+ return true;
+}
+
+bool GenPrinter::configure(){
+return true;
+}
+
+bool GenPrinter::finish(){
+return true;
+}
diff --git a/Generator/src/GenPrinter.h b/Generator/src/GenPrinter.h
new file mode 100644
index 00000000..9fe8a1c9
--- /dev/null
+++ b/Generator/src/GenPrinter.h
@@ -0,0 +1,19 @@
+#ifndef GenPrinter_h
+#define GenPrinter_h 1
+
+#include "GenEvent.h"
+#include "IGenTool.h"
+
+using namespace std;
+
+class GenPrinter: public IGenTool{
+
+ public:
+ GenPrinter(string name);
+ ~GenPrinter();
+ bool configure() override;
+ bool mutate(MyHepMC::GenEvent& event) override;
+ bool finish() override;
+};
+
+#endif
diff --git a/Generator/src/GenReader.cpp b/Generator/src/GenReader.cpp
new file mode 100644
index 00000000..04c38c53
--- /dev/null
+++ b/Generator/src/GenReader.cpp
@@ -0,0 +1,5 @@
+#include "GenReader.h"
+
+GenReader::~GenReader(){
+}
+
diff --git a/Generator/src/GenReader.h b/Generator/src/GenReader.h
new file mode 100644
index 00000000..b9abac0a
--- /dev/null
+++ b/Generator/src/GenReader.h
@@ -0,0 +1,19 @@
+#ifndef GenReader_h
+#define GenReader_h 1
+
+#include "GenEvent.h"
+#include "IGenTool.h"
+
+using namespace std;
+
+class GenReader: public IGenTool{
+
+ public:
+ ~GenReader();
+ virtual bool configure()=0;
+ virtual bool mutate(MyHepMC::GenEvent& event)=0;
+ virtual bool finish()=0;
+ virtual bool isEnd()=0;
+};
+
+#endif
diff --git a/Generator/src/GenWriter.cpp b/Generator/src/GenWriter.cpp
new file mode 100644
index 00000000..f1f93a57
--- /dev/null
+++ b/Generator/src/GenWriter.cpp
@@ -0,0 +1,44 @@
+#include "GenWriter.h"
+#include "GenEvent.h"
+
+#include "podio/EventStore.h" //podio
+#include "podio/ROOTWriter.h"
+
+
+#include "plcio/MCParticleCollection.h"//plico
+#include "plcio/EventHeaderCollection.h"
+
+
+GenWriter::GenWriter(string name){
+ m_output_name = name;
+ store = new podio::EventStore();
+ writer = new podio::ROOTWriter(m_output_name, store);
+ ehc = &store->create<plcio::EventHeaderCollection>("EvtHeaders");
+ mcc = &store->create<plcio::MCParticleCollection>("MCParticles");
+
+ writer->registerForWrite("EvtHeaders");
+ writer->registerForWrite("MCParticles");
+}
+
+GenWriter::~GenWriter(){
+}
+
+bool GenWriter::mutate(MyHepMC::GenEvent& event){
+ std::cout << "write mc info for event "<< event.getID() << ", mc size ="<< event.m_mc_vec.size() << std::endl;
+ mcc=&event.m_mc_vec;
+ auto header = plcio::EventHeader(event.getID(), event.getRun(), event.getTime(), event.getName());
+ ehc->push_back(header);
+ writer->writeEvent();
+ store->clearCollections();
+ return true;
+}
+
+bool GenWriter::configure(){
+return true;
+}
+
+bool GenWriter::finish(){
+ writer->finish();
+ std::cout<<"Saved root "<<m_output_name<<std::endl;
+ return true;
+}
diff --git a/Generator/src/GenWriter.h b/Generator/src/GenWriter.h
new file mode 100644
index 00000000..30e6db43
--- /dev/null
+++ b/Generator/src/GenWriter.h
@@ -0,0 +1,35 @@
+#ifndef GenWriter_h
+#define GenWriter_h 1
+
+#include "GenEvent.h"
+#include "IGenTool.h"
+
+#include "podio/EventStore.h" //podio
+#include "podio/ROOTWriter.h"
+
+
+#include "plcio/MCParticleCollection.h"//plico
+#include "plcio/EventHeaderCollection.h"
+
+
+
+using namespace std;
+
+class GenWriter: public IGenTool{
+
+ public:
+ GenWriter(string name);
+ ~GenWriter();
+ bool configure() override;
+ bool mutate(MyHepMC::GenEvent& event) override;
+ bool finish() override;
+ private:
+ string m_output_name;
+ podio::EventStore* store ;
+ podio::ROOTWriter* writer;
+ plcio::EventHeaderCollection* ehc ;
+ plcio::MCParticleCollection* mcc ;
+
+};
+
+#endif
diff --git a/Generator/src/HepMCRdr.cpp b/Generator/src/HepMCRdr.cpp
new file mode 100644
index 00000000..22acd3ff
--- /dev/null
+++ b/Generator/src/HepMCRdr.cpp
@@ -0,0 +1,113 @@
+#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 "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 <iostream>
+#include <vector>
+#include <fstream>
+
+
+using namespace plcio;
+using namespace std;
+
+
+HepMCRdr::HepMCRdr(string name){
+
+ascii_in = new HepMC::IO_GenEvent(name.c_str(),std::ios::in);
+
+m_processed_event=0;
+}
+
+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;
+ plcio::MCParticle 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 (plcio::DoubleThree (three));
+ }
+ else mcp.setVertex(plcio::DoubleThree());
+ 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 (plcio::DoubleThree (three));
+ }
+ else mcp.setEndpoint (plcio::DoubleThree());
+ mcp.setMomentum (plcio::FloatThree(float((*p)->momentum().px()), float((*p)->momentum().py()), float((*p)->momentum().pz()) ));
+ mcp.setMomentumAtEndpoint (plcio::FloatThree(float((*p)->momentum().px()), float((*p)->momentum().py()), float((*p)->momentum().pz()) ));
+ const HepMC::Polarization & polar = (*p)->polarization();
+ mcp.setSpin (plcio::FloatThree(polar.normal3d().x(), polar.normal3d().y(), polar.normal3d().z()) );
+ int two[2] = {1, (*p)->flow(1)};
+ mcp.setColorFlow (plcio::IntTwo (two) );
+ }
+ // second loop for setting parents and daughters
+ index = 0 ;
+ for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) {
+ plcio::MCParticle 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.addParent( 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.addDaughter( 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(){
+return true;
+}
+
+bool HepMCRdr::finish(){
+return true;
+}
diff --git a/Generator/src/HepMCRdr.h b/Generator/src/HepMCRdr.h
new file mode 100644
index 00000000..61c777bb
--- /dev/null
+++ b/Generator/src/HepMCRdr.h
@@ -0,0 +1,27 @@
+#ifndef HepMCRdr_h
+#define HepMCRdr_h 1
+
+#include "GenReader.h"
+#include "GenEvent.h"
+
+#include "HepMC/IO_GenEvent.h"//HepMC
+#include "HepMC/GenEvent.h"
+
+
+class HepMCRdr: public GenReader{
+
+ public:
+ HepMCRdr(string name);
+ ~HepMCRdr();
+ bool configure();
+ bool mutate(MyHepMC::GenEvent& event);
+ bool finish();
+ bool isEnd();
+ private:
+ HepMC::IO_GenEvent *ascii_in;
+ long m_total_event;
+ long m_processed_event;
+};
+
+#endif
+
diff --git a/Generator/src/HepevtRdr.cpp b/Generator/src/HepevtRdr.cpp
new file mode 100644
index 00000000..7f036e2e
--- /dev/null
+++ b/Generator/src/HepevtRdr.cpp
@@ -0,0 +1,110 @@
+#include "HepevtRdr.h"
+#include "GenEvent.h"
+
+#include "lcio.h" //LCIO
+#include "EVENT/LCIO.h"
+#include "LCAscHepRdr.h"
+#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 <iostream>
+#include <vector>
+#include <fstream>
+
+
+using namespace lcio;
+using namespace IMPL;
+using namespace plcio;
+using namespace std;
+
+typedef enum HEPFILEFORMATS
+{
+ stdhep = 0,
+ HEPEvt,
+ hepevt,
+ slcio
+} HEPFILEFORMAT;
+
+
+HepevtRdr::HepevtRdr(string name){
+
+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;
+}
+
+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();
+ 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();
+ pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i));
+ //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl;
+
+ mcp.setPDG (mc->getPDG());
+ mcp.setGeneratorStatus (mc->getGeneratorStatus());
+ mcp.setSimulatorStatus (mc->getSimulatorStatus());
+ mcp.setCharge (mc->getCharge());
+ mcp.setTime (mc->getTime());
+ 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.setSpin (mc->getSpin());
+ mcp.setColorFlow (mc->getColorFlow());
+ }
+ // second loop for setting parents and daughters
+
+ for (int i=0; i < n_mc; i++){
+ 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);
+ //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) ) );
+ }
+ //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) ) );
+ }
+ }
+
+ event.SetEventHeader( m_processed_event, -99, 9999, "Generator");
+ //std::cout<<"end event :"<< m_processed_event <<std::endl;
+ return true;
+}
+
+bool HepevtRdr::isEnd(){
+return false;
+}
+
+bool HepevtRdr::configure(){
+return true;
+}
+
+bool HepevtRdr::finish(){
+return true;
+}
diff --git a/Generator/src/HepevtRdr.h b/Generator/src/HepevtRdr.h
new file mode 100644
index 00000000..9ee2d7fd
--- /dev/null
+++ b/Generator/src/HepevtRdr.h
@@ -0,0 +1,27 @@
+#ifndef HepevtRdr_h
+#define HepevtRdr_h 1
+
+#include "GenReader.h"
+#include "GenEvent.h"
+
+#include "lcio.h"
+#include "EVENT/LCIO.h"
+#include "LCAscHepRdr.h"
+
+class HepevtRdr: public GenReader{
+
+ public:
+ HepevtRdr(string name);
+ ~HepevtRdr();
+ bool configure();
+ bool mutate(MyHepMC::GenEvent& event);
+ bool finish();
+ bool isEnd();
+ private:
+ UTIL::LCAscHepRdr* m_hepevt_rdr;
+ long m_total_event;
+ long m_processed_event;
+};
+
+#endif
+
diff --git a/Generator/src/IGenTool.cpp b/Generator/src/IGenTool.cpp
new file mode 100644
index 00000000..5562baf5
--- /dev/null
+++ b/Generator/src/IGenTool.cpp
@@ -0,0 +1,5 @@
+#include "IGenTool.h"
+
+IGenTool::~IGenTool() {
+
+}
diff --git a/Generator/src/IGenTool.h b/Generator/src/IGenTool.h
new file mode 100644
index 00000000..87af6bd8
--- /dev/null
+++ b/Generator/src/IGenTool.h
@@ -0,0 +1,15 @@
+#ifndef IGenTool_h
+#define IGenTool_h 1
+
+#include "GenEvent.h"
+
+
+class IGenTool {
+ public:
+ virtual bool mutate(MyHepMC::GenEvent& event)=0;
+ virtual bool finish()=0;
+ virtual bool configure()=0;
+ virtual ~IGenTool();
+};
+
+#endif
diff --git a/Generator/src/LCAscHepRdr.cc b/Generator/src/LCAscHepRdr.cc
new file mode 100644
index 00000000..eabca0bb
--- /dev/null
+++ b/Generator/src/LCAscHepRdr.cc
@@ -0,0 +1,333 @@
+#include "LCAscHepRdr.h"
+#include "EVENT/MCParticle.h"
+#include "IMPL/MCParticleImpl.h"
+#include "lcio.h"
+#include "EVENT/LCIO.h"
+#include <sstream>
+#include "Exceptions.h"
+
+using namespace EVENT ;
+using namespace IMPL ;
+
+//
+// the time being we leave it as #define, to avoid to depend
+// on Mokka business if it could be part of lcio package
+//
+#define HEPEvt 1
+#define hepevt 2
+
+namespace UTIL{
+
+ LCAscHepRdr::LCAscHepRdr(const char* evfile, int fileFormat)
+ : theFileFormat(fileFormat)
+ {
+ //
+ // Adapt the Geant4 G4HEPEvtInterface code
+ //
+ switch(fileFormat)
+ {
+ case HEPEvt :
+ case hepevt :
+ inputFile.open(evfile);
+ if (!inputFile)
+ {
+ std::stringstream description ;
+ description << "LCAscHepRdr, no ascii Hep file found: " << evfile << std::ends ;
+ throw IO::IOException( description.str() );
+ }
+ break;
+ default :
+ std::stringstream description ;
+ description << "LCAscHepRdr, bad suffix on file name: " << evfile << std::ends ;
+ throw IO::IOException( description.str() );
+ }
+ }
+
+ LCAscHepRdr::~LCAscHepRdr(){}
+
+//
+// Read an event and return a LCCollectionVec of MCParticles
+//
+ IMPL::LCCollectionVec * LCAscHepRdr::readEvent()
+ {
+
+ IMPL::LCCollectionVec * mcVec = nullptr;
+ double c_light = 299.792;// mm/ns
+ //
+ // 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;
+ if( inputFile.eof() )
+ {
+ //
+ // End of File :: ??? Exception ???
+ // -> FG: EOF is not an exception as it happens for every file at the end !
+ //
+ return mcVec;
+ }
+
+ //
+ // Create a Collection Vector
+ //
+ mcVec = new IMPL::LCCollectionVec(LCIO::MCPARTICLE);
+ MCParticleImpl* p;
+ MCParticleImpl* d;
+
+ //
+ // Loop over particles
+ //
+ int ISTHEP; // status code
+ int IDHEP; // PDG code
+ int JMOHEP1; // first mother
+ int JMOHEP2; // last mother
+ int JDAHEP1; // first daughter
+ int JDAHEP2; // last daughter
+ double PHEP1; // px in GeV/c
+ double PHEP2; // py in GeV/c
+ double PHEP3; // pz in GeV/c
+ double PHEP4; // energy in GeV
+ double PHEP5; // mass in GeV/c**2
+ double VHEP1; // x vertex position in mm
+ double VHEP2; // y vertex position in mm
+ double VHEP3; // z vertex position in mm
+ double VHEP4; // production time in mm/c
+
+ std::vector<int> *daughter1 = new std::vector<int> ();
+ std::vector<int> *daughter2 = new std::vector<int> ();
+
+ for( int IHEP=0; IHEP<NHEP; IHEP++ )
+ {
+ //if ( theFileFormat == HEPEvt)
+ if ( false)
+ inputFile >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
+ >> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5;
+ else
+ inputFile >> ISTHEP >> IDHEP
+ >> JMOHEP1 >> JMOHEP2
+ >> JDAHEP1 >> JDAHEP2
+ >> PHEP1 >> PHEP2 >> PHEP3
+ >> PHEP4 >> PHEP5
+ >> VHEP1 >> VHEP2 >> VHEP3
+ >> VHEP4;
+
+ if(inputFile.eof())
+ return nullptr;
+ //
+ // Create a MCParticle and fill it from stdhep info
+ //
+ MCParticleImpl* mcp = new MCParticleImpl();
+ //
+ // PDGID
+ //
+ mcp->setPDG(IDHEP);
+ //
+ // Momentum vector
+ //
+ float p0[3] = {PHEP1,PHEP2,PHEP3};
+ mcp->setMomentum(p0);
+ //
+ // Mass
+ //
+ mcp->setMass(PHEP5);
+ //
+ // Vertex
+ // (missing information in HEPEvt files)
+ double v0[3] = {0.,0.,0.};
+ mcp->setVertex(v0);
+ //
+ // Generator status
+ //
+ mcp->setGeneratorStatus(ISTHEP);
+ //
+ // Simulator status 0 until simulator acts on it
+ //
+ mcp->setSimulatorStatus(0);
+ //
+ // Creation time (note the units)
+ // (No information in HEPEvt files)
+ mcp->setTime(0./c_light);
+ //
+ // Add the particle to the collection vector
+ //
+ mcVec->push_back(mcp);
+ //
+ // Keep daughters information for later
+ //
+ daughter1->push_back(JDAHEP1);
+ daughter2->push_back(JDAHEP2);
+
+ // fg: comment out the mother relationships altogether ....
+
+// //
+// // Add the parent information. The implicit assumption here is that
+// // no particle is read in before its parents.
+// //
+// int fp = _reader->mother1(IHEP) - 1;
+// int lp = _reader->mother2(IHEP) - 1;
+// //
+// // If both first parent and second parent > 0, and second parent >
+// // first parent, assume a range
+// //
+// if( (fp > -1) && (lp > -1) )
+// {
+// if(lp >= fp)
+// {
+// for(int ip=fp;ip<lp+1;ip++)
+// {
+// p = dynamic_cast<MCParticleImpl*>
+// (mcVec->getElementAt(ip));
+// mcp->addParent(p);
+// }
+// }
+// //
+// // If first parent < second parent, assume 2 discreet parents
+// //
+// else
+// {
+// p = dynamic_cast<MCParticleImpl*>
+// (mcVec->getElementAt(fp));
+// mcp->addParent(p);
+// p = dynamic_cast<MCParticleImpl*>
+// (mcVec->getElementAt(lp));
+// mcp->addParent(p);
+// }
+// }
+// //
+// // Only 1 parent > 0, set it
+// //
+// else if(fp > -1)
+// {
+// p = dynamic_cast<MCParticleImpl*>
+// (mcVec->getElementAt(fp));
+// mcp->addParent(p);
+// }
+// else if(lp > -1)
+// {
+// p = dynamic_cast<MCParticleImpl*>
+// (mcVec->getElementAt(lp));
+// mcp->addParent(p);
+// }
+
+ }// End loop over particles
+//
+// Now make a second loop over the particles, checking the daughter
+// information. This is not always consistent with parent
+// information, and this utility assumes all parents listed are
+// parents and all daughters listed are daughters
+//
+ for( int IHEP=0; IHEP<NHEP; IHEP++ )
+ {
+ //
+ // Get the MCParticle
+ //
+ MCParticleImpl* mcp =
+ dynamic_cast<MCParticleImpl*>
+ (mcVec->getElementAt(IHEP));
+ //
+ // Get the daughter information, discarding extra information
+ // sometimes stored in daughter variables.
+ //
+
+ int fd = daughter1->operator[](IHEP) - 1;
+ int ld = daughter2->operator[](IHEP) - 1;
+
+ //
+ // As with the parents, look for range, 2 discreet or 1 discreet
+ // daughter.
+ //
+ if( (fd > -1) && (ld > -1) )
+ {
+ if(ld >= fd)
+ {
+ for(int id=fd;id<ld+1;id++)
+ {
+ //
+ // Get the daughter, and see if it already lists this particle as
+ // a parent.
+ //
+ d = dynamic_cast<MCParticleImpl*>
+ (mcVec->getElementAt(id));
+ int np = d->getParents().size();
+ bool gotit = false;
+ for(int ip=0;ip < np;ip++)
+ {
+ p = dynamic_cast<MCParticleImpl*>
+ (d->getParents()[ip]);
+ if(p == mcp)gotit = true;
+ }
+ //
+ // If not already listed, add this particle as a parent
+ //
+ if(!gotit)d->addParent(mcp);
+ }
+ }
+ //
+ // Same logic, discreet cases
+ //
+ else
+ {
+ d = dynamic_cast<MCParticleImpl*>
+ (mcVec->getElementAt(fd));
+ int np = d->getParents().size();
+ bool gotit = false;
+ for(int ip=0;ip < np;ip++)
+ {
+ p = dynamic_cast<MCParticleImpl*>
+ (d->getParents()[ip]);
+ if(p == mcp)gotit = true;
+ }
+ if(!gotit)d->addParent(mcp);
+ d = dynamic_cast<MCParticleImpl*>
+ (mcVec->getElementAt(ld));
+ np = d->getParents().size();
+ gotit = false;
+ for(int ip=0;ip < np;ip++)
+ {
+ p = dynamic_cast<MCParticleImpl*>
+ (d->getParents()[ip]);
+ if(p == mcp)gotit = true;
+ }
+ if(!gotit)d->addParent(mcp);
+ }
+ }
+ else if(fd > -1)
+ {
+ d = dynamic_cast<MCParticleImpl*>
+ (mcVec->getElementAt(fd));
+ int np = d->getParents().size();
+ bool gotit = false;
+ for(int ip=0;ip < np;ip++)
+ {
+ p = dynamic_cast<MCParticleImpl*>
+ (d->getParents()[ip]);
+ if(p == mcp)gotit = true;
+ }
+ if(!gotit)d->addParent(mcp);
+ }
+ else if(ld > -1)
+ {
+ d = dynamic_cast<MCParticleImpl*>
+ (mcVec->getElementAt(ld));
+ int np = d->getParents().size();
+ bool gotit = false;
+ for(int ip=0;ip < np;ip++)
+ {
+ p = dynamic_cast<MCParticleImpl*>
+ (d->getParents()[ip]);
+ if(p == mcp)gotit = true;
+ }
+ if(!gotit)d->addParent(mcp);
+ }
+ }// End second loop over particles
+ //
+ // Return the collection
+ //
+ return mcVec;
+ }
+
+
+} // namespace UTIL
diff --git a/Generator/src/LCAscHepRdr.h b/Generator/src/LCAscHepRdr.h
new file mode 100644
index 00000000..4d05d62b
--- /dev/null
+++ b/Generator/src/LCAscHepRdr.h
@@ -0,0 +1,40 @@
+#ifndef UTIL_LCAscHepRdr_H
+#define UTIL_LCAscHepRdr_H 1
+
+#include "IMPL/LCCollectionVec.h"
+#include <fstream>
+
+namespace UTIL{
+
+ /**Basic utility for reading a ASCII HEPEvt file and filling
+ * a LCCollectionVec with MCParticles containing the HEPEvt
+ * file information.
+ *
+ * @author Mora de Freitas
+ * @version $Id:
+ */
+ class LCAscHepRdr{
+
+ public:
+
+ /** Open the HEPEvt input file in the constructer
+ */
+ LCAscHepRdr(const char* evfile, int fileFormat) ;
+
+ /** noop
+ */
+ ~LCAscHepRdr() ;
+
+ /** Read an event and return a LCCollectionVec of MCParticles.
+ */
+ IMPL::LCCollectionVec * readEvent() ;
+
+ private:
+ std::ifstream inputFile;
+ int theFileFormat;
+
+ }; // class
+
+} // namespace UTIL
+
+#endif /* ifndef UTIL_LCAscHepRdr_H */
diff --git a/Generator/src/SLCIORdr.cpp b/Generator/src/SLCIORdr.cpp
new file mode 100644
index 00000000..c57a8e7d
--- /dev/null
+++ b/Generator/src/SLCIORdr.cpp
@@ -0,0 +1,169 @@
+#include "SLCIORdr.h"
+#include "GenEvent.h"
+
+#include "lcio.h" //LCIO
+#include "LCIOSTLTypes.h"
+#include "IOIMPL/LCFactory.h"
+#include "EVENT/LCIO.h"
+#include "EVENT/LCEvent.h"
+#include "EVENT/LCCollection.h"
+#include "IO/LCReader.h"
+#include "IMPL/MCParticleImpl.h"
+#include "IMPL/LCCollectionVec.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 <iostream>
+#include <vector>
+#include <fstream>
+
+
+using namespace lcio;
+using namespace IMPL;
+using namespace plcio;
+using namespace std;
+
+
+SLCIORdr::SLCIORdr(string name){
+
+m_slcio_rdr = IOIMPL::LCFactory::getInstance()->createLCReader();
+m_slcio_rdr->open(name.c_str());
+m_processed_event=0;
+}
+
+SLCIORdr::~SLCIORdr(){
+delete m_slcio_rdr;
+}
+
+bool SLCIORdr::mutate(MyHepMC::GenEvent& event){
+
+ IMPL::LCCollectionVec* lcMCVec = new IMPL::LCCollectionVec(LCIO::MCPARTICLE);
+ EVENT::LCEvent *lcEvent = m_slcio_rdr->readNextEvent(LCIO::UPDATE);
+ LCCollection *lcCol = NULL;
+ if(lcEvent) lcCol = lcEvent->getCollection("MCParticle");
+ else return false;
+ if(lcCol){
+ MCParticleImpl* p;
+ //MCParticleImpl* d;
+ int NHEP = lcCol->getNumberOfElements();
+ lcMCVec->parameters() = lcCol->parameters();
+ lcMCVec->resize(NHEP);
+ for( int IHEP=0; IHEP<NHEP; IHEP++ ){
+ MCParticleImpl* mcp = new MCParticleImpl();
+ EVENT::MCParticle* in = dynamic_cast<EVENT::MCParticle*>(lcCol->getElementAt(IHEP));
+ lcMCVec->at(IHEP) = mcp ;
+ mcp->setPDG(in->getPDG());
+ mcp->setCharge(in->getCharge()) ;
+ mcp->setMomentum(in->getMomentum());
+ mcp->setMass(in->getMass());
+ mcp->setVertex(in->getVertex());
+ mcp->setGeneratorStatus(in->getGeneratorStatus());
+ mcp->setSimulatorStatus(in->getSimulatorStatus());
+ mcp->setTime(in->getTime());
+ mcp->setSpin(in->getSpin());
+ }
+ for( int IHEP=0; IHEP<NHEP; IHEP++ ){
+ lcio::MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>(lcMCVec->getElementAt(IHEP));
+ EVENT::MCParticleVec parents = (dynamic_cast<EVENT::MCParticle*>(lcCol->getElementAt(IHEP)))->getParents();
+ int np = parents.size();
+ //cout << "DEBUG: the " << IHEP << "th particle has " << np << " parents ";
+ int i;
+ for(int ip=0;ip<np;ip++){
+ p = dynamic_cast<MCParticleImpl*>(parents[ip]);
+ for(i=0;i<NHEP;i++){
+ if(p==lcCol->getElementAt(i)) break;
+ }
+ if(i==NHEP) cout << "HepLCIOInterfaceNew: error" << endl;
+ mcp->addParent(dynamic_cast<MCParticleImpl*>(lcMCVec->getElementAt(i)));
+ }
+ //ignore daughter table, auto-regonize relationship while addParent()
+ /*
+ EVENT::MCParticleVec daughters = (dynamic_cast<EVENT::MCParticle*>(lcCol->getElementAt(IHEP)))->getDaughters();
+ int nd = daughters.size();
+ //cout << "and " << nd << " daughters." << endl;
+ for(int id=0;id<nd;id++){
+ d = dynamic_cast<MCParticleImpl*>(daughters[id]);
+ for(i=0;i<NHEP;i++){
+ if(d==lcCol->getElementAt(i)) break;
+ }
+ if(i==NHEP) cout << "HepLCIOInterfaceNew: error" << endl;
+ mcp->addDaughter(dynamic_cast<MCParticleImpl*>(lcMCVec->getElementAt(i)));
+ }
+ */
+
+ }
+ }
+ else{
+ cout << "Debug: no MCParticle Collection is read!" << endl;
+ return false;
+ }
+
+
+ m_processed_event ++;
+ int n_mc = lcMCVec->getNumberOfElements();
+ 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*) lcMCVec->getElementAt(i);
+ //std::cout<<"At mc :"<< i <<std::endl;
+ plcio::MCParticle 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;
+
+ mcp.setPDG (mc->getPDG());
+ mcp.setGeneratorStatus (mc->getGeneratorStatus());
+ mcp.setSimulatorStatus (mc->getSimulatorStatus());
+ mcp.setCharge (mc->getCharge());
+ mcp.setTime (mc->getTime());
+ 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.setSpin (mc->getSpin());
+ mcp.setColorFlow (mc->getColorFlow());
+ }
+ // second loop for setting parents and daughters
+
+ for (int i=0; i < n_mc; i++){
+ MCParticleImpl* mc = (MCParticleImpl*) lcMCVec->getElementAt(i);
+ const MCParticleVec & mc_parents = mc->getParents();
+ const MCParticleVec & mc_daughters = mc->getDaughters();
+ plcio::MCParticle 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) ) );
+ }
+ //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) ) );
+ }
+ }
+ event.SetEventHeader( m_processed_event, -99, 9999, "Generator");
+ //std::cout<<"end event :"<< m_processed_event <<std::endl;
+ //for(unsigned int i=0; i< lcMCVec->size(); i++) delete lcMCVec->at(i);
+ delete lcMCVec;
+ return true;
+}
+
+bool SLCIORdr::isEnd(){
+return false;
+}
+
+bool SLCIORdr::configure(){
+return true;
+}
+
+bool SLCIORdr::finish(){
+return true;
+}
diff --git a/Generator/src/SLCIORdr.h b/Generator/src/SLCIORdr.h
new file mode 100644
index 00000000..b429475f
--- /dev/null
+++ b/Generator/src/SLCIORdr.h
@@ -0,0 +1,32 @@
+#ifndef SLCIORdr_h
+#define SLCIORdr_h 1
+
+#include "GenReader.h"
+#include "GenEvent.h"
+
+#include "lcio.h"
+#include "LCIOSTLTypes.h"
+#include "IOIMPL/LCFactory.h"
+#include "EVENT/LCIO.h"
+#include "EVENT/LCEvent.h"
+#include "EVENT/LCCollection.h"
+#include "IO/LCReader.h"
+
+
+class SLCIORdr: public GenReader{
+
+ public:
+ SLCIORdr(string name);
+ ~SLCIORdr();
+ bool configure();
+ bool mutate(MyHepMC::GenEvent& event);
+ bool finish();
+ bool isEnd();
+ private:
+ IO::LCReader* m_slcio_rdr;
+ long m_total_event;
+ long m_processed_event;
+};
+
+#endif
+
diff --git a/Generator/src/StdHepRdr.cpp b/Generator/src/StdHepRdr.cpp
new file mode 100644
index 00000000..ec3bc120
--- /dev/null
+++ b/Generator/src/StdHepRdr.cpp
@@ -0,0 +1,104 @@
+#include "StdHepRdr.h"
+#include "GenEvent.h"
+
+#include "lcio.h" //LCIO
+#include "EVENT/LCIO.h"
+#include "UTIL/LCStdHepRdrNew.h"
+#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 <iostream>
+#include <vector>
+#include <fstream>
+
+
+using namespace lcio;
+using namespace IMPL;
+using namespace plcio;
+using namespace std;
+
+
+StdHepRdr::StdHepRdr(string name){
+
+m_stdhep_rdr = new LCStdHepRdrNew(name.c_str());
+m_stdhep_rdr->printHeader();
+m_total_event = m_stdhep_rdr->getNumberOfEvents() - 1 ;
+m_processed_event=0;
+}
+
+StdHepRdr::~StdHepRdr(){
+delete m_stdhep_rdr;
+}
+
+bool StdHepRdr::mutate(MyHepMC::GenEvent& event){
+ if(isEnd()) return false;
+ LCCollectionVec* mc_vec = m_stdhep_rdr->readEvent();
+ m_processed_event ++;
+ int n_mc = mc_vec->getNumberOfElements();
+ //std::cout<<"Debug: 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();
+ pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i));
+ //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl;
+
+ mcp.setPDG (mc->getPDG());
+ mcp.setGeneratorStatus (mc->getGeneratorStatus());
+ mcp.setSimulatorStatus (mc->getSimulatorStatus());
+ mcp.setCharge (mc->getCharge());
+ mcp.setTime (mc->getTime());
+ 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.setSpin (mc->getSpin());
+ mcp.setColorFlow (mc->getColorFlow());
+ }
+ // second loop for setting parents and daughters
+
+ for (int i=0; i < n_mc; i++){
+ 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);
+ //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) ) );
+ }
+ //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) ) );
+ }
+ }
+
+ event.SetEventHeader( m_processed_event, -99, 9999, "Generator");
+ //std::cout<<"end event :"<< m_processed_event <<std::endl;
+ return true;
+}
+
+bool StdHepRdr::isEnd(){
+if(m_processed_event == m_total_event) {std::cout<<"Have read all events, end now."<<std::endl; return true;}
+else return false;
+}
+
+bool StdHepRdr::configure(){
+return true;
+}
+
+bool StdHepRdr::finish(){
+return true;
+}
diff --git a/Generator/src/StdHepRdr.h b/Generator/src/StdHepRdr.h
new file mode 100644
index 00000000..3d92510b
--- /dev/null
+++ b/Generator/src/StdHepRdr.h
@@ -0,0 +1,28 @@
+#ifndef StdHepRdr_h
+#define StdHepRdr_h 1
+
+#include "GenReader.h"
+#include "GenEvent.h"
+
+#include "lcio.h"
+#include "EVENT/LCIO.h"
+#include "UTIL/LCStdHepRdrNew.h"
+
+
+class StdHepRdr: public GenReader{
+
+ public:
+ StdHepRdr(string name);
+ ~StdHepRdr();
+ bool configure();
+ bool mutate(MyHepMC::GenEvent& event);
+ bool finish();
+ bool isEnd();
+ private:
+ lcio::LCStdHepRdrNew* m_stdhep_rdr;
+ long m_total_event;
+ long m_processed_event;
+};
+
+#endif
+
--
GitLab