Generator/src/GenPrinter.cpp

@@ -3,7 +3,7 @@
 
 DECLARE_COMPONENT(GenPrinter)
 
-bool GenPrinter::mutate(MyHepMC::GenEvent& event){
+bool GenPrinter::mutate(Gen::GenEvent& event){
     auto msglevel = msgLevel();
 
     // only print when current msglevel is MSG::DEBUG/VERBOSE

Generator/src/GenPrinter.h

@@ -5,8 +5,6 @@
 #include "GenEvent.h"
 #include "IGenTool.h"
 
-using namespace std;
-
 class GenPrinter: public extends<AlgTool, IGenTool> {
 public:
     using extends::extends;
@@ -17,7 +15,7 @@ public:
 
 public:
     bool configure_gentool() override;               
-    bool mutate(MyHepMC::GenEvent& event) override;    
+    bool mutate(Gen::GenEvent& event) override;    
     bool finish() override;
 };
 

Generator/src/GenReader.h

 #ifndef GenReader_h
 #define GenReader_h 1
 
+/*
+ * GenReader is a base class for the gentools which load events from file
+ */
+
 #include "GenEvent.h"
 #include "IGenTool.h"
 
-using namespace std;
-
 class GenReader: virtual public IGenTool{
 
     public:
-        ~GenReader();
+        virtual ~GenReader() = 0;
         virtual bool configure_gentool()=0;               
-        virtual bool mutate(MyHepMC::GenEvent& event)=0;    
+        virtual bool mutate(Gen::GenEvent& event)=0;    
         virtual bool finish()=0;
         virtual bool isEnd()=0;
+
+        // The start index is used to skip the first n events.
+        // if start index is 0, then the first event is the first event in the file
+        virtual int startIndex() { return 0; };
 };
 
 #endif

Generator/src/GtBeamBackgroundTool.cpp

@@ -2,8 +2,13 @@
 #include "IBeamBackgroundFileParser.h"
 
 #include "BeamBackgroundFileParserV0.h"
+#include "BeamBackgroundFileParserV1.h"
+#include "BeamBackgroundFileParserV2.h"
 #include "GuineaPigPairsFileParser.h"
 
+#include "CLHEP/Random/RandPoisson.h"
+#include "CLHEP/Random/RandFlat.h"
+
 #include "TVector3.h" // for rotation
 DECLARE_COMPONENT(GtBeamBackgroundTool)
 
@@ -22,9 +27,17 @@ StatusCode GtBeamBackgroundTool::initialize() {
 
         if (format == "BeamBackgroundFileParserV0") {
             init_BeamBackgroundFileParserV0(label, inputfn);
-        } else if (format == "GuineaPigPairsFileParser") {
+        } 
+        else if (format == "GuineaPigPairsFileParser") {
             init_GuineaPigPairsFileParser(label, inputfn);
-        } else {
+        }
+        else if(format == "BeamBackgroundFileParserV1"){
+            init_BeamBackgroundFileParserV1(label, inputfn);
+        }
+        else if(format == "BeamBackgroundFileParserV2"){
+            init_BeamBackgroundFileParserV2(label, inputfn);
+        } 
+        else {
             init_BeamBackgroundFileParserV0(label, inputfn);
         }
         
@@ -44,44 +57,51 @@ StatusCode GtBeamBackgroundTool::finalize() {
 }
 
 
-bool GtBeamBackgroundTool::mutate(MyHepMC::GenEvent& event) {
+bool GtBeamBackgroundTool::mutate(Gen::GenEvent& event) {
     if (m_beaminputs.empty()) {
         return false;
     }
-    // TODO: should sample according to the rates
-    // dummy: get one and stop to generate
-    for (auto& [label, parser]: m_beaminputs) {
-        IBeamBackgroundFileParser::BeamBackgroundData beamdata;
-        auto isok = parser->load(beamdata);
-        if (not isok) {
-            error() << "Failed to load beam background data from the parser " << label << endmsg;
-            return false;
-        }
-        // fill the value
-        float charge = beamdata.pdgid == 11 ? -1: 1;
 
-        TVector3 pos(beamdata.x,beamdata.y,beamdata.z);
-        TVector3 mom(beamdata.px,beamdata.py,beamdata.pz);
+    for (auto& [label, parser]: m_beaminputs) {
 
-        auto itrot = m_rotYmaps.find(label);
-        if (itrot  != m_rotYmaps.end() ) {
-            info() << "Apply rotation along Y " << itrot->second << endmsg;
+        int nparticles = 0;
+        int startNo = 0;
+        parser->SampleParticleNum(nparticles, startNo);
+        for(auto ipart = 0; ipart < nparticles; ++ipart){
+            IBeamBackgroundFileParser::BeamBackgroundData beamdata;
+            auto isok = parser->load(beamdata, ipart+startNo);
+            if (not isok) {
+                error() << "Failed to load beam background data from the parser " << label << endmsg;
+                return false;
+            }
+
+            float charge = beamdata.charge;
+
+            TVector3 pos(beamdata.x,beamdata.y,beamdata.z);
+            TVector3 mom(beamdata.px,beamdata.py,beamdata.pz);
+
+
+            auto itrot = m_rotYMap.find(label);
+            if (itrot  != m_rotYMap.end() ) {
+                //info() << "Apply rotation along Y " << itrot->second << endmsg;
+
+                pos.RotateY(itrot->second);
+                mom.RotateY(itrot->second);
+            }            
+
+            // create the MC particle
+            auto mcp = event.m_mc_vec.create();
+            mcp.setPDG(beamdata.pdgid);
+            mcp.setGeneratorStatus(1);
+            mcp.setSimulatorStatus(1);
+            mcp.setCharge(static_cast<float>(charge));
+            mcp.setTime(beamdata.t);
+            mcp.setMass(beamdata.mass);
+            mcp.setVertex(edm4hep::Vector3d(pos.X(), pos.Y(), pos.Z()));
+            mcp.setMomentum(edm4hep::Vector3f(mom.X(), mom.Y(), mom.Z()));
 
-            pos.RotateY(itrot->second);
-            mom.RotateY(itrot->second);
         }
 
-        // create the MC particle
-        auto mcp = event.m_mc_vec.create();
-        mcp.setPDG(beamdata.pdgid);
-        mcp.setGeneratorStatus(1);
-        mcp.setSimulatorStatus(1);
-        mcp.setCharge(static_cast<float>(charge));
-        mcp.setTime(beamdata.t);
-        mcp.setMass(beamdata.mass);
-        mcp.setVertex(edm4hep::Vector3d(pos.X(), pos.Y(), pos.Z())); 
-        mcp.setMomentum(edm4hep::Vector3f(mom.X(), mom.Y(), mom.Z()));
-
     }
     
     return true;
@@ -98,17 +118,49 @@ bool GtBeamBackgroundTool::configure_gentool() {
 
 bool GtBeamBackgroundTool::init_BeamBackgroundFileParserV0(const std::string& label,
                                                            const std::string& inputfn) {
-    double beamE = 120.;
-    auto itBeamE = m_Ebeammaps.find(label);
-    if (itBeamE != m_Ebeammaps.end()) {
-        beamE = itBeamE->second;
-    }
     info() << "Initializing beam background ... "
            << label << " "
-           << beamE << " "
+           << m_Ebeam << " "
            << inputfn
            << endmsg;
-    m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV0>(inputfn, 11, beamE);
+    m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV0>(inputfn, 11, m_Ebeam);
+
+    return true;
+}
+
+bool GtBeamBackgroundTool::init_BeamBackgroundFileParserV1(const std::string& label,
+                                                           const std::string& inputfn) {
+    info() << "Initializing beam background ... " << label << inputfn <<endmsg;
+    auto itRate = m_ratemaps.find(label);
+    if(itRate == m_ratemaps.end()){ 
+      error() << "Did not find the beam process rate for: " << label << endmsg;
+      return false;
+    }
+    auto itNmcp = m_Nmcpmaps.find(label);
+    if(itNmcp == m_Nmcpmaps.end()){ 
+      error() << "Did not find the beam process Nmcp for: " << label << endmsg;
+      return false;
+    }
+    m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV1>(inputfn, "BeamTree", m_Ebeam, itRate->second, m_timewindow, itNmcp->second);
+
+    return true;
+}
+
+bool GtBeamBackgroundTool::init_BeamBackgroundFileParserV2(const std::string& label,
+                                                           const std::string& inputfn) {
+
+    info() << "Initializing beam background ... " << label << inputfn <<endmsg;
+    auto itTime = m_timebkgmaps.find(label);
+    if(itTime == m_timebkgmaps.end()){ 
+      error() << "Did not find the beam process time for: " << label << endmsg;
+      return false;
+    }
+    auto itNmcp = m_Nmcpmaps.find(label);
+    if(itNmcp == m_Nmcpmaps.end()){ 
+      error() << "Did not find the beam process Nmcp for: " << label << endmsg;
+      return false;
+    }
+    m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV2>(inputfn, "BeamTree", m_Ebeam, itTime->second, itNmcp->second);
 
     return true;
 }

Generator/src/GtBeamBackgroundTool.h

@@ -42,25 +42,36 @@ public:
     StatusCode finalize() override;
 
     // IGenTool
-    bool mutate(MyHepMC::GenEvent& event) override;
+    bool mutate(Gen::GenEvent& event) override;
     bool finish() override;
     bool configure_gentool() override;
 
 
 private:
     bool init_BeamBackgroundFileParserV0(const std::string& label, const std::string& inputfn);
+    bool init_BeamBackgroundFileParserV1(const std::string& label, const std::string& inputfn);
+    bool init_BeamBackgroundFileParserV2(const std::string& label, const std::string& inputfn);
     bool init_GuineaPigPairsFileParser(const std::string& label, const std::string& inputfn);
 
 private:
     Gaudi::Property<std::map<std::string, std::string>> m_inputmaps{this, "InputFileMap"};
     Gaudi::Property<std::map<std::string, std::string>> m_formatmaps{this, "InputFormatMap"};
-    Gaudi::Property<std::map<std::string, double>>      m_ratemaps {this, "InputRateMap"};
+    Gaudi::Property<std::map<std::string, double>>      m_ratemaps {this, "InputRateMap"}; // unit: Hz
+
+    // Detector time window. Should be different for different sub-Det. Unit: s
+    Gaudi::Property<double>      m_timewindow{this, "TimeWindow"}; 
 
     // unit of beam energy: GeV
-    Gaudi::Property<std::map<std::string, double>>      m_Ebeammaps{this, "InputBeamEnergyMap"};
+    Gaudi::Property<double>      m_Ebeam{this, "InputBeamEnergy"};
+
+    // Rotation along Y for single beam background. Unit: rad
+    Gaudi::Property<std::map<std::string, double>>      m_rotYMap {this, "RotationAlongYMap"};
+
+    // Time of bunch crossing. Unit: ns, consistent with simulation
+    Gaudi::Property<std::map<std::string, double>>      m_timebkgmaps{this, "TimeBkgMap"};
 
-    // unit of the rotation along Y: rad
-    Gaudi::Property<std::map<std::string, double>>      m_rotYmaps {this, "RotationAlongYMap"};
+    // Number of McParticles in different beambkg. -1: pair use one file and single beam use rate*time; >=0: fixed number
+    Gaudi::Property<std::map<std::string, int>>      m_Nmcpmaps{this, "NumberMcParticle"};
 
 private:
     std::map<std::string, std::shared_ptr<IBeamBackgroundFileParser>> m_beaminputs;

Generator/src/GtGunTool.cpp

@@ -73,17 +73,6 @@ GtGunTool::initialize() {
     }
     
     // others should be empty or specify
-    if (m_thetamins.value().size()
-        && m_thetamins.value().size() != m_particles.value().size()) {
-        error() << "Mismatched thetamins and particles." << endmsg;
-        return StatusCode::FAILURE;
-    }
-    if (m_thetamaxs.value().size()
-        && m_thetamaxs.value().size() != m_particles.value().size()) {
-        error() << "Mismatched thetamaxs and particles." << endmsg;
-        return StatusCode::FAILURE;
-    }
-
     if (m_phimins.value().size()
         && m_phimins.value().size() != m_particles.value().size()) {
         error() << "Mismatched phimins and particles." << endmsg;
@@ -95,6 +84,47 @@ GtGunTool::initialize() {
         return StatusCode::FAILURE;
     }
 
+    if (m_usecostheta.value()) {
+        if (m_costhetamins.value().size() != m_particles.value().size()) {
+            error() << "Mismatched CosthetaMins and particles." << endmsg;
+            return StatusCode::FAILURE;
+        }
+        if (m_costhetamaxs.value().size() != m_particles.value().size()) {
+            error() << "Mismatched CosthetaMaxs and particles." << endmsg;
+            return StatusCode::FAILURE;
+        }
+        for (int i=0; i<m_thetamins.value().size(); ++i) {
+            if ((m_thetamins.value()[i] < -1) || (m_thetamins.value()[i] > 1)) {
+                error() << "UseCostheta: ThetaMins has values outside the range [-1, 1]." << endmsg;
+                return StatusCode::FAILURE;
+            }
+            if ((m_thetamaxs.value()[i] < -1) || (m_thetamaxs.value()[i] > 1)) {
+                error() << "UseCostheta: ThetaMaxs has values outside the range [-1, 1]." << endmsg;
+                return StatusCode::FAILURE;
+            }
+        }
+    } else {
+        if (m_thetamins.value().size()
+            && m_thetamins.value().size() != m_particles.value().size()) {
+            error() << "Mismatched thetamins and particles." << endmsg;
+            return StatusCode::FAILURE;
+        }
+        if (m_thetamaxs.value().size()
+            && m_thetamaxs.value().size() != m_particles.value().size()) {
+            error() << "Mismatched thetamaxs and particles." << endmsg;
+            return StatusCode::FAILURE;
+        }
+    }
+
+    // Time
+    if (m_times.value().size()==0){
+      for(int i=0; i<m_particles.value().size(); i++) m_times.value().push_back(0);
+    }
+    else if (m_times.value().size() != m_particles.value().size()) {
+        error() << "Mismatched times and particles." << endmsg;
+        return StatusCode::FAILURE;
+    }
+
     return sc;
 }
 
@@ -105,7 +135,7 @@ GtGunTool::finalize() {
 }
 
 bool
-GtGunTool::mutate(MyHepMC::GenEvent& event) {
+GtGunTool::mutate(Gen::GenEvent& event) {
 
     TDatabasePDG* db_pdg = TDatabasePDG::Instance();
 
@@ -163,7 +193,9 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
         mcp.setGeneratorStatus(1);
         mcp.setSimulatorStatus(1);
         mcp.setCharge(static_cast<float>(charge));
-        mcp.setTime(0.0);
+        mcp.setTime(m_times.value()[i]);
+        //if (i<m_times.value().size()) { mcp.setTime(m_times.value()[i]); }
+        //else { mcp.setTime(0.0); }
         mcp.setMass(mass);
 
         // Unit is mm
@@ -232,16 +264,26 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
             return false;
         }
 
+        double costheta = 0;
+        double theta = 0;
+        if (m_usecostheta.value()) {
+            costheta = m_costhetamins.value()[i]==m_costhetamaxs.value()[i] ? m_costhetamins.value()[i] : CLHEP::RandFlat::shoot(m_costhetamins.value()[i], m_costhetamaxs.value()[i]);
+        } else {
+            theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]);
+            costheta = cos(theta*acos(-1)/180);
+        }
+        double phi =   m_phimins.value()[i]==m_phimaxs.value()[i] ? m_phimins.value()[i] : CLHEP::RandFlat::shoot(m_phimins.value()[i], m_phimaxs.value()[i]);
 
-        double theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]);
-        double phi =   m_phimins  .value()[i]==m_phimaxs  .value()[i] ? m_phimins  .value()[i] : CLHEP::RandFlat::shoot(m_phimins  .value()[i], m_phimaxs  .value()[i]);
-        double costheta = cos(theta*acos(-1)/180);
         double phi_  = phi*acos(-1)/180;
         double sintheta = sqrt(1.-costheta*costheta);
         double px = p*sintheta*cos(phi_);
         double py = p*sintheta*sin(phi_);
         double pz = p*costheta;
-        std::cout<<"GenGt p="<<p<<", px="<<px<<",py="<<py<<",pz="<<pz<<",theta="<<theta<<",phi="<<phi<<std::endl;
+        if(m_usecostheta.value()) {
+            std::cout<<"GenGt p="<<p<<", px="<<px<<",py="<<py<<",pz="<<pz<<",costheta="<<costheta<<",phi="<<phi<<std::endl;
+        } else {
+            std::cout<<"GenGt p="<<p<<", px="<<px<<",py="<<py<<",pz="<<pz<<",theta="<<theta<<",phi="<<phi<<std::endl;
+        }
         mcp.setMomentum(edm4hep::Vector3f(px,py,pz));
         // mcp.setMomentumAtEndpoint();
         // mcp.setSpin();

Generator/src/GtGunTool.h

@@ -26,7 +26,7 @@ public:
     StatusCode finalize() override;
 
     // IGenTool
-    bool mutate(MyHepMC::GenEvent& event) override;
+    bool mutate(Gen::GenEvent& event) override;
     bool finish() override;
     bool configure_gentool() override;
 
@@ -57,6 +57,12 @@ private:
     Gaudi::Property<std::vector<double>> m_phimins{this, "PhiMins"};
     Gaudi::Property<std::vector<double>> m_phimaxs{this, "PhiMaxs"};
 
+    Gaudi::Property<bool> m_usecostheta{this, "UseCostheta", false};
+    Gaudi::Property<std::vector<double>> m_costhetamins{this, "CosthetaMins"};
+    Gaudi::Property<std::vector<double>> m_costhetamaxs{this, "CosthetaMaxs"};
+    // For time
+    Gaudi::Property<std::vector<double>> m_times{this, "Times"};
+
 };
 
 

Generator/src/GtPythiaTool.cpp

+#include "GtPythiaTool.hh"
+
+DECLARE_COMPONENT(GtPythiaTool)
+
+StatusCode GtPythiaTool::initialize() {
+    StatusCode sc;
+    m_pythia = std::make_unique<Pythia8::Pythia>();
+    // configure with the card file first
+    m_pythia->readFile(m_card.value());
+    // tune with the additional commands
+    for (auto cmd: m_cmds.value()) {
+        m_pythia->readString(cmd);
+    }
+    // initialize pythia
+    m_pythia->init();
+    return sc;
+}
+
+StatusCode GtPythiaTool::finalize() {
+    StatusCode sc;
+    return sc;
+}
+
+bool GtPythiaTool::mutate(Gen::GenEvent& event) {
+    // generate the event
+    while(!m_pythia->next()) {
+        // if failed, try again
+    }
+
+    // get the particles
+    auto& pythia_particles = m_pythia->event;
+    // loop over the particles
+    for (int i = 0; i < pythia_particles.size(); ++i) {
+        auto& p = pythia_particles[i];
+        // create the MCParticle
+        auto mcp = event.getMCVec().create();
+        // set the properties
+        mcp.setPDG(p.id());
+        int status = 0;
+        if (p.isFinal()) {
+            status = 1;
+        } else {
+            status = 0;
+        }
+        mcp.setGeneratorStatus(status);
+        mcp.setCharge(p.charge());
+        mcp.setTime(p.tProd());
+        mcp.setMass(p.m());
+        mcp.setVertex(edm4hep::Vector3d(p.xProd(), p.yProd(), p.zProd()));
+        // update the endpoint later
+        mcp.setEndpoint(edm4hep::Vector3d(p.xProd(), p.yProd(), p.zProd()));
+        mcp.setMomentum(edm4hep::Vector3f(p.px(), p.py(), p.pz()));
+    }
+    // setup the relationships (mother and daughter)
+    for (int i = 0; i < pythia_particles.size(); ++i) {
+        auto& p = pythia_particles[i];
+        auto mcp = event.getMCVec()[i];
+
+        auto mother_list = p.motherList();
+        for (auto idx: mother_list) {
+            auto mother = event.getMCVec()[idx];
+            mcp.addToParents(mother);
+        }
+
+        auto daughter_list = p.daughterList();
+        bool need_endpoint = true;
+        for (auto idx: daughter_list) {
+            auto daughter = event.getMCVec()[idx];
+            mcp.addToDaughters(daughter);
+            // Update the endpoint to the daughter's vertex
+            if (need_endpoint) {
+                mcp.setEndpoint(daughter.getVertex());
+                need_endpoint = false;
+            }
+        }
+    }
+    return true;
+}
+
+bool GtPythiaTool::finish() {
+    return true;
+}
+
+bool GtPythiaTool::configure_gentool() {
+    return true;
+}
\ No newline at end of file

Generator/src/GtPythiaTool.hh

+#ifndef GtPythiaTool_hh
+#define GtPythiaTool_hh
+
+/*
+ * Description:
+ *   This tool is used to generate particles using Pythia.
+ *   User need to specify the card.
+ */
+
+#include <GaudiKernel/AlgTool.h>
+#include <Gaudi/Property.h>
+#include "IGenTool.h"
+
+#include "Pythia8/Pythia.h"
+
+#include <vector>
+#include <memory>
+
+class GtPythiaTool: public extends<AlgTool, IGenTool> {
+public:
+    using extends::extends;
+
+    // Overriding initialize and finalize
+    StatusCode initialize() override;
+    StatusCode finalize() override;
+
+    // IGenTool
+    bool mutate(Gen::GenEvent& event) override;
+    bool finish() override;
+    bool configure_gentool() override;
+
+private:
+    std::unique_ptr<Pythia8::Pythia> m_pythia;
+
+    // below properties will be used by Pythia
+    Gaudi::Property<std::vector<std::string>> m_cmds{this, "Commands"};
+    Gaudi::Property<std::string> m_card{this, "Card"};
+};
+
+#endif
\ No newline at end of file

Generator/src/GuineaPigPairsFileParser.h

@@ -24,6 +24,8 @@ public:
     GuineaPigPairsFileParser(const std::string& filename);
 
     bool load(IBeamBackgroundFileParser::BeamBackgroundData&);
+    bool load(IBeamBackgroundFileParser::BeamBackgroundData&, int) { return 0; }
+    bool SampleParticleNum(int&, int&) { return true; }
 
 private:
     std::ifstream m_input;

Generator/src/HepMCRdr.cpp

@@ -20,19 +20,19 @@
 
 
 using namespace edm4hep;
-using namespace std;
+//using namespace std;
 
 DECLARE_COMPONENT(HepMCRdr)
 
 HepMCRdr::~HepMCRdr(){
-delete ascii_in;
+    delete ascii_in;
 }
 
-bool HepMCRdr::mutate(MyHepMC::GenEvent& event){
+bool HepMCRdr::mutate(Gen::GenEvent& event){
 
+    ++m_processed_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;
@@ -93,13 +93,13 @@ bool HepMCRdr::mutate(MyHepMC::GenEvent& event){
 }
 
 bool HepMCRdr::isEnd(){
-return false;
+    return false;
 }
 
 bool HepMCRdr::configure_gentool(){
     ascii_in = new HepMC::IO_GenEvent(m_filename.value().c_str(),std::ios::in);
 
-    m_processed_event=0;
+    m_processed_event=-1;
     return true;
 }
 
@@ -115,6 +115,17 @@ HepMCRdr::initialize() {
         return StatusCode::FAILURE;
     }
 
+    // skip the first n events if startIndex is not 0.
+    if (startIndex() > 0) {
+        for (int i=0; i<startIndex(); ++i) {
+            ++m_processed_event;
+            HepMC::GenEvent* evt = ascii_in->read_next_event();
+            if(!evt) break;
+            delete evt;
+        }
+        info() << "Skip the first " << startIndex() << " events." << endmsg;
+    }
+
     return sc;
 }
 
@@ -127,3 +138,7 @@ HepMCRdr::finalize() {
     }
     return sc;
 }
+
+int HepMCRdr::startIndex(){
+    return m_startIndex.value();
+}
\ No newline at end of file

Generator/src/HepMCRdr.h

@@ -20,9 +20,11 @@ class HepMCRdr: public extends<AlgTool, GenReader> {
         StatusCode finalize() override;    
 
         bool configure_gentool();               
-        bool mutate(MyHepMC::GenEvent& event);    
+        bool mutate(Gen::GenEvent& event);    
         bool finish();
         bool isEnd();
+
+        int startIndex() override;
     private:
         HepMC::IO_GenEvent *ascii_in{nullptr};
         long m_total_event{-1};
@@ -30,6 +32,7 @@ class HepMCRdr: 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"};
 
 };
 

Generator/src/HepevtRdr.cpp

@@ -7,12 +7,10 @@
 #include "IMPL/MCParticleImpl.h"
 
 
-#include "plcio/MCParticle.h" //plcio
-#include "plcio/MCParticleObj.h"
-#include "plcio/MCParticleCollection.h"
-#include "plcio/DoubleThree.h"
-#include "plcio/FloatThree.h"
-#include "plcio/EventHeaderCollection.h"
+#include "edm4hep/MCParticle.h" //edm4hep
+#include "edm4hep/MCParticleObj.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/EventHeaderCollection.h"
 
 
 
@@ -23,8 +21,7 @@
 
 using namespace lcio;
 using namespace IMPL;
-using namespace plcio;
-using namespace std;
+using namespace edm4hep;
 
 typedef enum HEPFILEFORMATS
 {
@@ -35,28 +32,70 @@ typedef enum HEPFILEFORMATS
 }  HEPFILEFORMAT;
 
 
-HepevtRdr::HepevtRdr(string name){
+DECLARE_COMPONENT(HepevtRdr)
 
-m_hepevt_rdr = new UTIL::LCAscHepRdr(name.c_str(), hepevt);
-m_processed_event=0;
-std::cout<<"initial hepevt_rdr"<<std::endl;
+HepevtRdr::~HepevtRdr(){
+    delete m_hepevt_rdr;
 }
 
-HepevtRdr::~HepevtRdr(){
-delete m_hepevt_rdr;
+StatusCode HepevtRdr::initialize() {
+    StatusCode sc;
+    if (not configure_gentool()) {
+        error() << "failed to initialize." << endmsg;
+        return StatusCode::FAILURE;
+    }
+
+    // skip the first n events if startIndex is not 0.
+    if (startIndex() > 0) {
+        for (int i=0; i<startIndex(); ++i) {
+            ++m_processed_event;
+            LCCollectionVec* mc_vec = m_hepevt_rdr->readEvent();
+            if(mc_vec==nullptr) {
+                break;
+            }
+            delete mc_vec;
+        }
+        info() << "Skip the first " << startIndex() << " events." << endmsg;
+    }
+
+    return sc;
 }
 
-bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
+StatusCode HepevtRdr::finalize() {
+    StatusCode sc;
+    if (not finish()) {
+        error() << "Failed to finalize." << endmsg;
+        return StatusCode::FAILURE;
+    }
+    return sc;
+}
+
+
+bool HepevtRdr::configure_gentool(){
+    int format = hepevt;
+    if (m_format == "HEPEvt") {
+        format = HEPEvt;
+    } else if (m_format == "hepevt") {
+        format = hepevt;
+    }
+
+    m_hepevt_rdr = new UTIL::LCAscHepRdr(m_filename.value().c_str(), format);
+    m_processed_event=-1;
+    std::cout<<"initial hepevt_rdr"<<std::endl;
+    return true;
+}
+
+bool HepevtRdr::mutate(Gen::GenEvent& event){
+    ++m_processed_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;
+    int n_mc = mc_vec->size();
+    info()<<"Read event :"<< m_processed_event <<", mc size :"<< n_mc <<endmsg;
     std::map<int, int> pmcid_lmcid;
     for (int i=0; i < n_mc; i++){
         MCParticleImpl* mc = (MCParticleImpl*) mc_vec->getElementAt(i);
-        //std::cout<<"At mc :"<< i <<std::endl;
-        plcio::MCParticle mcp = event.m_mc_vec.create();
+        // std::cout<<"At mc :"<< i <<std::endl;
+        auto mcp = event.m_mc_vec.create();
         pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i));
         //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl;
                                  
@@ -68,8 +107,8 @@ bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
         mcp.setMass               (mc->getMass());
         mcp.setVertex             (mc->getVertex()); 
         mcp.setEndpoint           (mc->getEndpoint());
-        mcp.setMomentum           (FloatThree(float(mc->getMomentum()[0]), float(mc->getMomentum()[1]), float(mc->getMomentum()[2]) ));
-        mcp.setMomentumAtEndpoint (FloatThree(float(mc->getMomentumAtEndpoint()[0]), float(mc->getMomentumAtEndpoint()[1]), float(mc->getMomentumAtEndpoint()[2]) ));
+        mcp.setMomentum           (edm4hep::Vector3f(float(mc->getMomentum()[0]), float(mc->getMomentum()[1]), float(mc->getMomentum()[2]) ));
+        mcp.setMomentumAtEndpoint (edm4hep::Vector3f(float(mc->getMomentumAtEndpoint()[0]), float(mc->getMomentumAtEndpoint()[1]), float(mc->getMomentumAtEndpoint()[2]) ));
         mcp.setSpin               (mc->getSpin());
         mcp.setColorFlow          (mc->getColorFlow());
     }
@@ -79,32 +118,35 @@ bool HepevtRdr::mutate(MyHepMC::GenEvent& event){
         MCParticleImpl* mc = (MCParticleImpl*) mc_vec->getElementAt(i);
         const MCParticleVec & mc_parents = mc->getParents();
         const MCParticleVec & mc_daughters = mc->getDaughters();
-        plcio::MCParticle pmc = event.m_mc_vec.at(i);
+        auto pmc = event.m_mc_vec.at(i);
         //std::cout<<"mc at "<< i<<", parent size "<<mc_parents.size() <<std::endl;
-        for(unsigned int j=0; j< mc_parents.size(); j++){int p_id = mc_parents.at(j)->id();
-                                                 //std::cout<<"parent id "<<p_id<<std::endl;
-                                                 pmc.addParent( event.m_mc_vec.at( pmcid_lmcid.at(p_id) ) );
-                                                }
+        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.addDaughter( event.m_mc_vec.at( pmcid_lmcid.at(d_id) ) );
-                                                }
+        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;
+    delete mc_vec;
     return true;
 }
 
 bool HepevtRdr::isEnd(){
-return false;
-}
-
-bool HepevtRdr::configure(){
-return true;
+    return false;
 }
 
 bool HepevtRdr::finish(){
-return true;
+    return true;
 }
+
+int HepevtRdr::startIndex(){
+    return m_startIndex.value();
+}
\ No newline at end of file