4 added + 15 modified, total 19 files
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioManager.cc 2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioManager.cc 2014-01-15 00:17:32 UTC (rev 2876)
@@ -336,18 +336,19 @@
lcevt->setRunNumber(m_runHdr->getRunNumber());
lcevt->setDetectorName(m_runHdr->getDetectorName());
+ // TODO: Add this back!!!
// Set the event weight and idrup from the StdHep file.
- EventSourceManager* genMgr = EventSourceManager::instance();
- if (genMgr->getCurrentSourceType() == EventSourceManager::eStdHep) {
- lStdHep* stdhep = ((StdHepEventSource*) (genMgr->getCurrentSource()))->getStdHepReader();
- if (stdhep->isStdHepEv4()) {
- lcevt->parameters().setValue("_weight", (float) stdhep->eventweight());
- lcevt->parameters().setValue("idrup", (int) stdhep->idrup());
- }
- } else {
- lcevt->parameters().setValue("_weight", (float) 1.0);
- lcevt->parameters().setValue("idrup", (float) 0.0);
- }
+ //EventSourceManager* genMgr = EventSourceManager::instance();
+ //if (genMgr->getCurrentSourceType() == EventSourceManager::eStdHep) {
+ // lStdHep* stdhep = ((StdHepEventSource*) (genMgr->getCurrentSource()))->getStdHepReader();
+ // if (stdhep->isStdHepEv4()) {
+ // lcevt->parameters().setValue("_weight", (float) stdhep->eventweight());
+ // lcevt->parameters().setValue("idrup", (int) stdhep->idrup());
+ // }
+ //} else {
+ // lcevt->parameters().setValue("_weight", (float) 1.0);
+ // lcevt->parameters().setValue("idrup", (float) 0.0);
+ //}
// Write Geant4 version into event header.
lcevt->parameters().setValue("GEANT4_VERSION", Geant4VersionInfo::getVersion());
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- 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
--- 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
--- 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
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepEventSource.cc 2014-01-14 01:09:17 UTC (rev 2875)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepEventSource.cc 2014-01-15 00:17:32 UTC (rev 2876)
@@ -7,65 +7,75 @@
namespace slic {
StdHepEventSource::StdHepEventSource(const std::string& fname) :
EventSourceWithInputFile("StdHepEventSource", fname),
- m_eventGenerator(new StdHepGenerator(fname)) {
-
- // create loader with internal reader
- m_loader = new StdHepLoader();
+ m_eventGenerator(0) {
}
void StdHepEventSource::generate(G4Event* anEvent) {
+ //std::cout << "StdHepEventSource::generate" << std::endl;
+
m_eventGenerator->generateEvent(anEvent);
}
// open the current file
void StdHepEventSource::open() {
- m_loader->openStdHepFile(m_filename);
+
+ //std::cout << "StdHepEventSource::open" << std::endl;
+
+ /* Initialize the event generator using the current file name. */
+ m_eventGenerator = new StdHepGenerator(getFilename());
+
m_fileIsOpen = true;
}
// close the current file
void StdHepEventSource::close() {
- m_loader->closeStdHepFile();
+
+ //std::cout << "StdHepEventSource::close" << std::endl;
+
+ delete m_eventGenerator;
+
m_fileIsOpen = false;
}
// read the next event
void StdHepEventSource::readNextEvent() {
- m_loader->readNextEvent();
- if (m_loader->isEndOfInput()) {
- m_eof = true;
+ //std::cout << "StdHepEventSource::readNextEvent" << std::endl;
+
+ m_eventGenerator->readNextEvent();
+ if (m_eventGenerator->getCurrentParticleCollection() == 0) {
+ m_eof = true;
}
}
void StdHepEventSource::dumpCurrentEvent() {
- log() << LOG::okay << "Dumping StdHep event info ..." << LOG::done;
- lStdHep* rdr = m_loader->getStdHepReader();
+ // TODO: Reimplement using MCParticle dump.
- if (rdr) {
- rdr->printEventTable();
- } else {
- log() << LOG::error << "No current StdHep reader." << LOG::done;
- }
-
- log() << LOG::okay << LOG::endl;
+ //log() << LOG::okay << "Dumping StdHep event info ..." << LOG::done;
+ //lStdHep* rdr = m_loader->getStdHepReader();
+ //if (rdr) {
+ // rdr->printEventTable();
+ //} else {
+ // log() << LOG::error << "No current StdHep reader." << LOG::done;
+ //}
+ //log() << LOG::okay << LOG::endl;
}
void StdHepEventSource::beginRun(const G4Run* aRun) {
+
+ //std::cout << "StdHepEventSource::beginRun" << std::endl;
+
// do superclass setup
EventSourceWithInputFile::beginRun(aRun);
-
- // setup convertor
- assert(m_loader->getStdHepReader());
}
void StdHepEventSource::beginEvent(const G4Event* anEvent) {
+
+ //std::cout << "StdHepEventSource::beginEvent" << std::endl;
+
// read an event
EventSourceWithInputFile::beginEvent(anEvent);
}
-lStdHep* StdHepEventSource::getStdHepReader() {
- return m_loader->getStdHepReader();
}
-}
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- 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
--- 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
--- 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. */
SVNspam 0.1