diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index 84370551108eb4573207681008b44298368ec14a..3e8ef6111d24a1ad4afe6873d96ead5542e66ca8 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -108,8 +108,11 @@ def run():
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_1.stdhep"
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/qq_gen_128_999.stdhep"
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/smuonLR_PointK_3TeV_BS_noBkg_run0001.stdhep"
- gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbbb_3TeV.stdhep"
+ #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbbb_3TeV.stdhep"
#gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_pi-_5GeV.slcio"
+ gen.Input = "LCIOFileReader|/home/frankm/SW/data/bbbb_3TeV.slcio"
+ #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/FCC-eh.stdhep"
+ #gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt"
gen.OutputLevel = 4 # generator_output_level
gen.MomentumScale = 0.1
gen.Mask = 1
@@ -123,7 +126,7 @@ def run():
gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns)
gen.enableUI()
kernel.generatorAction().adopt(gen)
-
+ """
# Second particle file reader
gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO2");
gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_2.stdhep"
@@ -141,7 +144,7 @@ def run():
gen.Sigma = (2*mm, 1*mm, 1*mm, 0*ns)
gen.enableUI()
kernel.generatorAction().adopt(gen)
-
+ """
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
# Merge all existing interaction records
diff --git a/DDG4/lcio/LCIOInputAction.h b/DDG4/lcio/LCIOInputAction.h
deleted file mode 100644
index 035dee40a783c1983495e56f7db08ad71c31d04b..0000000000000000000000000000000000000000
--- a/DDG4/lcio/LCIOInputAction.h
+++ /dev/null
@@ -1,85 +0,0 @@
-// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
-//====================================================================
-// AIDA Detector description implementation for LCD
-//--------------------------------------------------------------------
-//
-//====================================================================
-#ifndef DD4HEP_DDG4_LCIOINPUTACTION_H
-#define DD4HEP_DDG4_LCIOINPUTACTION_H
-
-// Framework include files
-#include "DD4hep/Primitives.h"
-#include "DDG4/Geant4Particle.h"
-#include "DDG4/Geant4GeneratorAction.h"
-
-// C/C++ include files
-#include <set>
-#include <map>
-#include <vector>
-
-// Forward declarations
-class G4Event;
-namespace IO { class LCReader; }
-namespace UTIL { class LCStdHepRdr; }
-namespace EVENT { class MCParticle; }
-namespace EVENT { class LCCollection; }
-namespace IMPL { class MCParticleImpl; }
-namespace IMPL { class LCCollectionVec; }
-
-/*
- * DD4hep namespace declaration
- */
-namespace DD4hep {
-
- /*
- * Simulation namespace declaration
- */
- namespace Simulation {
-
- // Forward declarations
- class Geant4Particle;
- class LCIOEventReader;
-
- /// Input action reading LCIO files
- /**
- * Concrete implementation of the Geant4 generator action base class
- * populating Geant4 primaries from LCIO and HepStd files.
- *
- * @author P.Kostka (main author)
- * @author M.Frank (code reshuffeling into new DDG4 scheme)
- * @version 1.0
- */
- class LCIOInputAction : public Simulation::Geant4GeneratorAction {
-
- public:
- typedef Geant4Particle Particle;
-
- protected:
- /// Property: input file
- std::string m_input;
- /// Property: SYNCEVT
- int m_SYNCEVT;
- /// Property; interaction mask
- int m_mask;
- /// Property: Momentum downscaler for debugging
- double m_momScale;
- /// Event reader object
- LCIOEventReader* m_reader;
-
- /// Read an event and return a LCCollectionVec of MCParticles.
- EVENT::LCCollection* readParticles(int which);
- /// helper to report Geant4 exceptions
- std::string issue(int i) const;
-
- public:
- /// Standard constructor
- LCIOInputAction(Simulation::Geant4Context* context, const std::string& name);
- /// Default destructor
- virtual ~LCIOInputAction();
- /// Callback to generate primary particles
- virtual void operator()(G4Event*);
- };
- } /* End namespace lcio */
-} /* End namespace DD4hep */
-
-#endif /* DD4HEP_DDG4_LCIOINPUTACTION_H */
diff --git a/DDG4/plugins/Geant4EscapeCounter.cpp b/DDG4/plugins/Geant4EscapeCounter.cpp
index 8bc974034b4befef5c829cf00b1463732bec0134..0cb563facd3324db413e30894594dcd359ea1c85 100644
--- a/DDG4/plugins/Geant4EscapeCounter.cpp
+++ b/DDG4/plugins/Geant4EscapeCounter.cpp
@@ -67,9 +67,7 @@ namespace DD4hep {
#include "G4VProcess.hh"
using namespace std;
-using namespace CLHEP;
using namespace DD4hep;
-using namespace DD4hep::Geometry;
using namespace DD4hep::Simulation;
/// Standard constructor
@@ -90,7 +88,7 @@ Geant4EscapeCounter::~Geant4EscapeCounter() {
/// G4VSensitiveDetector interface: Method for generating hit(s) using the information of G4Step object.
bool Geant4EscapeCounter::process(G4Step* step, G4TouchableHistory* /* history */) {
- Geant4StepHandler h(step);
+ Geant4StepHandler h(step);
Geant4TrackHandler th(h.track);
Geant4TouchableHandler handler(step);
string path = handler.path();
@@ -106,7 +104,7 @@ bool Geant4EscapeCounter::process(G4Step* step, G4TouchableHistory* /* history *
mark(h.track);
print("+++ Track:%4d %8.2f MeV [%s] %s Geant4 path:%s",
- h.trkID(),h.trkEnergy()/MeV,th.name().c_str(),
+ h.trkID(),h.trkEnergy()/CLHEP::MeV,th.name().c_str(),
th.creatorName().c_str(),path.c_str());
// Kill track, so that it does no longer participate in the propagation
h.track->SetTrackStatus(fStopAndKill);
diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp
index 07b35190ce204afc19024b0971e20b2ebdc5bc83..62d95d535d2ee709dc9c3e51a9b8b474a6ae5104 100644
--- a/DDG4/plugins/Geant4EventReaderHepMC.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp
@@ -101,39 +101,39 @@ namespace DD4hep {
class EventStream {
public:
+ typedef std::map<int,Geant4Vertex*> Vertices;
+ typedef std::map<int,Geant4Particle*> Particles;
+
istream& instream;
+
// io information
string key;
double mom_unit, pos_unit;
int io_type;
- float xsection, xsection_err;
- //WeightContainer weights;
- typedef std::map<int,Geant4Vertex*> Vertices;
+
+ float xsection, xsection_err;
+ EventHeader header;
Vertices m_vertices;
- typedef std::map<int,Geant4Particle*> Particles;
Particles m_particles;
- EventHeader header;
EventStream(istream& in) : instream(in), io_type(0)
{ use_default_units(); }
+ /// Check if data stream is in proper state and has data
+ bool ok() const;
+ Geant4Vertex* vertex(int i);
Particles& particles() { return m_particles; }
Vertices& vertices() { return m_vertices; }
- void add_vertex(int id, Geant4Vertex* v)
- { m_vertices.insert(make_pair(id,v)); }
- void add_particle(Geant4Particle* p);
void set_io(int typ, const string& k)
{ io_type = typ; key = k; }
void use_default_units()
- { mom_unit = MeV; pos_unit = mm; }
- Geant4Vertex* vertex(int i);
- void fix_particles();
+ { mom_unit = MeV; pos_unit = mm; }
bool read();
void clear();
};
char get_input(istream& is, istringstream& iline);
- int find_event_end(istream & is);
+ int read_until_event_end(istream & is);
int read_weight_names(EventStream &, istringstream& iline);
int read_particle(EventStream &info, istringstream& iline, Geant4Particle * p);
int read_vertex(EventStream &info, istream& is, istringstream & iline);
@@ -142,6 +142,8 @@ namespace DD4hep {
int read_units(EventStream &info, istringstream & input);
int read_heavy_ion(EventStream &, istringstream & input);
int read_pdf(EventStream &, istringstream & input);
+ Geant4Vertex* vertex(EventStream& info, int i);
+ void fix_particles(EventStream &info);
}
}
}
@@ -166,22 +168,14 @@ Geant4EventReaderHepMC::~Geant4EventReaderHepMC() {
m_events = 0;
m_input.close();
}
-namespace DD4hep{
- template <typename M> ReferenceObjects<typename M::value_type> __reference2nd(M&) {
- return ReferenceObjects<typename M::value_type>();
- }
-}
+
/// Read an event and fill a vector of MCParticles.
Geant4EventReaderHepMC::EventReaderStatus
Geant4EventReaderHepMC::readParticles(int /* ev_id */, Particles& output) {
- // make sure the stream is good
- if ( m_input.eof() || m_input.fail() ) {
- cerr << "streaming input: end of stream found setting badbit." << endl;
- m_input.clear(ios::badbit);
+ if ( !m_events->ok() ) {
return EVENT_READER_IO_ERROR;
}
-
- if ( m_events->read() ) {
+ else if ( m_events->read() ) {
EventStream::Particles& p = m_events->particles();
output.reserve(p.size());
transform(p.begin(),p.end(),back_inserter(output),reference2nd(p));
@@ -202,15 +196,16 @@ Geant4EventReaderHepMC::readParticles(int /* ev_id */, Particles& output) {
return EVENT_READER_IO_ERROR;
}
-
-void HepMC::EventStream::fix_particles() {
+void HepMC::fix_particles(EventStream& info) {
+ EventStream::Particles& parts = info.particles();
+ EventStream::Vertices& verts = info.vertices();
EventStream::Particles::iterator i;
std::set<int>::const_iterator id, ip;
- for(i=m_particles.begin(); i != m_particles.end(); ++i) {
+ for(i=parts.begin(); i != parts.end(); ++i) {
Geant4ParticleHandle p((*i).second);
int end_vtx_id = p->secondaries;
p->secondaries = 0;
- Geant4Vertex* v = vertex(end_vtx_id);
+ Geant4Vertex* v = vertex(info,end_vtx_id);
if ( v ) {
p->vex = v->x;
p->vey = v->y;
@@ -221,17 +216,22 @@ void HepMC::EventStream::fix_particles() {
}
}
EventStream::Vertices::iterator j;
- for(j=m_vertices.begin(); j != m_vertices.end(); ++j) {
+ for(j=verts.begin(); j != verts.end(); ++j) {
Geant4Vertex* v = (*j).second;
for (id=v->in.begin(); id!=v->in.end();++id) {
for (ip=v->out.begin(); ip!=v->out.end();++ip) {
- EventStream::Particles::iterator ipp = m_particles.find(*ip);
+ EventStream::Particles::iterator ipp = parts.find(*ip);
(*ipp).second->parents.insert(*id);
}
}
}
}
+Geant4Vertex* HepMC::vertex(EventStream& info, int i) {
+ EventStream::Vertices::iterator it=info.vertices().find(i);
+ return (it==info.vertices().end()) ? 0 : (*it).second;
+}
+
char HepMC::get_input(istream& is, istringstream& iline) {
char value = is.peek();
if ( !is ) { // make sure the stream is valid
@@ -251,7 +251,7 @@ char HepMC::get_input(istream& is, istringstream& iline) {
return iline ? value : -1;
}
-int HepMC::find_event_end(istream & is) {
+int HepMC::read_until_event_end(istream & is) {
string line;
while ( is ) {
char val = is.peek();
@@ -281,7 +281,7 @@ int HepMC::read_weight_names(EventStream&, istringstream&) {
cout << "debug: attempting to read past the end of the named weight line " << endl;
cout << "debug: We should never get here" << endl;
cout << "debug: Looking for the end of this event" << endl;
- find_event_end(is);
+ read_until_event_end(is);
}
i2 = line.find("\"",i1+1);
name = line.substr(i1+1,i2-i1-1);
@@ -297,63 +297,73 @@ int HepMC::read_weight_names(EventStream&, istringstream&) {
return 1;
}
-int HepMC::read_particle(EventStream &info, istringstream& iline, Geant4Particle * p) {
+int HepMC::read_particle(EventStream &info, istringstream& input, Geant4Particle * p) {
float ene = 0., theta = 0., phi = 0;
int size = 0;
// check that the input stream is still OK after reading item
- iline >> p->id >> p->pdgID >> p->psx >> p->psy >> p->psz >> ene;
- if ( !iline )
+ input >> p->id >> p->pdgID >> p->psx >> p->psy >> p->psz >> ene;
+ p->psx /= info.mom_unit;
+ p->psy /= info.mom_unit;
+ p->psz /= info.mom_unit;
+ ene /= info.mom_unit;
+ if ( !input )
return 0;
- else if( info.io_type != ascii )
- iline >> p->mass;
- else
+ else if ( info.io_type != ascii ) {
+ input >> p->mass;
+ p->mass /= info.mom_unit;
+ }
+ else {
p->mass = sqrt(fabs(ene*ene - p->psx*p->psx + p->psy*p->psy + p->psz*p->psz));
-
+ }
// Reuse here the secondaries to store the end-vertex ID
- iline >> p->status >> theta >> phi >> p->secondaries >> size;
- if(!iline)
+ input >> p->status >> theta >> phi >> p->secondaries >> size;
+ if(!input)
return 0;
// read flow patterns if any exist
for (int i = 0; i < size; ++i ) {
- iline >> p->colorFlow[0] >> p->colorFlow[1];
- if(!iline) return 0;
+ input >> p->colorFlow[0] >> p->colorFlow[1];
+ if(!input) return 0;
}
return 1;
}
-int HepMC::read_vertex(EventStream &info, istream& is, istringstream & iline) {
+int HepMC::read_vertex(EventStream &info, istream& is, istringstream & input) {
int id=0, dummy=0, num_orphans_in=0, num_particles_out=0, weights_size=0;
vector<float> weights;
Geant4Vertex* v = new Geant4Vertex();
Geant4Particle* p;
- iline >> id >> dummy >> v->x >> v->y >> v->z >> v->time
+ input >> id >> dummy >> v->x >> v->y >> v->z >> v->time
>> num_orphans_in >> num_particles_out >> weights_size;
- if(!iline)
+ if(!input)
return 0;
+ v->x /= info.pos_unit;
+ v->y /= info.pos_unit;
+ v->z /= info.pos_unit;
weights.resize(weights_size);
for (int i1 = 0; i1 < weights_size; ++i1) {
- iline >> weights[i1];
- if(!iline)
+ input >> weights[i1];
+ if(!input)
return 0;
}
- info.add_vertex(id,v);
+ info.vertices().insert(make_pair(id,v));
//cout << "Add Vertex:" << id << endl;
for(char value = is.peek(); value=='P'; value=is.peek()) {
- value = get_input(is,iline);
- if( !iline || value < 0 )
+ value = get_input(is,input);
+ if( !input || value < 0 )
return 0;
- read_particle(info, iline,p = new Geant4Particle());
- if(!iline) {
+ read_particle(info, input,p = new Geant4Particle());
+ if(!input) {
cerr << "Failed to read particle!" << endl;
delete p;
return 0;
}
- info.add_particle(p);
+ p->id = info.particles().size();
+ info.particles().insert(make_pair(p->id,p));
p->pex = p->psx;
p->pey = p->psy;
p->pez = p->psz;
@@ -502,14 +512,14 @@ int HepMC::read_pdf(EventStream &, istringstream & input) {
return input.fail() ? 0 : 1;
}
-void HepMC::EventStream::add_particle(Geant4Particle* p) {
- p->id = m_particles.size();
- m_particles.insert(make_pair(p->id,p));
-}
-
-Geant4Vertex* HepMC::EventStream::vertex(int i) {
- Vertices::iterator it=m_vertices.find(i);
- return (it==m_vertices.end()) ? 0 : (*it).second;
+/// Check if data stream is in proper state and has data
+bool HepMC::EventStream::ok() const {
+ // make sure the stream is good
+ if ( instream.eof() || instream.fail() ) {
+ instream.clear(ios::badbit);
+ return false;
+ }
+ return true;
}
void HepMC::EventStream::clear() {
@@ -519,20 +529,21 @@ void HepMC::EventStream::clear() {
bool HepMC::EventStream::read() {
EventStream& info = *this;
+ bool event_read = false;
- // OK - now ready to start reading the event, so set the header flag
- // The flag will be set to false when we reach the end of the header
- m_particles.clear();
- m_vertices.clear();
+ releaseObjects(vertices())();
+ releaseObjects(particles())();
- bool event_read = false;
while( instream.good() ) {
char value = instream.peek();
- if ( value == 'E' && event_read ) break;
- else if ( event_read ) break;
- else if ( instream.eof() ) return false;
-
istringstream input_line;
+
+ if ( value == 'E' && event_read ) break;
+ else if ( instream.eof() ) return false;
+ else if ( value=='#' || ::isspace(value) ) {
+ get_input(instream,input_line);
+ continue;
+ }
value = get_input(instream,input_line);
// On failure switch to end
@@ -577,6 +588,9 @@ bool HepMC::EventStream::read() {
instream.clear(ios::badbit);
return false;
}
+ else if ( iotype != 0 ) {
+ break;
+ }
continue;
}
case 'E': // deal with the event line
@@ -622,11 +636,11 @@ bool HepMC::EventStream::read() {
printout(WARNING,"HepMC::EventStream","+++ Skip event with ID: %d",this->header.id);
releaseObjects(vertices())();
releaseObjects(particles())();
- find_event_end(instream);
+ read_until_event_end(instream);
event_read = false;
if ( instream.eof() ) return false;
}
- this->fix_particles();
+ fix_particles(info);
releaseObjects(vertices())();
return true;
}