diff --git a/CMakeLists.txt b/CMakeLists.txt
index b5aff3852bae9599d85aa429ad9da92400449c9f..8f17382b31c4fff3f93b48257d822be58f31a2b5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,10 +15,10 @@ ENDIF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
#
#---Options-------------------------------------------------------------------------
option(DD4HEP_USE_XERCESC "Enable 'Detector Builders' based on XercesC" OFF)
-option(DD4HEP_USE_PYROOT "Enable 'Detector Builders' based on PyROOT" OFF) # does not work (compile error)
-option(DD4HEP_USE_GEANT4 "Enable the simulation part based on Geant4" OFF)
-option(DD4HEP_USE_GEAR "Build gear wrapper for backward compatibility" OFF)
-option(DD4HEP_USE_LCIO "Build lcio extensions" OFF)
+option(DD4HEP_USE_PYROOT "Enable 'Detector Builders' based on PyROOT" ON) # does not work (compile error)
+option(DD4HEP_USE_GEANT4 "Enable the simulation part based on Geant4" ON)
+option(DD4HEP_USE_GEAR "Build gear wrapper for backward compatibility" ON)
+option(DD4HEP_USE_LCIO "Build lcio extensions" ON)
option(BUILD_TESTING "Enable and build tests" ON)
option(DD4HEP_USE_CXX14 "Build DD4hep using c++14" OFF)
option(CMAKE_MACOSX_RPATH "Build with rpath on macos" ON)
diff --git a/DDG4/include/DDG4/Geant4InputAction.h b/DDG4/include/DDG4/Geant4InputAction.h
index 3a63bdfd9046d1899b6bfb67a67110b2642ea9a8..a1b37eac670f65fddc26aa4b5f29f4b3fb97679a 100644
--- a/DDG4/include/DDG4/Geant4InputAction.h
+++ b/DDG4/include/DDG4/Geant4InputAction.h
@@ -47,6 +47,7 @@ namespace DD4hep {
typedef Geant4Vertex Vertex;
typedef Geant4Particle Particle;
typedef std::vector<Particle*> Particles;
+ typedef std::vector<Vertex*> Vertices;
/// Status codes of the event reader object. Anything with NOT low-bit set is an error.
enum EventReaderStatus {
EVENT_READER_ERROR=0,
@@ -89,7 +90,7 @@ namespace DD4hep {
/** The additional argument
*/
virtual EventReaderStatus readParticles(int event_number,
- Vertex& primary_vertex,
+ Vertices& vertices,
Particles& particles) = 0;
};
@@ -109,6 +110,7 @@ namespace DD4hep {
typedef Geant4Vertex Vertex;
typedef Geant4Particle Particle;
typedef std::vector<Particle*> Particles;
+ typedef std::vector<Vertex*> Vertices;
protected:
/// Property: input file
std::string m_input;
@@ -127,7 +129,9 @@ namespace DD4hep {
public:
/// Read an event and return a LCCollectionVec of MCParticles.
- int readParticles(int event_number, Vertex& primary_vertex, Particles& particles);
+ int readParticles(int event_number,
+ Vertices& vertices,
+ Particles& particles);
/// helper to report Geant4 exceptions
std::string issue(int i) const;
diff --git a/DDG4/lcio/LCIOEventReader.cpp b/DDG4/lcio/LCIOEventReader.cpp
index 1f02e5e9e16cfa1ab7e8f9d75c9ffb02398947e3..edf45ff62b1fdc5dc231532a00e695adce7b6de2 100644
--- a/DDG4/lcio/LCIOEventReader.cpp
+++ b/DDG4/lcio/LCIOEventReader.cpp
@@ -61,7 +61,7 @@ LCIOEventReader::~LCIOEventReader() {
/// Read an event and fill a vector of MCParticles.
LCIOEventReader::EventReaderStatus
LCIOEventReader::readParticles(int event_number,
- Vertex& /* primary_vertex */,
+ Vertices& vertices,
vector<Particle*>& particles)
{
EVENT::LCCollection* primaries = 0;
@@ -75,6 +75,17 @@ LCIOEventReader::readParticles(int event_number,
// check if there is at least one particle
if ( NHEP == 0 ) return EVENT_READER_NO_PRIMARIES;
+ //fg: for now we create exactly one event vertex here ( as before )
+ // and fill it below from the first final state particle
+ Geant4Vertex* vtx = new Geant4Vertex ;
+ vertices.push_back( vtx );
+ vtx->x = 0;
+ vtx->y = 0;
+ vtx->z = 0;
+ vtx->time = 0;
+ bool haveVertex = true ;
+
+
mcpcoll.resize(NHEP,0);
for(int i=0; i<NHEP; ++i ) {
EVENT::MCParticle* p = dynamic_cast<EVENT::MCParticle*>(primaries->getElementAt(i));
@@ -128,6 +139,16 @@ LCIOEventReader::readParticles(int event_number,
cout << " #### WARNING - LCIOInputAction : unknown generator status : "
<< genStatus << " -> ignored ! " << endl;
}
+
+ // fill vertex information from first stable particle
+ if( !haveVertex && genStatus == 1 ){
+ vtx->x = p->vsx ;
+ vtx->y = p->vsy ;
+ vtx->z = p->vsz ;
+ vtx->time = p->time ;
+ haveVertex = false ;
+ }
+
if ( mcp->isCreatedInSimulation() ) status.set(G4PARTICLE_SIM_CREATED);
if ( mcp->isBackscatter() ) status.set(G4PARTICLE_SIM_BACKSCATTER);
if ( mcp->vertexIsNotEndpointOfParent() ) status.set(G4PARTICLE_SIM_PARENT_RADIATED);
diff --git a/DDG4/lcio/LCIOEventReader.h b/DDG4/lcio/LCIOEventReader.h
index 55b3bd56e4c8466ad0aa60fafb5ed96e9830dc66..77790f3469c01bc48392e8efd6a7e6b746ff6dce 100644
--- a/DDG4/lcio/LCIOEventReader.h
+++ b/DDG4/lcio/LCIOEventReader.h
@@ -42,7 +42,7 @@ namespace DD4hep {
/// Read an event and fill a vector of MCParticles.
virtual EventReaderStatus readParticles(int event_number,
- Vertex& primary_vertex,
+ Vertices& vertices,
std::vector<Particle*>& particles);
/// Read an event and return a LCCollectionVec of MCParticles.
virtual EventReaderStatus readParticleCollection(int event_number, EVENT::LCCollection** particles) = 0;
diff --git a/DDG4/plugins/Geant4EventReaderHepEvt.cpp b/DDG4/plugins/Geant4EventReaderHepEvt.cpp
index 73b2a228213b365332fff7ddcd37178f44b483ef..4e19d4af78c2568be9d5da9b169f2778be5f8dbb 100644
--- a/DDG4/plugins/Geant4EventReaderHepEvt.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepEvt.cpp
@@ -47,7 +47,7 @@ namespace DD4hep {
virtual ~Geant4EventReaderHepEvt();
/// Read an event and fill a vector of MCParticles.
virtual EventReaderStatus readParticles(int event_number,
- Vertex& primary_vertex,
+ Vertices& vertices,
std::vector<Particle*>& particles);
virtual EventReaderStatus moveToEvent(int event_number);
virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; }
@@ -129,9 +129,10 @@ Geant4EventReaderHepEvt::moveToEvent(int event_number) {
printout(INFO,"EventReaderHepEvt::moveToEvent","Skipping the first %d events ", event_number );
printout(INFO,"EventReaderHepEvt::moveToEvent","Event number before skipping: %d", m_currEvent );
while ( m_currEvent < event_number ) {
- Geant4Vertex vertex;
std::vector<Particle*> particles;
- EventReaderStatus sc = readParticles(m_currEvent,vertex,particles);
+ Vertices vertices ;
+ EventReaderStatus sc = readParticles(m_currEvent,vertices,particles);
+ for_each(vertices.begin(),vertices.end(),deleteObject<Vertex>);
for_each(particles.begin(),particles.end(),deleteObject<Particle>);
if ( sc != EVENT_READER_OK ) return sc;
//Current event is increased in readParticles already!
@@ -145,9 +146,10 @@ Geant4EventReaderHepEvt::moveToEvent(int event_number) {
/// Read an event and fill a vector of MCParticles.
Geant4EventReader::EventReaderStatus
Geant4EventReaderHepEvt::readParticles(int /* event_number */,
- Vertex& /* primary_vertex */,
+ Vertices& vertices,
vector<Particle*>& particles) {
+
// First check the input file status
if ( !m_input.good() || m_input.eof() ) {
return EVENT_READER_IO_ERROR;
@@ -165,23 +167,33 @@ Geant4EventReaderHepEvt::readParticles(int /* event_number */,
// -> FG: EOF is not an exception as it happens for every file at the end !
return EVENT_READER_IO_ERROR;
}
+
+ //fg: for now we create exactly one event vertex here ( as before )
+ // and fill it below from the first final state particle
+ Geant4Vertex* vtx = new Geant4Vertex ;
+ vertices.push_back( vtx );
+ vtx->x = 0;
+ vtx->y = 0;
+ vtx->z = 0;
+ vtx->time = 0;
+ bool haveVertex = false ;
//
// 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
+ int ISTHEP(0); // status code
+ int IDHEP(0); // PDG code
+ int JMOHEP1(0); // first mother
+ int JMOHEP2(0); // last mother
+ int JDAHEP1(0); // first daughter
+ int JDAHEP2(0); // last daughter
+ double PHEP1(0); // px in GeV/c
+ double PHEP2(0); // py in GeV/c
+ double PHEP3(0); // pz in GeV/c
+ double PHEP4(0); // energy in GeV
+ double PHEP5(0); // mass in GeV/c**2
+ double VHEP1(0); // x vertex position in mm
+ double VHEP2(0); // y vertex position in mm
+ double VHEP3(0); // z vertex position in mm
+ double VHEP4(0); // production time in mm/c
vector<int> daughter1;
vector<int> daughter2;
@@ -201,8 +213,9 @@ Geant4EventReaderHepEvt::readParticles(int /* event_number */,
if(m_input.eof())
return EVENT_READER_IO_ERROR;
+
//
- // Create a MCParticle and fill it from stdhep info
+ // create a MCParticle and fill it from stdhep info
Particle* p = new Particle(IHEP);
PropertyMask status(p->status);
//
@@ -218,14 +231,18 @@ Geant4EventReaderHepEvt::readParticles(int /* event_number */,
p->mass = PHEP5*GeV;
//
// Vertex
- // (missing information in HEPEvt files)
- p->vsx = 0.0;
- p->vsy = 0.0;
- p->vsz = 0.0;
+ p->vsx = VHEP1*mm;
+ p->vsy = VHEP2*mm;
+ p->vsz = VHEP3*mm;
+ // endpoint (missing information in HEPEvt files)
p->vex = 0.0;
p->vey = 0.0;
p->vez = 0.0;
//
+ // Creation time (note the units [1/c_light])
+ p->time = VHEP4*ns;
+ p->properTime = VHEP4*ns;
+ //
// Generator status
// Simulator status 0 until simulator acts on it
p->status = 0;
@@ -234,11 +251,16 @@ Geant4EventReaderHepEvt::readParticles(int /* event_number */,
else if ( ISTHEP == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
else if ( ISTHEP == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
else status.set(G4PARTICLE_GEN_DOCUMENTATION);
- //
- // Creation time (note the units [1/c_light])
- // (No information in HEPEvt files)
- p->time = 0.0;
- p->properTime = 0.0;
+
+ // fill vertex information from first stable particle
+ if( !haveVertex && ISTHEP == 1 ){
+ vtx->x = p->vsx ;
+ vtx->y = p->vsy ;
+ vtx->z = p->vsz ;
+ vtx->time = p->time ;
+ haveVertex = true ;
+ }
+
//
// Keep daughters information for later
daughter1.push_back(JDAHEP1);
@@ -309,6 +331,20 @@ Geant4EventReaderHepEvt::readParticles(int /* event_number */,
if ( !part.findParent(mcp) ) part.addParent(mcp);
}
} // End second loop over particles
+
+
+ //--- need a third loop to add particles to the vertex
+ for(size_t i=0; i<particles.size(); ++i ) {
+ Geant4ParticleHandle p(particles[i]);
+ if ( p->parents.size() == 0 ) {
+ PropertyMask status(p->status);
+ if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
+ vtx->in.insert(p->id); // Beam particles and primary quarks etc.
+ else
+ vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
+ }
+ }
+
++m_currEvent;
return EVENT_READER_OK;
}
diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp
index c06566eeb47a115c2a17e1b5d0730fb5bf2c6a94..454fba2802ab973bb4a30cea9e41c646e98e0a78 100644
--- a/DDG4/plugins/Geant4EventReaderHepMC.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp
@@ -48,7 +48,7 @@ namespace DD4hep {
virtual ~Geant4EventReaderHepMC();
/// Read an event and fill a vector of MCParticles.
virtual EventReaderStatus readParticles(int event_number,
- Vertex& primary_vertex,
+ Vertices& vertices,
std::vector<Particle*>& particles) override;
virtual EventReaderStatus moveToEvent(int event_number) override;
virtual EventReaderStatus skipEvent() override { return EVENT_READER_OK; }
@@ -199,14 +199,25 @@ Geant4EventReaderHepMC::moveToEvent(int event_number) {
/// Read an event and fill a vector of MCParticles.
Geant4EventReaderHepMC::EventReaderStatus
Geant4EventReaderHepMC::readParticles(int /* ev_id */,
- Vertex& primary_vertex,
+ Vertices& vertices,
Particles& output) {
+
+ //fg: for now we create exactly one event vertex here ( as before )
+ // this needs revisiting as HepMC allows to have more than one vertex ...
+ Geant4Vertex* primary_vertex = new Geant4Vertex ;
+ vertices.push_back( primary_vertex );
+ primary_vertex->x = 0;
+ primary_vertex->y = 0;
+ primary_vertex->z = 0;
+
if ( !m_events->ok() ) {
return EVENT_READER_IO_ERROR;
}
else if ( m_events->read() ) {
EventStream::Particles& parts = m_events->particles();
- Position pos(primary_vertex.x,primary_vertex.y,primary_vertex.z);
+
+ Position pos(primary_vertex->x,primary_vertex->y,primary_vertex->z);
+
output.reserve(parts.size());
transform(parts.begin(),parts.end(),back_inserter(output),reference2nd(parts));
m_events->clear();
@@ -231,6 +242,16 @@ Geant4EventReaderHepMC::readParticles(int /* ev_id */,
p->daughters.size(),
p->parents.size());
//output.push_back(p);
+
+ //ad particles to the 'primary vertex'
+ if ( p->parents.size() == 0 ) {
+ PropertyMask status(p->status);
+ if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
+ primary_vertex->in.insert(p->id); // Beam particles and primary quarks etc.
+ else
+ primary_vertex->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
+ }
+
}
++m_currEvent;
return EVENT_READER_OK;
diff --git a/DDG4/src/Geant4InputAction.cpp b/DDG4/src/Geant4InputAction.cpp
index f396abb1225596f8e6544f35a3a6e72c5a10f484..216b93e518b8d9514d287142685c06cc1be5c95e 100644
--- a/DDG4/src/Geant4InputAction.cpp
+++ b/DDG4/src/Geant4InputAction.cpp
@@ -24,6 +24,8 @@
using namespace std;
using namespace DD4hep::Simulation;
typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+typedef Geant4InputAction::Vertices Vertices ;
+
/// Initializing constructor
Geant4EventReader::Geant4EventReader(const std::string& nam)
@@ -42,9 +44,10 @@ Geant4EventReader::EventReaderStatus Geant4EventReader::skipEvent() {
return EVENT_READER_OK;
}
std::vector<Particle*> particles;
- Geant4Vertex vertex;
+ Vertices vertices ;
+
++m_currEvent;
- EventReaderStatus sc = readParticles(m_currEvent,vertex,particles);
+ EventReaderStatus sc = readParticles(m_currEvent,vertices,particles);
for_each(particles.begin(),particles.end(),deleteObject<Particle>);
return sc;
}
@@ -100,7 +103,7 @@ string Geant4InputAction::issue(int i) const {
/// Read an event and return a LCCollection of MCParticles.
int Geant4InputAction::readParticles(int evt_number,
- Vertex& prim_vertex,
+ Vertices& vertices,
std::vector<Particle*>& particles)
{
int evid = evt_number + m_firstEvent;
@@ -140,7 +143,9 @@ int Geant4InputAction::readParticles(int evt_number,
except("Error when reading file %s.", m_input.c_str());
return status;
}
- status = m_reader->readParticles(evid, prim_vertex, particles);
+ status = m_reader->readParticles(evid, vertices, particles);
+
+
if ( Geant4EventReader::EVENT_READER_OK != status ) {
string msg = issue(evid)+"Error when moving to event - may be end of file.";
if ( m_abort ) {
@@ -158,14 +163,10 @@ void Geant4InputAction::operator()(G4Event* event) {
vector<Particle*> primaries;
Geant4Event& evt = context()->event();
Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
- dd4hep_ptr<Geant4Vertex> vertex(new Geant4Vertex());
+ Vertices vertices ;
int result;
- vertex->x = 0;
- vertex->y = 0;
- vertex->z = 0;
- vertex->time = 0;
- result = readParticles(m_currentEventNumber, *(vertex.get()), primaries);
+ result = readParticles(m_currentEventNumber, vertices, primaries);
event->SetEventID(m_firstEvent + m_currentEventNumber);
++m_currentEventNumber;
@@ -176,13 +177,21 @@ void Geant4InputAction::operator()(G4Event* event) {
Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
prim->add(m_mask, inter);
+
// check if there is at least one particle
if ( primaries.empty() ) return;
- print("+++ Particle interaction with %d generator particles ++++++++++++++++++++++++",
- int(primaries.size()));
- Geant4Vertex* vtx = vertex.get();
- inter->vertices.insert(make_pair(m_mask,vertex.release())); // Move vertex ownership
+ // check if there is at least one primary vertex
+ if ( vertices.empty() ) return;
+
+ print("+++ Particle interaction with %d generator particles and %d vertices ++++++++++++++++++++++++",
+ int(primaries.size()), int(vertices.size()) );
+
+
+ for(size_t i=0; i<vertices.size(); ++i ) {
+ inter->vertices.insert(make_pair(m_mask,vertices[i]));
+ }
+
// build collection of MCParticles
for(size_t i=0; i<primaries.size(); ++i ) {
Geant4ParticleHandle p(primaries[i]);
@@ -192,13 +201,19 @@ void Geant4InputAction::operator()(G4Event* event) {
p->psy = mom_scale*p->psy;
p->psz = mom_scale*p->psz;
- if ( p->parents.size() == 0 ) {
- if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
- vtx->in.insert(p->id); // Beam particles and primary quarks etc.
- else
- vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
- }
+ //FIXME: this needs to be done now in the readers ...
+ // // if ( p->parents.size() == 0 ) {
+ // // if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
+ // // vtx->in.insert(p->id); // Beam particles and primary quarks etc.
+ // // else
+ // // vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
+ // // }
+
+
inter->particles.insert(make_pair(p->id,p));
p.dumpWithMomentumAndVertex(outputLevel()-1,name(),"->");
}
+
+
+
}
diff --git a/DDTest/src/test_EventReaders.cc b/DDTest/src/test_EventReaders.cc
index c8a43f8fd0d026360762a62dd88d8e45f860b754..2dc456299f8eacd1c60d84b31f73a960e7c6d07c 100644
--- a/DDTest/src/test_EventReaders.cc
+++ b/DDTest/src/test_EventReaders.cc
@@ -60,12 +60,9 @@ int main(int argc, char** argv ){
thisReader->moveToEvent(1);
test( thisReader->currentEventNumber() == 1 , readerType + std::string("Event Number after Skip") );
std::vector<Particle*> particles;
- Vertex vertex;
- vertex.x = 0;
- vertex.y = 0;
- vertex.z = 0;
- vertex.time = 0;
- DD4hep::Simulation::Geant4EventReader::EventReaderStatus sc = thisReader->readParticles(3,vertex,particles);
+ std::vector<Vertex*> vertices ;
+
+ DD4hep::Simulation::Geant4EventReader::EventReaderStatus sc = thisReader->readParticles(3,vertices,particles);
std::for_each(particles.begin(),particles.end(),DD4hep::deleteObject<Particle>);
test( thisReader->currentEventNumber() == 2 && sc == DD4hep::Simulation::Geant4EventReader::EVENT_READER_OK,
readerType + std::string("Event Number Read") );
diff --git a/doc/release.notes b/doc/release.notes
index 2dba0b336f7787d98806807d443bcde8c1b3e131..bc7bd73844a9dc40feebc8b70810021b69f5b458 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -3,6 +3,18 @@
DD4hep ---- Release Notes
=================================
+
+
+Frank Gaede 2017-02-10
+ - allow event readers to create more than one vertex per event
+ this should be possible as most generator formats allow to specify
+ more than one event vertex
+ - changed signature of Geant4EventReader::readParticles(int,Vertex*, Particles&)
+ to Geant4EventReader::readParticles(int,Vertices&, Particles& )
+ - implement in LCIOEventReader, Geant4EventReaderHepEvt and Geant4EventReaderHepMC
+ - for now still one vertex only is created using the first final state particle
+ for HepEvt and LCIO
+
--------
| v00-20 |
--------