diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..80e594d92eed792026e0bdeaf384522832ff6a76 --- /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 0000000000000000000000000000000000000000..001dd825a84c3f6ec04c3d8a61cfeeabaa73fb36 --- /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 ef5f37288904999d8fff0953e54169bc519b7ca7..a9a5cd9f9e5946684d60058219d2b28725528c08 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 0000000000000000000000000000000000000000..0d8c4698a1fe5a619d41bd88b084adeaafbc0c27 --- /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 0000000000000000000000000000000000000000..43b4339b71eed42a04dc862bc0cf539956b26b60 --- /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 0000000000000000000000000000000000000000..1cd5ccb3c6bcc02df464106e391b846a77bbd159 --- /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 0000000000000000000000000000000000000000..dae1b3772d3044ad3434898ec1c216ab41493d57 --- /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 0000000000000000000000000000000000000000..827f5b4055dec886e91169b5407ba2ddda026e98 --- /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 0000000000000000000000000000000000000000..9fe8a1c9d9fa3306db67471fc686885fbbfa2c8b --- /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 0000000000000000000000000000000000000000..04c38c534fcffd6283b14842824df2450b867a91 --- /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 0000000000000000000000000000000000000000..b9abac0a8a6c994872e73d094a07a65c1a99d87b --- /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 0000000000000000000000000000000000000000..f1f93a571dabc98f710c660754ae6507248c0212 --- /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 0000000000000000000000000000000000000000..30e6db4365343c44f0f7de97781774f0e5700db3 --- /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 0000000000000000000000000000000000000000..22acd3ff013603ace42a4086c61d2621fe5e2fdf --- /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 0000000000000000000000000000000000000000..61c777bba02c8f8a807940162322691c04d127ae --- /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 0000000000000000000000000000000000000000..7f036e2e970d5af1292d0ac22bcfc8769748862c --- /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 0000000000000000000000000000000000000000..9ee2d7fd411f554fdc2114ca854f2fc52f8e7a62 --- /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 0000000000000000000000000000000000000000..5562baf5de1088c292ce37a4a48adbdb6226ef61 --- /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 0000000000000000000000000000000000000000..87af6bd8167d7f1d57415b537de53923c136a34a --- /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 0000000000000000000000000000000000000000..eabca0bb170238264887b321795fcc401b36f777 --- /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 0000000000000000000000000000000000000000..4d05d62b07cb1ec7959622f2e2c49b5eb3c6dee0 --- /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 0000000000000000000000000000000000000000..c57a8e7ded7b0ba100c4d26cc4725537c6f2e78d --- /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 0000000000000000000000000000000000000000..b429475fee3f53e80c48b881a54fc1bbea324077 --- /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 0000000000000000000000000000000000000000..ec3bc12096ed9f5b718e8adec6d74386e91759a0 --- /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 0000000000000000000000000000000000000000..3d92510bb90e9abe2dbc60f29b2770b563d35598 --- /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 +