diff --git a/Generator/src/SLCIORdr.cpp b/Generator/src/SLCIORdr.cpp index 14dcec3af8c57a27dede1e28f92ec5d1a279a159..7478f3408a9f26fbdb940d66e333bc490590143f 100644 --- a/Generator/src/SLCIORdr.cpp +++ b/Generator/src/SLCIORdr.cpp @@ -27,162 +27,177 @@ using namespace lcio; using namespace IMPL; using namespace edm4hep; -using namespace std; DECLARE_COMPONENT(SLCIORdr) SLCIORdr::~SLCIORdr(){ -delete m_slcio_rdr; + delete m_slcio_rdr; } bool SLCIORdr::mutate(MyHepMC::GenEvent& event){ + ++m_processed_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 = nullptr; + //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()); + } - 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 << "Heedm4hepInterfaceNew: 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 << "Heedm4hepInterfaceNew: error" << endl; - mcp->addToDaughters(dynamic_cast<MCParticleImpl*>(lcMCVec->getElementAt(i))); - } - */ - - } + 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) error() << "SLCIORdr: error" << endmsg; + mcp->addParent(dynamic_cast<MCParticleImpl*>(lcMCVec->getElementAt(i))); } - else{ - cout << "Debug: no MCParticle Collection is read!" << endl; - return false; + //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 << "Heedm4hepInterfaceNew: error" << endl; + mcp->addToDaughters(dynamic_cast<MCParticleImpl*>(lcMCVec->getElementAt(i))); } + */ + + } + } else { + error() << "Debug: no MCParticle Collection is read!" << endmsg; + 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; - auto mcp = event.m_mc_vec.create(); - pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i)); - //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl; - - 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 (Vector3f(float(mc->getMomentum()[0]), float(mc->getMomentum()[1]), float(mc->getMomentum()[2]) )); - mcp.setMomentumAtEndpoint (Vector3f(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(); - auto pmc = event.m_mc_vec.at(i); - //std::cout<<"mc at "<< i<<", parent size "<<mc_parents.size() <<std::endl; - for(unsigned int j=0; j< mc_parents.size(); j++){int p_id = mc_parents.at(j)->id(); - //std::cout<<"parent id "<<p_id<<std::endl; - pmc.addToParents( event.m_mc_vec.at( pmcid_lmcid.at(p_id) ) ); - } - //std::cout<<"mc at "<< i<<", daughter size "<<mc_daughters.size() <<std::endl; - for(unsigned int j=0; j< mc_daughters.size(); j++){int d_id = mc_daughters.at(j)->id(); - //std::cout<<"daughter id "<<d_id<<std::endl; - pmc.addToDaughters( 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; + 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; + auto mcp = event.m_mc_vec.create(); + pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i)); + //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl; + + 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 (Vector3f(float(mc->getMomentum()[0]), float(mc->getMomentum()[1]), float(mc->getMomentum()[2]) )); + mcp.setMomentumAtEndpoint (Vector3f(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(); + auto pmc = event.m_mc_vec.at(i); + //std::cout<<"mc at "<< i<<", parent size "<<mc_parents.size() <<std::endl; + for(unsigned int j=0; j< mc_parents.size(); j++){ + int p_id = mc_parents.at(j)->id(); + //std::cout<<"parent id "<<p_id<<std::endl; + pmc.addToParents( event.m_mc_vec.at( pmcid_lmcid.at(p_id) ) ); + } + //std::cout<<"mc at "<< i<<", daughter size "<<mc_daughters.size() <<std::endl; + for(unsigned int j=0; j< mc_daughters.size(); j++){ + int d_id = mc_daughters.at(j)->id(); + //std::cout<<"daughter id "<<d_id<<std::endl; + pmc.addToDaughters( 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; + return false; } bool SLCIORdr::configure_gentool(){ - m_slcio_rdr = IOIMPL::LCFactory::getInstance()->createLCReader(); - m_slcio_rdr->open(m_filename.value().c_str()); - m_processed_event=0; + m_slcio_rdr = IOIMPL::LCFactory::getInstance()->createLCReader(); + m_slcio_rdr->open(m_filename.value().c_str()); + m_processed_event=-1; - return true; + return true; } bool SLCIORdr::finish(){ -return true; + return true; } StatusCode SLCIORdr::initialize() { - StatusCode sc; - if (not configure_gentool()) { - error() << "failed to initialize." << endmsg; - return StatusCode::FAILURE; - } - - return sc; + StatusCode sc; + if (not configure_gentool()) { + error() << "failed to initialize." << endmsg; + return StatusCode::FAILURE; + } + + // skip the first n events if startIndex is not 0. + if (startIndex() > 0) { + m_slcio_rdr->skipNEvents(startIndex()); + info() << "Skip the first " << startIndex() << " events." << endmsg; + } + + return sc; } StatusCode SLCIORdr::finalize() { - StatusCode sc; - if (not finish()) { - error() << "Failed to finalize." << endmsg; - return StatusCode::FAILURE; - } - return sc; + StatusCode sc; + if (not finish()) { + error() << "Failed to finalize." << endmsg; + return StatusCode::FAILURE; + } + return sc; } + +int SLCIORdr::startIndex() { + return m_startIndex.value(); +} \ No newline at end of file diff --git a/Generator/src/SLCIORdr.h b/Generator/src/SLCIORdr.h index a37fa5edb610093a984f4a907f308d5f1fad50bc..7640d0780b54dd0612444e494e723a5972afe8bf 100644 --- a/Generator/src/SLCIORdr.h +++ b/Generator/src/SLCIORdr.h @@ -30,6 +30,8 @@ class SLCIORdr: public extends<AlgTool, GenReader> { bool mutate(MyHepMC::GenEvent& event) override; bool finish() override; bool isEnd() override; + + int startIndex() override; private: IO::LCReader* m_slcio_rdr{nullptr}; long m_total_event{-1}; @@ -37,6 +39,7 @@ class SLCIORdr: public extends<AlgTool, GenReader> { // input file name Gaudi::Property<std::string> m_filename{this, "Input"}; + Gaudi::Property<int> m_startIndex{this, "StartIndex", 0, "Default start index"}; };