--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/EventSourceManager.hh	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/EventSourceManager.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -19,6 +19,7 @@

 #include "G4UnitsTable.hh"
 
 namespace slic {

+

 /**
  * @class EventSourceManager
  * @brief This singleton is responsible for managing physics event sources.

@@ -31,7 +32,11 @@

      * Enum for different source types.
      */
     enum ESourceType {

-        eUnknown = 0, eStdHep = 1, eLCIO, eGPS, eParticleGun

+        eUnknown = 0,
+        eStdHep = 1,
+        eLCIO,
+        eGPS,
+        eParticleGun

     };
 
 public:

@@ -286,7 +291,7 @@

 
     public:
 

-        virtual void GeneratePrimaryVertex(G4Event* evt);

+        virtual void generateNextEvent(G4Event* evt);

 
     private:
 

@@ -314,11 +319,11 @@

         bool m_newSource;
 
         // gen strings

-        static std::string m_stdhepStr;
-        static std::string m_lcioStr;
-        static std::string m_gpsStr;
-        static std::string m_gunStr;
-        static std::string m_unknownStr;

+        static std::string STDHEP;
+        static std::string LCIO;
+        static std::string GPS;
+        static std::string PARTICLE_GUN;
+        static std::string UNKNOWN;

 
         // source type
         ESourceType m_sourceType;

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioEventSource.hh	                        (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioEventSource.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -0,0 +1,56 @@

+// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioEventSource.hh,v 1.9 2012-11-27 19:32:17 jeremy Exp $
+
+#ifndef SLIC_LCIOEVENTSOURCE_HH
+#define SLIC_LCIOEVENTSOURCE_HH 1
+
+// slic
+#include "EventSourceWithInputFile.hh"
+#include "MCParticleGenerator.hh"
+#include "LcioManager.hh"
+
+// lcio
+#include "EVENT/LCEvent.h"
+#include "EVENT/LCCollection.h"
+#include "IOIMPL/LCFactory.h"
+#include "IO/LCReader.h"
+
+namespace slic {
+
+/**
+ * @class LcioEventSource
+ * @brief Generates events from the MCParticle collection of an LCIO file
+ */
+class LcioEventSource: public EventSourceWithInputFile {
+
+public:
+
+	LcioEventSource(const std::string& fname = "");
+
+	virtual ~LcioEventSource();
+
+public:
+
+	// open the current file
+	void open();
+
+	// close the current file
+	void close();
+
+	// read the next event
+	void readNextEvent();
+
+	void dumpCurrentEvent();
+
+	void generate(G4Event* anEvent);
+
+	void beginRun(const G4Run* aRun);
+
+	void beginEvent(const G4Event* anEvent);
+
+private:
+	IO::LCReader* m_reader;
+	MCParticleGenerator* m_eventGenerator;
+};
+}
+
+#endif

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/MCParticleGenerator.hh	                        (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/MCParticleGenerator.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -0,0 +1,54 @@

+#ifndef slic_MCParticleGenerator_hh_
+#define slic_MCParticleGenerator_hh_ 1
+
+// LCIO
+#include "EVENT/LCEvent.h"
+#include "EVENT/MCParticle.h"
+#include "IO/LCReader.h"
+#include "IMPL/LCCollectionVec.h"
+
+// Geant4
+#include "G4Event.hh"
+#include "G4PrimaryParticle.hh"
+
+// STL
+#include <set>
+
+using EVENT::MCParticle;
+using EVENT::LCCollection;
+using EVENT::LCEvent;
+using IMPL::LCCollectionVec;
+using IO::LCReader;
+
+namespace slic {
+
+class MCParticleManager;
+
+/**
+ * @class MCParticle
+ * @brief Converts from MCParticle collection into a G4Event
+ */
+class MCParticleGenerator {
+
+public:
+
+    MCParticleGenerator(G4String eventFile);
+
+    virtual ~MCParticleGenerator();
+
+    void generateEvent(G4Event*);
+
+    void readNextEvent();
+
+    LCEvent* getCurrentEvent();
+
+private:
+
+    LCEvent* _lcevent;
+    LCReader* _reader;
+    MCParticleManager* _particleManager;
+};
+
+}
+
+#endif

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/MCParticleManager.hh	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/MCParticleManager.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -2,15 +2,18 @@

 #define slic_MCParticleManager_hh_ 1
 
 // SLIC

+#include "Module.hh"

 #include "Singleton.hh"
 
 // LCIO

+#include "EVENT/LCCollection.h"

 #include "IMPL/LCCollectionVec.h"
 #include "IMPL/MCParticleImpl.h"
 #include "EVENT/MCParticle.h"
 
 // Geant4
 #include "globals.hh"

+#include "G4Event.hh"

 #include "G4PrimaryParticle.hh"
 
 // STL

@@ -20,58 +23,82 @@

 using std::map;
 using std::set;
 

+using EVENT::LCCollection;

 using EVENT::MCParticle;
 using IMPL::MCParticleImpl;
 using IMPL::LCCollectionVec;
 
 namespace slic {
 

-typedef map<G4int, MCParticle*> MCParticleTrackIDMap;
-typedef map<MCParticle*, G4PrimaryParticle*> MCParticleGeant4Map;
-typedef MCParticleGeant4Map::iterator MCParticleGeant4MapIterator;
-typedef set<MCParticle*> VisitedMCParticles;

+typedef map<G4int, MCParticle*> TrackIDMap;
+typedef map<MCParticle*, G4PrimaryParticle*> PrimaryParticleMap;
+typedef PrimaryParticleMap::iterator PrimaryParticleMapIterator;
+typedef set<MCParticle*> MCParticleSet;

-class MCParticleManager : public Singleton<MCParticleManager> {

+/**
+ * @class MCParticleManager
+ *
+ * @note
+ * This class is responsible for management of MCParticle information
+ * and provides central utilities for use by event generation classes
+ * such as the MCParticleGenerator and StdHepGenerator.
+ */
+class MCParticleManager : public Singleton<MCParticleManager>, Module {

 
 public:
 

+    /**
+     * Class constructor.
+     */

     MCParticleManager()

-        : _mcpVec(0) {

+        : Module("MCParticleManager"), _mcpVec(0) {

     }
 

+    /**
+     * Class destructor.
+     */

     virtual ~MCParticleManager() {}
 

+    /**
+     * Get the current MCParticle collection.
+     */

     LCCollectionVec* getMCParticleCollection() {
         return _mcpVec;
     }
 

-    void setMCParticleCollection(LCCollectionVec* mcpVec) {
-        _mcpVec = mcpVec;
-    }
-

     void addMCParticleTrackID(G4int trackId, MCParticle* particle) {
         if (trackId <= 0) {
             return;
         }

-        _mcpMap[trackId] = particle;

+        _trackIDMap[trackId] = particle;

     }
 
     MCParticle* findMCParticle(G4int trackId) {

-        return _mcpMap[trackId];

+        return _trackIDMap[trackId];

     }
 

-    void clearMaps() {
-        _mcpMap.clear();
-        _particleMap.clear();
-        _visitedParticles.clear();

+    PrimaryParticleMap* getPrimaryParticleMap() {
+        return &_primaryMap;

     }
 

-    static MCParticleGeant4Map* getMCParticleToGeant4Map() {
-        return &_particleMap;

+    void reset() {
+        _trackIDMap.clear();
+        _primaryMap.clear();
+        _visitedParticles.clear();

     }
 

-    static std::set<MCParticle*>* getVisitedParticles() {
-        return &_visitedParticles;

+    void generateEvent(LCCollectionVec* mcparticles, G4Event* event);
+
+    void applyLorentzTransformation(LCCollection*, const G4double);
+
+    void applyZSmearing(LCCollection*, const G4double);
+
+private:
+
+    std::set<G4PrimaryParticle*> getRelevantParticles(MCParticle* p);
+
+    void setMCParticleCollection(LCCollectionVec* mcpVec) {
+        _mcpVec = mcpVec;

     }
 
 private:

@@ -79,20 +106,22 @@

     /**
      * Map of MCParticle to G4PrimaryParticle.
      */

-    static MCParticleGeant4Map _particleMap;

+    PrimaryParticleMap _primaryMap;

 
     /**
      * Set keeping track of which MCParticles have been processed in generator.
      */

-    static VisitedMCParticles _visitedParticles;

+    MCParticleSet _visitedParticles;

 
     /**
      * Map of MCParticle to track ID.
      */

-    MCParticleTrackIDMap _mcpMap;

+    TrackIDMap _trackIDMap;

 
     /**

-     * The output MCParticle collection.

+     * The output MCParticle collection.  This pointer is "taken" from the input reader,
+     * because it is added to the output event, which then automatically deletes it once the
+     * event is written (I think?).

      */
     LCCollectionVec* _mcpVec;
 };

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/PrimaryGeneratorAction.hh	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/PrimaryGeneratorAction.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -3,10 +3,11 @@

 #ifndef SLIC_PRIMARYGENERATORACTION_HH
 #define SLIC_PRIMARYGENERATORACTION_HH 1
 

-// slic

+// SLIC
+#include "EventSourceManager.hh"

 #include "Module.hh"
 

-// G4

+// Geant4

 #include "G4VUserPrimaryGeneratorAction.hh"
 #include "G4GeneralParticleSource.hh"
 #include "G4ParticleGun.hh"

@@ -17,7 +18,6 @@

 /**
  * @class PrimaryGeneratorAction
  * @brief Implementation of G4VUserPrimaryGeneratorAction.

- * @note  Delegates all control to EventSourceManager.

  */
 class PrimaryGeneratorAction: public G4VUserPrimaryGeneratorAction, public Module {
 

@@ -41,26 +41,16 @@

      */
     virtual void GeneratePrimaries(G4Event *anEvent);
 

-    /**
-     * Apply Lorentz transformation to event.  The lorentz transformation
-     * is retrieved from the singleton EventSourceManager.
-     * @param[in] anEvent The current Geant4 event.
-     */
-    //void applyLorentzTransformation(G4Event *anEvent);
-
-    /**
-     * Apply Z smearing to event.  The Z spread parameter is retrieved
-     * from the singleton EventSourceManager.
-     * @param[in] anEvent The current Geant4 event.
-     */
-    //void zsmear(G4Event* anEvent);
-

 private:
 
     /**
      * Begin the beginning of event message.
      */
     void printBeginEventMessage(G4Event* anEvent);

+
+private:
+
+    EventSourceManager* _manager;

 };
 }
 

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepEventSource.hh	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepEventSource.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -64,14 +64,7 @@

      */
     void beginEvent(const G4Event* anEvent);
 

-    /**
-     * Get the lStdHep reader associated with this source.
-     * @return The StdHep reader.
-     */
-    lStdHep* getStdHepReader();
-

 private:

-    StdHepLoader* m_loader;

     StdHepGenerator* m_eventGenerator;
 };
 }

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepGenerator.hh	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepGenerator.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -7,6 +7,7 @@

 // LCIO
 #include "UTIL/LCStdHepRdr.h"
 #include "EVENT/MCParticle.h"

+#include "IMPL/LCCollectionVec.h"

 
 // Geant4
 #include "G4Event.hh"

@@ -18,6 +19,7 @@

 using UTIL::LCStdHepRdr;
 using EVENT::MCParticle;
 using EVENT::LCCollection;

+using IMPL::LCCollectionVec;

 
 namespace slic {
 

@@ -36,14 +38,13 @@

 
     void generateEvent(G4Event*);
 

-private:

+    void readNextEvent();

-    std::set<G4PrimaryParticle*> getRelevantParticles(MCParticle* p);

+    LCCollectionVec* getCurrentParticleCollection();

 
 private:
     LCStdHepRdr* _reader;

-    MCParticleManager* _particleManager;
-    VisitedMCParticles* _visitedParticles;

+    LCCollectionVec* _particles;

 };
 
 }

--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackManager.hh	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackManager.hh	2014-01-15 00:17:32 UTC (rev 2876)
@@ -31,6 +31,7 @@

 
     virtual ~TrackManager() {
         delete _trackSummaries;

+        delete _trackSummaryMap;

     }
 
 public:

@@ -39,13 +40,9 @@

         return _trackSummaries;
     }
 

-    void resetTrackSummaries() {
-        _trackSummaries->clear();
-    }
-

     void saveTrackSummaries(const G4Event*, LCEvent*);
 

-    void clearTrackSummaries() {

+    void reset() {

         /* Clear the vector. */
         for (TrackSummaryVector::iterator it = _trackSummaries->begin();
                 it != _trackSummaries->end();

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventAction.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventAction.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -30,12 +30,6 @@

 
 void EventAction::BeginOfEventAction(const G4Event *anEvent) {
 

-    // Reset the global TrackSummaryVector used for MCParticle output.
-    TrackManager::instance()->resetTrackSummaries();
-
-    /* Clear the MCParticle map from last event. */
-    MCParticleManager::instance()->clearMaps();
-

     /* Create an empty LCEvent. */
     LcioManager::instance()->createLCEvent();
 }

@@ -54,9 +48,6 @@

 	/* Begin of event processing for the current event generator. */
 	EventSourceManager::instance()->endEvent(anEvent);
 

-	/* Clear the track information for this event. */
-	TrackManager::instance()->clearTrackSummaries();
-

 	/* Stop the event timer. */
 	stopEventTimer();
 

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventSourceManager.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventSourceManager.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -1,13 +1,16 @@

 // $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/EventSourceManager.cc,v 1.26 2013-06-26 01:57:02 jeremy Exp $

+#include "EventSourceManager.hh"

 
 // SLIC

-#include "EventAction.hh"
-#include "EventSourceManager.hh"

 #include "FileUtil.hh"

-#include "GPSEventSource.hh"
-#include "ParticleGunEventSource.hh"

+
+#include "EventAction.hh"

 #include "RunManager.hh"
 #include "SlicApplication.hh"

+
+#include "GPSEventSource.hh"
+#include "LcioEventSource.hh"
+#include "ParticleGunEventSource.hh"

 #include "StdHepEventSource.hh"
 
 // LCDD

@@ -15,11 +18,11 @@

 
 namespace slic {
 

-std::string EventSourceManager::m_stdhepStr = std::string("stdhep");
-std::string EventSourceManager::m_lcioStr = std::string("lcio");
-std::string EventSourceManager::m_gpsStr = std::string("gps");
-std::string EventSourceManager::m_gunStr = std::string("gun");
-std::string EventSourceManager::m_unknownStr = std::string("unknown");

+std::string EventSourceManager::STDHEP = std::string("stdhep");
+std::string EventSourceManager::LCIO = std::string("lcio");
+std::string EventSourceManager::GPS = std::string("gps");
+std::string EventSourceManager::PARTICLE_GUN = std::string("gun");
+std::string EventSourceManager::UNKNOWN = std::string("unknown");

 
 EventSourceManager::EventSourceManager() :
 		Module("EventSourceManager"), m_currentEventSource(0), m_filename(""), m_fileIsSet(false),

@@ -97,16 +100,16 @@

 
 const std::string& EventSourceManager::getSourceNameFromType(ESourceType egt) const {
 	if (egt == eStdHep) {

-		return m_stdhepStr;

+		return STDHEP;

 	} else if (egt == eLCIO) {

-		return m_lcioStr;

+		return LCIO;

 	} else if (egt == eGPS) {

-		return m_gpsStr;

+		return GPS;

 	} else if (egt == eParticleGun) {

-		return m_gunStr;

+		return PARTICLE_GUN;

 	}
 

-	return m_unknownStr;

+	return UNKNOWN;

 }
 
 EventSourceManager::ESourceType EventSourceManager::getSourceTypeFromName(const std::string& s) const {

@@ -114,13 +117,13 @@

 
 	ESourceType egt = eUnknown;
 

-	if (sl == m_stdhepStr) {

+	if (sl == STDHEP) {

 		egt = eStdHep;

-	} else if (sl == m_gpsStr) {

+	} else if (sl == GPS) {

 		egt = eGPS;

-	} else if (sl == m_gunStr) {

+	} else if (sl == PARTICLE_GUN) {

 		egt = eParticleGun;

-	} else if (sl == m_lcioStr) {

+	} else if (sl == LCIO) {

 		egt = eLCIO;
 	}
 

@@ -174,10 +177,9 @@

 
 EventSource* EventSourceManager::createEventSource(ESourceType st) {
 	EventSource* src = 0;

-	/*if (st == eLCIO) {

+	if (st == eLCIO) {

 		src = new LcioEventSource(getFilename());

-	} else*/
-	if (st == eStdHep) {

+	} else if (st == eStdHep) {

 		src = new StdHepEventSource(getFilename());
 	} else if (st == eGPS) {
 		src = new GPSEventSource();

@@ -261,10 +263,10 @@

 
 		// Reset event counter.
 		m_ngen = 0;

-	} /*else {

+	} else {

 		log() << LOG::fatal << "No event generator was setup." << LOG::done;
 		RunManager::instance()->abortRun(SlicApplication::GENERATOR_NOT_SETUP);

-	}*/

+	}

 }
 
 void EventSourceManager::setupEventSource() {

@@ -286,11 +288,8 @@

 	bool eof = false;
 	if (isFileSource()) {
 		EventSourceWithInputFile* src = dynamic_cast<EventSourceWithInputFile*>(m_currentEventSource);

-

 		if (src) {
 			eof = src->isEOF();

-		} else {
-			G4Exception("", "", JustWarning, "Cast to EventSourceWithInputFile failed.");

 		}
 	}
 	return eof;

@@ -338,7 +337,7 @@

 	return gun;
 }
 

-void EventSourceManager::GeneratePrimaryVertex(G4Event* evt) {

+void EventSourceManager::generateNextEvent(G4Event* evt) {

 	generate(evt);
 }
 

@@ -346,10 +345,7 @@

 
 	m_lorentzTransformationAngle = lorentzTransformationAngle;
 #ifdef SLIC_LOG

-	log() << LOG::always
-	<< "Set Lorentz transformation angle to "
-	<< G4BestUnit(m_lorentzTransformationAngle, "Angle")
-	<< LOG::done;

+	log() << LOG::always << "Set Lorentz transformation angle to " << G4BestUnit(m_lorentzTransformationAngle, "Angle") << LOG::done;

 #endif
 
 }

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioEventSource.cc	                        (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioEventSource.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -0,0 +1,79 @@

+// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioEventSource.cc,v 1.9 2012-11-27 19:32:19 jeremy Exp $
+#include "LcioEventSource.hh"
+
+// slic
+//#include "LcioMcpFilter.hh"
+
+namespace slic {
+
+LcioEventSource::LcioEventSource(const std::string& fname) :
+		EventSourceWithInputFile("LcioEventSource", fname), m_eventGenerator(0) {
+}
+
+LcioEventSource::~LcioEventSource() {
+    delete m_eventGenerator;
+}
+
+// open the current file
+void LcioEventSource::open() {
+    m_eventGenerator = new MCParticleGenerator(m_filename);
+    m_fileIsOpen = true;
+}
+
+// close the current file
+void LcioEventSource::close() {
+    delete m_eventGenerator;
+    m_eventGenerator = 0;
+    m_fileIsOpen = false;
+}
+
+// read the next event
+void LcioEventSource::readNextEvent() {
+    m_eventGenerator->readNextEvent();
+    if (m_eventGenerator->getCurrentEvent() == 0) {
+        m_eof = true;
+    }
+}
+
+void LcioEventSource::dumpCurrentEvent() {
+	//log() << LOG::debug << "LcioEventSource - Initial MCParticle Collection" << LOG::done;
+	//LcioMcpManager::instance()->printMcpCollection("initial");
+}
+
+void LcioEventSource::generate(G4Event* anEvent) {
+    m_eventGenerator->generateEvent(anEvent);
+}
+
+void LcioEventSource::beginRun(const G4Run* aRun) {
+	// do superclass action
+	EventSourceWithInputFile::beginRun(aRun);
+}
+
+void LcioEventSource::beginEvent(const G4Event* anEvent) {
+	// read an event
+	EventSourceWithInputFile::beginEvent(anEvent);
+
+	// set initial collection to filtered coll w/ gen parts only
+	//LcioMcpManager::instance()->setInitialMcpCollection(
+	//		m_filter->filterSimParticles(getCurrentMcpLCCollection()));
+}
+
+/*
+EVENT::LCCollection* LcioEventSource::getCurrentMcpLCCollection() {
+	EVENT::LCCollection* mcpColl = 0;
+
+	if (m_event) {
+
+		mcpColl = LcioMcpManager::instance()->findMcpCollection(m_event);
+
+		if (mcpColl == 0) {
+			G4Exception("", "", FatalException, "No Mcp collection found in LCIO input event.");
+		}
+	} else {
+		log() << LOG::error << "LCEvent is null" << LOG::done;
+	}
+	return mcpColl;
+}
+*/
+
+}

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/MCParticleGenerator.cc	                        (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/MCParticleGenerator.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -0,0 +1,57 @@

+#include "MCParticleGenerator.hh"
+
+// SLIC
+#include "MCParticleManager.hh"
+
+// LCIO
+#include "IOIMPL/LCFactory.h"
+#include "IMPL/MCParticleImpl.h"
+#include "EVENT/LCEvent.h"
+#include "EVENT/LCIO.h"
+
+// Geant4
+#include "G4SystemOfUnits.hh"
+
+using EVENT::MCParticle;
+using EVENT::LCCollection;
+using EVENT::LCEvent;
+using IMPL::MCParticleImpl;
+using IOIMPL::LCFactory;
+
+namespace slic {
+
+MCParticleGenerator::MCParticleGenerator(G4String eventFile)
+    : _lcevent(0) {
+
+    _particleManager = MCParticleManager::instance();
+
+    /**
+     * Setup the reader.
+     */
+    LCFactory* factory = LCFactory::getInstance();
+    _reader = factory->createLCReader();
+    _reader->open(eventFile);
+}
+
+MCParticleGenerator::~MCParticleGenerator() {
+    _reader->close();
+}
+
+void MCParticleGenerator::readNextEvent() {
+    _lcevent = _reader->readNextEvent(EVENT::LCIO::UPDATE);
+}
+
+LCEvent* MCParticleGenerator::getCurrentEvent() {
+    return _lcevent;
+}
+
+void MCParticleGenerator::generateEvent(G4Event* event) {
+
+    /* Get collection of MCParticles. */
+    LCCollectionVec* particles = (LCCollectionVec*)_lcevent->takeCollection("MCParticle");
+
+    /* Generate the Geant4 event using the MCParticleManager. */
+    _particleManager->generateEvent(particles, event);
+}
+
+} /* namespace slic */

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/MCParticleManager.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/MCParticleManager.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -1,8 +1,251 @@

 #include "MCParticleManager.hh"
 

+// SLIC
+#include "EventSourceManager.hh"
+
+// Geant4
+#include "G4SystemOfUnits.hh"
+#include "Randomize.hh"
+

 namespace slic {
 

-MCParticleGeant4Map MCParticleManager::_particleMap;
-VisitedMCParticles MCParticleManager::_visitedParticles;

+void MCParticleManager::generateEvent(LCCollectionVec* particles, G4Event* event) {

+    /* Apply the Lorentz Transformation to input particles. */
+    applyLorentzTransformation(particles, EventSourceManager::instance()->getLorentzTransformationAngle());
+
+    /* Apply Z smearing to input particles. */
+    applyZSmearing(particles, EventSourceManager::instance()->getZSmearing());
+
+    /*
+     * Process particles and add them to the primary vertex of the event.
+     */
+    G4ThreeVector particlePosition;
+    G4double particleTime;
+    G4PrimaryVertex* vertex = new G4PrimaryVertex(particlePosition, particleTime);
+    for(size_t i=0; i < particles->size(); i++) {
+        MCParticle* mcp = (MCParticle*)(*particles)[i];
+        if (mcp->getParents().size()==0) {
+            //std::cout << "construct tree of relevant G4PrimaryParticles for initial MCParticle " << i << " " << std::flush;
+            std::set<G4PrimaryParticle*> primaries = getRelevantParticles(mcp);
+            //std::cout << "... finished" << std::endl;
+            std::set<G4PrimaryParticle*>::iterator primaryIt;
+            for (primaryIt = primaries.begin(); primaryIt != primaries.end(); primaryIt++) {
+                vertex->SetPrimary(*primaryIt);
+                //std::cout << "G4PrimaryParticle ("<< (*primaryIt)->GetPDGcode() << ") added to G4PrimaryVertex." << std::endl;
+            }
+        }
+    }
+
+    /* Add the primary vertex to Geant4. */
+    event->AddPrimaryVertex(vertex);
+
+    /* Register the collection of MCParticles that was created. */
+    setMCParticleCollection(particles);
+}
+
+
+std::set<G4PrimaryParticle*> MCParticleManager::getRelevantParticles(MCParticle* p) {
+
+    //log each particle which has been called, to avoid double counting and increase efficiency
+    _visitedParticles.insert(p);
+    PrimaryParticleMap::iterator mcpIT;
+    std::set<G4PrimaryParticle*> relevantParticlesSet; //holds all relevant decay particles of p
+
+    /*
+     * Case #1:
+     *
+     * If p is a stable particle: search for it in LCIO2Geant4Map
+     *
+     * If exists in map: get the G4PrimaryParticle version of p from LCIO2Geant4Map
+     * Else: create G4PrimaryParticle version write it to LCIO2Geant4Map
+     *
+     * Finally: insert the G4PrimaryParticle verion of p into the relevant particle list
+     * and return the list.
+     */
+    if (p->getGeneratorStatus() == 1) {
+        G4PrimaryParticle* g4p;
+        mcpIT = _primaryMap.find(p);
+        if (mcpIT != _primaryMap.end()) {
+            g4p = mcpIT->second;
+        } else {
+            G4int IDHEP = p->getPDG();
+            G4double PHEP1 = p->getMomentum()[0];
+            G4double PHEP2 = p->getMomentum()[1];
+            G4double PHEP3 = p->getMomentum()[2];
+            G4double PHEP5 = p->getMass();
+            // create G4PrimaryParticle object
+            g4p = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
+            g4p->SetMass(PHEP5 * GeV);
+            _primaryMap[p] = g4p;
+            //std::cout << "*" << std::flush;
+        }
+        //std::cout << g4p->GetPDGcode() << std::flush;
+        relevantParticlesSet.insert(g4p);
+        return relevantParticlesSet;
+    }
+
+    /*
+     * CASE #2:
+     *
+     * If p is not stable: get first decay daughter and calculate the proper time of p.
+     *
+     * If proper time is not zero: particle is "relevant", since it carries vertex information.
+     *
+     * If p is already in LCIO2Geant4Map: take it.
+     * Else: create G4PrimaryParticle version of p and write it to LCIO2Geant4Map
+     *
+     * Collect all relevant particles of all daughters and setup "relevant mother-
+     * daughter-relations" between relevant decay particles and G4PrimaryParticle version of p
+     *
+     * Finally: insert only the G4PrimaryParticle version of p to the relevant particle list and return this "list".
+     * The added particle has now the correct pre-assigned decay products and (hopefully) the right lifetime.
+     */
+    else if (p->getDaughters().size() > 0) {  //fg: require that there is at least one daughter - otherwise forget the particle
+        // calculate proper time
+        MCParticle* dp = p->getDaughters()[0];
+
+        double proper_time = fabs((dp->getTime() - p->getTime()) * p->getMass()) / p->getEnergy();
+        //  fix by S.Morozov for real != 0
+        double aPrecision = dp->getTime() * p->getMass() / p->getEnergy();
+        double bPrecision = p->getTime() * p->getMass() / p->getEnergy();
+
+        double proper_time_Precision = pow(10, -DBL_DIG)*fmax(fabs(aPrecision),fabs(bPrecision));
+
+        bool isProperTimeZero = (proper_time <= proper_time_Precision);
+
+        // -- remove original --- if (proper_time != 0) {
+        if (isProperTimeZero == false) {
+
+            G4PrimaryParticle* g4p;
+            mcpIT = _primaryMap.find(p);
+            if (mcpIT != _primaryMap.end()) {
+                g4p = mcpIT->second;
+            } else {
+                G4int IDHEP = p->getPDG();
+                G4double PHEP1 = p->getMomentum()[0];
+                G4double PHEP2 = p->getMomentum()[1];
+                G4double PHEP3 = p->getMomentum()[2];
+                G4double PHEP5 = p->getMass();
+                // create G4PrimaryParticle object
+                g4p = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
+                g4p->SetMass(PHEP5 * GeV);
+                g4p->SetProperTime(proper_time * ns);
+                _primaryMap[p] = g4p;
+                std::set<G4PrimaryParticle*> vec3;
+                for (size_t i = 0; i < p->getDaughters().size(); i++) {
+                    if (_visitedParticles.count(p->getDaughters()[i]) == 0) {
+                        std::set<G4PrimaryParticle*> vec2 = getRelevantParticles(p->getDaughters()[i]);
+                        std::set<G4PrimaryParticle*>::iterator setit;
+                        for (setit = vec2.begin(); setit != vec2.end(); setit++) {
+                            vec3.insert(*setit);
+                        }
+                    }
+                }
+                std::set<G4PrimaryParticle*>::iterator setit;
+                for (setit = vec3.begin(); setit != vec3.end(); setit++) {
+                    g4p->SetDaughter(*setit);
+                }
+                //std::cout << "*" << std::flush;
+            }
+            //std::cout << g4p->GetPDGcode() << std::flush;
+            relevantParticlesSet.insert(g4p);
+            return relevantParticlesSet;
+        }
+        /* CASE #3:
+         *
+         * If p is not stable AND has zero lifetime: forget about p and retrieve all relevant
+         * decay particles of all daughters of p.
+         *
+         * In this case this step of recursion is just there for collecting all relevant decay
+         * products of daughters and returning them.
+         */
+        else {
+            for (size_t i = 0; i < p->getDaughters().size(); i++) {
+                if (_visitedParticles.count(p->getDaughters()[i]) == 0) {
+                    std::set<G4PrimaryParticle*> vec2 = getRelevantParticles(p->getDaughters()[i]);
+                    std::set<G4PrimaryParticle*>::iterator setit;
+                    for (setit = vec2.begin(); setit != vec2.end(); setit++) {
+                        relevantParticlesSet.insert(*setit);
+                    }
+                }
+            }
+            return relevantParticlesSet;
+        }
+    }
+    return relevantParticlesSet;
+}
+
+void MCParticleManager::applyLorentzTransformation(LCCollection* particles, const G4double alpha) {
+
+    if (alpha == 0)
+        return; // nothing to do
+
+#ifdef SLIC_LOG
+        log() << LOG::name << LOG::always << "Applying Lorentz Transformation angle " << alpha << " to MCParticles." << LOG::done;
+#endif
+
+    // parameters of the Lorentz transformation matrix
+    const G4double gamma = sqrt(1 + sqr(tan(alpha)));
+    const G4double betagamma = tan(alpha);
+
+    if (particles != 0) {
+
+        int nMCP = particles->getNumberOfElements();
+
+        for (int i = 0; i < nMCP; ++i) {
+
+            IMPL::MCParticleImpl* mcp = dynamic_cast<IMPL::MCParticleImpl*>(particles->getElementAt(i));
+
+            const double* p = mcp->getMomentum();
+
+            // before the transformation
+            double pPrime[3];
+
+            const G4double m = mcp->getMass();
+
+            // after the transformation (boost in x-direction)
+            pPrime[0] = betagamma * sqrt(sqr(p[0]) + sqr(p[1]) + sqr(p[2]) + sqr(m)) + gamma * p[0];
+
+            pPrime[1] = p[1];
+            pPrime[2] = p[2];
+
+            // py and pz remain the same, E changes implicitly with px
+            mcp->setMomentum(pPrime);
+        }
+    }
+}
+
+void MCParticleManager::applyZSmearing(LCCollection* particles, const G4double z) {
+
+    //G4cout << "zsmear" << G4endl;
+
+    //G4cout << "got zParam = " << zParam << G4endl;
+
+    if (z == 0)
+        return;
+
+    // Generate smeared Z position.
+    double zspread = (z == 0.0 ? 0.0 : G4RandGauss::shoot(0, z/mm));
+
+    // Apply Z smearing to MCParticle collection.
+    if (particles != 0) {
+        int nMCP = particles->getNumberOfElements();
+        //G4cout << "nMCP = " << nMCP << G4endl;
+        for (int i = 0; i < nMCP; ++i) {
+            IMPL::MCParticleImpl* mcp = dynamic_cast<IMPL::MCParticleImpl*>(particles->getElementAt(i));
+            if (z != 0) {
+                const double* v = mcp->getVertex();
+                double vPrime[3];
+                vPrime[0] = v[0];
+                vPrime[1] = v[1];
+                vPrime[2] = v[2] + zspread;
+                //G4cout << "MCP.preZSmear = " << mcp->getVertex()[2] << G4endl;
+                mcp->setVertex(vPrime);
+                //G4cout << "MCP.postZSmear = " << mcp->getVertex()[2] << G4endl;
+            }
+        }
+    }
+}
+

 };

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/PrimaryGeneratorAction.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/PrimaryGeneratorAction.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -1,13 +1,13 @@

 // $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/PrimaryGeneratorAction.cc,v 1.38 2013-06-26 01:57:02 jeremy Exp $
 
 // SLIC

-#include "PrimaryGeneratorAction.hh"

+#include "EventSourceManager.hh"

 #include "LcioManager.hh"

+#include "MCParticleManager.hh"
+#include "PrimaryGeneratorAction.hh"
+#include "TrackManager.hh"

 #include "RunManager.hh"
 #include "SlicApplication.hh"

-//#include "LcioMcpManager.hh"
-//#include "LcioMcpFilter.hh"
-#include "EventSourceManager.hh"

 
 // LCDD
 #include "lcdd/util/StringUtil.hh"

@@ -16,13 +16,11 @@

 #include "IMPL/MCParticleImpl.h"
 #include "EVENT/LCCollection.h"
 

-// Geant4
-#include "G4Types.hh"
-#include "G4StateManager.hh"
-

 namespace slic {

+

 PrimaryGeneratorAction::PrimaryGeneratorAction() :
 		Module("PrimaryGeneratorAction", false) {

+    _manager = EventSourceManager::instance();

 }
 
 PrimaryGeneratorAction::~PrimaryGeneratorAction() {

@@ -33,165 +31,31 @@

     /* Print begin of event message. */
     printBeginEventMessage(anEvent);
 

-    /* Get the event source manager. */
-    EventSourceManager* mgr = EventSourceManager::instance();

+    /* Reset the TrackManager's state from previous event. */
+    TrackManager::instance()->reset();

-    /* Begin event hook of event source manager. */
-    mgr->beginEvent(anEvent);

+    /* Reset the MCParticleManager's state from previous event. */
+    MCParticleManager::instance()->reset();

-    /* Generate primary vertices using the current event generator. */
-    mgr->GeneratePrimaryVertex(anEvent);

+    /* Check if run needs to be aborted before generating events. */
+    if (RunManager::instance()->isRunAborted()) {
+        G4RunManager::GetRunManager()->AbortRun();
+        log() << LOG::warning << LOG::name << LOG::sep << "Run was aborted.  Events will not be generated." << LOG::done;
+    } else {
+        /* Begin event hook of EventSourceManager. */
+        _manager->beginEvent(anEvent);

-    // DEBUG: stop program here!!!
-    //assert(false);

+        /* Generate primary vertices using the current event generator. */
+        _manager->generateNextEvent(anEvent);
+    }

-    /*
-	// Check if run needs to be aborted, in which case we immediately return without executing anything.
-	if (RunManager::instance()->isRunAborted()) {
-		G4RunManager::GetRunManager()->AbortRun();
-		log() << LOG::warning << LOG::name << LOG::sep << "Run was already aborted.  Will not generate events." << LOG::done;
-		return;
-	}
-
-	// Apply Lorentz transformation to event, used for non-zero beam crossing angle.
-	applyLorentzTransformation(anEvent);
-
-	// Apply Z smearing to event vertices.
-	zsmear(anEvent);
-
-	// If source has no more events then abort the run.
-	if (mgr->isEOF()) {

+	/* If the event source hit the end of file, then abort the run. */
+	if (_manager->isEOF()) {

 		SlicApplication::instance()->setAborting(true);
 	}

-	*/

 }
 
 void PrimaryGeneratorAction::printBeginEventMessage(G4Event* anEvent) {

-	log() << LOG::okay << ">>>> BeginEvent <" << StringUtil::toString(anEvent->GetEventID()) << ">"
-			<< LOG::done;

+	log() << LOG::okay << ">>>> BeginEvent <" << StringUtil::toString(anEvent->GetEventID()) << ">" << LOG::done;

 }

-
-// TODO: The lorentz transform and z smearing functions below should be moved elsewhere!!!
-
-/*
-void PrimaryGeneratorAction::applyLorentzTransformation(G4Event *evt) {
-	const G4double alpha = EventSourceManager::instance()->getLorentzTransformationAngle();
-
-	if (alpha == 0)
-		return; // nothing to do
-
-//#ifdef SLIC_LOG
-//        log() << LOG::always << "applying Lorentz Transformation angle <" << alpha << LOG::done;
-//#endif
-
-	// parameters of the Lorentz transformation matrix
-	const G4double gamma = sqrt(1 + sqr(tan(alpha)));
-	const G4double betagamma = tan(alpha);
-
-	// scan through all vertices and all valid primary particles
-	for (int iVertex = 0; iVertex < evt->GetNumberOfPrimaryVertex(); iVertex++) {
-		G4PrimaryVertex *vertex = evt->GetPrimaryVertex(iVertex);
-		for (int iParticle = 0; iParticle < vertex->GetNumberOfParticle(); iParticle++) {
-			G4PrimaryParticle *particle = vertex->GetPrimary(iParticle);
-			// does the particle have a valid particle definition attached?
-			if (particle->GetG4code()) {
-				// before the transformation
-				const G4double px = particle->GetPx();
-				const G4double py = particle->GetPy();
-				const G4double pz = particle->GetPz();
-
-				//FG: take generator mass ( not necessarily euqal to PDG mass ) 
-				const G4double m = particle->GetMass();
-
-				// after the transformation (boost in x-direction)
-				const G4double pxPrime = betagamma * sqrt(sqr(px) + sqr(py) + sqr(pz) + sqr(m)) + gamma * px;
-
-				// py and pz remain the same, E changes implicitly with px
-				particle->SetMomentum(pxPrime, py, pz);
-			}
-		}
-		// the position of the vertex is not transformed here
-	}
-
-	EVENT::LCCollection* col = LcioMcpManager::instance()->getInitialMcpCollection();
-
-	if (col != 0) { // if particle gun is used no collection exists ...
-
-		int nMCP = col->getNumberOfElements();
-
-		for (int i = 0; i < nMCP; ++i) {
-
-			IMPL::MCParticleImpl* mcp = dynamic_cast<IMPL::MCParticleImpl*>(col->getElementAt(i));
-
-			const double* p = mcp->getMomentum();
-
-			// before the transformation
-			double pPrime[3];
-
-			const G4double m = mcp->getMass();
-
-			// after the transformation (boost in x-direction)
-			pPrime[0] = betagamma * sqrt(sqr(p[0]) + sqr(p[1]) + sqr(p[2]) + sqr(m)) + gamma * p[0];
-
-			pPrime[1] = p[1];
-			pPrime[2] = p[2];
-
-			// py and pz remain the same, E changes implicitly with px
-			mcp->setMomentum(pPrime);
-		}
-	}

 }

-*/
-
-/**
- * Apply a Z smear to the position of every G4PrimaryVertex in an event.
- */
-/*
-void PrimaryGeneratorAction::zsmear(G4Event* anEvent) {
-
-	//G4cout << "zsmear" << G4endl;
-
-	G4double zParam = EventSourceManager::instance()->getZSmearing();
-
-	//G4cout << "got zParam = " << zParam << G4endl;
-
-	if (zParam == 0)
-		return;
-
-	// Generate smeared Z position.
-	double zspread = (zParam == 0.0 ? 0.0 : G4RandGauss::shoot(0, zParam/mm));
-
-	//G4cout << "zspread = " << zspread << G4endl;
-
-	// Apply Z smearing to Geant4 event.
-	int nVertices = anEvent->GetNumberOfPrimaryVertex();
-	//G4cout << "nVertices = " << nVertices << G4endl;
-	for (int iVertex = 0; iVertex < nVertices; iVertex++) {
-		G4PrimaryVertex* vertex = anEvent->GetPrimaryVertex(iVertex);
-		//G4cout << "G4.preZSmear = " << vertex->GetZ0() << G4endl;
-		vertex->SetPosition(vertex->GetX0(), vertex->GetY0(), vertex->GetZ0() + zspread);
-		//G4cout << "G4.postZSmear = " << vertex->GetZ0() << G4endl;
-	}
-
-	// Apply Z smearing to MCParticle collection.
-	EVENT::LCCollection* mcpColl = LcioMcpManager::instance()->getInitialMcpCollection();
-	if (mcpColl != 0) {
-		int nMCP = mcpColl->getNumberOfElements();
-		//G4cout << "nMCP = " << nMCP << G4endl;
-		for (int i = 0; i < nMCP; ++i) {
-			IMPL::MCParticleImpl* mcp = dynamic_cast<IMPL::MCParticleImpl*>(mcpColl->getElementAt(i));
-			if (zParam != 0) {
-				const double* v = mcp->getVertex();
-				double vPrime[3];
-				vPrime[0] = v[0];
-				vPrime[1] = v[1];
-				vPrime[2] = v[2] + zspread;
-				//G4cout << "MCP.preZSmear = " << mcp->getVertex()[2] << G4endl;
-				mcp->setVertex(vPrime);
-				//G4cout << "MCP.postZSmear = " << mcp->getVertex()[2] << G4endl;
-			}
-		}
-	}
-}
-*/
-}

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepGenerator.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepGenerator.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -1,5 +1,8 @@

 #include "StdHepGenerator.hh"
 

+// SLIC
+#include "MCParticleManager.hh"
+

 // LCIO
 #include "IMPL/MCParticleImpl.h"
 

@@ -13,193 +16,30 @@

 
 namespace slic {
 

-StdHepGenerator::StdHepGenerator(G4String eventFile) {

+StdHepGenerator::StdHepGenerator(G4String eventFile)
+    : _particles(0) {

     _reader = new LCStdHepRdr(eventFile.data());

-    _particleManager = MCParticleManager::instance();
-    _visitedParticles = _particleManager->getVisitedParticles();

 }
 
 StdHepGenerator::~StdHepGenerator() {
     delete _reader;
 }
 

-void StdHepGenerator::generateEvent(G4Event* event) {

+void StdHepGenerator::readNextEvent() {
+    /* Read in the next StdHep event which is converted to an MCParticle collection. */
+    _particles = _reader->readEvent();
+}

-    /* Clear previous state before starting event generation. */
-    MCParticleManager::instance()->clearMaps();
-
-    /* Get collection of MCParticles from the LCStdHepRdr. */
-    LCCollectionVec* mcpVec = _reader->readEvent();
-    std::vector<MCParticle*> collection(mcpVec->getNumberOfElements());
-    G4int n = collection.size();
-

+void StdHepGenerator::generateEvent(G4Event* event) {

     /*

-     * Copy input MCParticles to a vector in order to work with them easily.

+     * Use the MCParticleManager utility method to create the Geant4 event using the current
+     * collection of MCParticles.

      */

-    for(G4int i=0; i<n; i++) {
-        MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>(mcpVec->getElementAt(i));
-        collection[i] = mcp;
-    }
-
-    /*
-     * Process particles and add them to the primary vertex of the event.
-     */
-    G4ThreeVector particlePosition; // FIXME: add Z smear
-    G4double particleTime;
-    G4PrimaryVertex* vertex = new G4PrimaryVertex(particlePosition, particleTime);
-    for(size_t i=0; i < collection.size(); i++) {
-        MCParticle* mcp = collection[i];
-        if (mcp->getParents().size()==0) {
-            //std::cout << "construct tree of relevant G4PrimaryParticles for initial MCParticle " << i << " " << std::flush;
-            std::set<G4PrimaryParticle*> primaries = getRelevantParticles(mcp);
-            //std::cout << "... finished" << std::endl;
-            std::set<G4PrimaryParticle*>::iterator primaryIt;
-            for (primaryIt = primaries.begin(); primaryIt != primaries.end(); primaryIt++) {
-                vertex->SetPrimary(*primaryIt);
-                //std::cout << "G4PrimaryParticle ("<< (*primaryIt)->GetPDGcode() << ") added to G4PrimaryVertex." << std::endl;
-            }
-        }
-    }
-
-    /* Add the primary vertex to Geant4. */
-    event->AddPrimaryVertex(vertex);
-
-    /* Register the collection of MCParticles that was created. */
-    _particleManager->setMCParticleCollection(mcpVec);

+    MCParticleManager::instance()->generateEvent(_particles, event);

 }
 

-std::set<G4PrimaryParticle*> StdHepGenerator::getRelevantParticles(MCParticle* p) {
-
-    //log each particle which has been called, to avoid double counting and increase efficiency
-    _visitedParticles->insert(p);
-    MCParticleGeant4Map::iterator mcpIT;
-    std::set<G4PrimaryParticle*> relevantParticlesSet; //holds all relevant decay particles of p
-
-    MCParticleGeant4Map* particleMap = _particleManager->getMCParticleToGeant4Map();
-
-    /*
-     * Case #1:
-     *
-     * If p is a stable particle: search for it in LCIO2Geant4Map
-     *
-     * If exists in map: get the G4PrimaryParticle version of p from LCIO2Geant4Map
-     * Else: create G4PrimaryParticle version write it to LCIO2Geant4Map
-     *
-     * Finally: insert the G4PrimaryParticle verion of p into the relevant particle list
-     * and return the list.
-     */
-    if (p->getGeneratorStatus() == 1) {
-        G4PrimaryParticle* g4p;
-        mcpIT = particleMap->find(p);
-        if (mcpIT != particleMap->end()) {
-            g4p = mcpIT->second;
-        } else {
-            G4int IDHEP = p->getPDG();
-            G4double PHEP1 = p->getMomentum()[0];
-            G4double PHEP2 = p->getMomentum()[1];
-            G4double PHEP3 = p->getMomentum()[2];
-            G4double PHEP5 = p->getMass();
-            // create G4PrimaryParticle object
-            g4p = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
-            g4p->SetMass(PHEP5 * GeV);
-            (*particleMap)[p] = g4p;
-            //std::cout << "*" << std::flush;
-        }
-        //std::cout << g4p->GetPDGcode() << std::flush;
-        relevantParticlesSet.insert(g4p);
-        return relevantParticlesSet;
-    }
-
-    /*
-     * CASE #2:
-     *
-     * If p is not stable: get first decay daughter and calculate the proper time of p.
-     *
-     * If proper time is not zero: particle is "relevant", since it carries vertex information.
-     *
-     * If p is already in LCIO2Geant4Map: take it.
-     * Else: create G4PrimaryParticle version of p and write it to LCIO2Geant4Map
-     *
-     * Collect all relevant particles of all daughters and setup "relevant mother-
-     * daughter-relations" between relevant decay particles and G4PrimaryParticle version of p
-     *
-     * Finally: insert only the G4PrimaryParticle version of p to the relevant particle list and return this "list".
-     * The added particle has now the correct pre-assigned decay products and (hopefully) the right lifetime.
-     */
-    else if (p->getDaughters().size() > 0) {  //fg: require that there is at least one daughter - otherwise forget the particle
-        // calculate proper time
-        MCParticle* dp = p->getDaughters()[0];
-
-        double proper_time = fabs((dp->getTime() - p->getTime()) * p->getMass()) / p->getEnergy();
-        //  fix by S.Morozov for real != 0
-        double aPrecision = dp->getTime() * p->getMass() / p->getEnergy();
-        double bPrecision = p->getTime() * p->getMass() / p->getEnergy();
-
-        double proper_time_Precision = pow(10, -DBL_DIG)*fmax(fabs(aPrecision),fabs(bPrecision));
-
-        bool isProperTimeZero = (proper_time <= proper_time_Precision);
-
-        // -- remove original --- if (proper_time != 0) {
-        if (isProperTimeZero == false) {
-
-            G4PrimaryParticle* g4p;
-            mcpIT = particleMap->find(p);
-            if (mcpIT != particleMap->end()) {
-                g4p = mcpIT->second;
-            } else {
-                G4int IDHEP = p->getPDG();
-                G4double PHEP1 = p->getMomentum()[0];
-                G4double PHEP2 = p->getMomentum()[1];
-                G4double PHEP3 = p->getMomentum()[2];
-                G4double PHEP5 = p->getMass();
-                // create G4PrimaryParticle object
-                g4p = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
-                g4p->SetMass(PHEP5 * GeV);
-                g4p->SetProperTime(proper_time * ns);
-                (*particleMap)[p] = g4p;
-                std::set<G4PrimaryParticle*> vec3;
-                for (size_t i = 0; i < p->getDaughters().size(); i++) {
-                    if (_visitedParticles->count(p->getDaughters()[i]) == 0) {
-                        std::set<G4PrimaryParticle*> vec2 = getRelevantParticles(p->getDaughters()[i]);
-                        std::set<G4PrimaryParticle*>::iterator setit;
-                        for (setit = vec2.begin(); setit != vec2.end(); setit++) {
-                            vec3.insert(*setit);
-                        }
-                    }
-                }
-                std::set<G4PrimaryParticle*>::iterator setit;
-                for (setit = vec3.begin(); setit != vec3.end(); setit++) {
-                    g4p->SetDaughter(*setit);
-                }
-                //std::cout << "*" << std::flush;
-            }
-            //std::cout << g4p->GetPDGcode() << std::flush;
-            relevantParticlesSet.insert(g4p);
-            return relevantParticlesSet;
-        }
-        /* CASE #3:
-         *
-         * If p is not stable AND has zero lifetime: forget about p and retrieve all relevant
-         * decay particles of all daughters of p.
-         *
-         * In this case this step of recursion is just there for collecting all relevant decay
-         * products of daughters and returning them.
-         */
-        else {
-            for (size_t i = 0; i < p->getDaughters().size(); i++) {
-                if (_visitedParticles->count(p->getDaughters()[i]) == 0) {
-                    std::set<G4PrimaryParticle*> vec2 = getRelevantParticles(p->getDaughters()[i]);
-                    std::set<G4PrimaryParticle*>::iterator setit;
-                    for (setit = vec2.begin(); setit != vec2.end(); setit++) {
-                        relevantParticlesSet.insert(*setit);
-                    }
-                }
-            }
-            return relevantParticlesSet;
-        }
-    }
-    return relevantParticlesSet;

+LCCollectionVec* StdHepGenerator::getCurrentParticleCollection() {
+    return _particles;

 }
 

-

 } /* namespace slic */

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackManager.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackManager.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -18,10 +18,10 @@

 
 void TrackManager::saveTrackSummaries(const G4Event* anEvent, LCEvent* lcEvent) {
 

-    MCParticleGeant4MapIterator it;

+    PrimaryParticleMapIterator it;

 
     /* Get the map of MCParticles to G4PrimaryParticles.  The latter has the track IDs. */

-    MCParticleGeant4Map* particleMap = MCParticleManager::instance()->getMCParticleToGeant4Map();

+    PrimaryParticleMap* particleMap = MCParticleManager::instance()->getPrimaryParticleMap();

 
     /* Copy MCParticle information from generated primary particles into the corresponding TrackSummary objects. */
     for (it = particleMap->begin(); it != particleMap->end(); it++) {

--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackSummary.cc	2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackSummary.cc	2014-01-15 00:17:32 UTC (rev 2876)
@@ -193,10 +193,10 @@

 
     /* Set momentum. */
     float p[3];

-    G4ThreeVector theMomentum = getMomentum();
-    p[0] = theMomentum(0) / GeV;
-    p[1] = theMomentum(1) / GeV;
-    p[2] = theMomentum(2) / GeV;

+    G4ThreeVector momentum = getMomentum();
+    p[0] = momentum(0) / GeV;
+    p[1] = momentum(1) / GeV;
+    p[2] = momentum(2) / GeV;

     _mcparticle->setMomentum(p);
 
     /* Set vertex. */